
[Sponsors] 
Problem with trimmed cells when calculating \nabla \cdot \nabla \cdot (U^TU) 

LinkBack  Thread Tools  Display Modes 
August 3, 2009, 03:59 
Problem with trimmed cells when calculating \nabla \cdot \nabla \cdot (U^TU)

#1 
Member
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8 
Dear all,
I would like to calculate the divergence of the convective term in the momentum equation which is as above, \nabla \cdot \nabla \cdot (U^T U) on a mesh converted from StarCD, which includes trimmed cells. As a test I have calculated \nabla \cdot \nabla \cdot (X^T X), where X here denotes (x,y,z), which analytically should be 12 on a 3D mesh as follows, volVectorField cellC(mesh.C()); surfaceScalarField cellCsurf(linearInterpolate(cellC) & mesh.Sf()); volScalarField divdivCellC2 ( IOobject ( "divdivCellC2", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::div(fvc::div(cellCsurf,cellC)) ); The figure below shows results for a slice of the first bend of the double bended pipe. Any suggestions of what could be wrong or how to solve the problem? Thanks, /NW 

August 4, 2009, 06:11 

#2 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
hmm.
Since most of your cells do have a value of 12, I would guess you are getting it almost right. One thing that immediately springs to mind is that most people do not use skewness correction in their interpolation. Thus linearInterpolate(cellC) will not provide a faceCentred value for distorted elements. It would be interesting to see what the results would be if you replace the interpolation with mesh.Cf() (face centres). 

August 4, 2009, 07:22 

#3 
Member
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8 
Thanks for your reply Eugene!
Unfortunately, surfaceScalarField cellCsurf(mesh.Cf() & mesh.Sf()); didn't change anything, in fact it gives exactly the same result!? To give some more info, when I'm calculating div(cellC) the result is 3 which is as expected. But when calculating div(cellCsurf,cellC) some distortion can already be seen, see figure. 

August 4, 2009, 08:22 

#4 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
I think you can focus on finding the issue in div(cellCsurf,cellC). The fact that div(cellC) produces 3 means, div(div(cellCsurf, cellC)) should also be fine provided div(cellCsurf, cellC) is fine.
The code for this stuff can all be found in src/finiteVolume/lnInclude/... fvcDiv.* convectionSheme* gaussConvectionScheme* It is really really simple mathematically, so I am rather curious to find out what the issue is. If I were you I would focus on one problem cell and walk through the code step by step using a spreadsheet/debugger or similar. 

August 6, 2009, 07:35 

#5 
Member
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8 
Finally I got the debug version to work!
But, it's showing that fvcDiv.C is called and then a lot of underlying files which I believe are not directly related to the calculation. However, I can't find where the actual calculations are performed. Could you please briefly explain the steps from when the software calls fvc::div(...) or give me some other hints on how to proceed, please. Thanks! 

August 6, 2009, 07:43 

#6 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
My previous post already gives the answer: convectionScheme* and gaussConvectionScheme* also found in src/finiteVolume/lnInclude/
The only outstanding bits are the interpolation schemes, but since you are probably using linear interpolation its not really necessary to look at the code for these. If however, you find you do need to look at the interpolation code, let me know and I will provide the details for the interpolation scheme of you choice. 

August 10, 2009, 04:22 

#7 
Member
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8 
I've been trying to solve the issue but without any real progress since the problem is beyond my programming skills.
Is it possible to get some more help to solve the problem? And, is this an OpenFoam bug which should be reported? 

August 10, 2009, 08:38 

#8 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
I doubt very much it is a "bug" in the normal sense of the word. It is obviously a result of the numerical schemes being used. A drop in accuracy on a nonuniform mesh is not surprising.
Unfortunately I don't have the time to look at this with any degree of rigour at the moment. I will see what I can do simply because my curiosity is piqued, but the chances of a result in the near future would have to be considered slim. If you have any specific questions relating to the code, I would be happy to answer those. 

August 12, 2009, 02:54 

#9 
Member
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8 
Eugene, as you commented in the beginning and from another input it certainly looks like it is the interpolation scheme that causes me trouble due to that the problematic cells are skewed. Instead of using linearInterpolate to obtain the values at the face centers I've tried fvc:interpolate together with the upwind and QUICK schemes which of course give me different results but are still not correct. What I would like to try is a scheme which corrects for skewness as e.g. skewLinear according to table 4.6 in UG, but I can not get it to work. Do you know the syntax or any other possibility to include skewness correction into my interpolation?
Thanks! 

August 12, 2009, 06:23 

#10 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
The syntax for a skewCorrected scheme would be something like:
div(phiX,X) Gauss skewCorrected linear; At least this is how it works in 1.5. 

August 14, 2009, 04:19 

#11 
Member
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8 
Together with skewCorrection the performance of the calculation becomes considerably improved. However, on the StarCD mesh when calculating div(div(phiX,X)) the errors are still too large, getting results from 10 to 34 where it is expected to be 12. The rmsvalue when calculating div(phiX,X) is approximately 10^4 compared to 10^3 without the correction.
It was also tried on the 2D pitzDaily case with mesh refinement. The relative error decreases when the mesh is refined but the maximum error is more or less constant. Is there a way to further improve the skewCorrection? 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Import netgen mesh to OpenFOAM  hsieh  Open Source Meshers: Gmsh, Netgen, CGNS, ...  32  September 13, 2011 05:50 
How to update polyPatchbs localPoints  liu  OpenFOAM Running, Solving & CFD  6  December 30, 2005 18:27 
physical boundary error!!  kris  CDadapco  2  August 3, 2005 00:32 