# Problem with trimmed cells when calculating \nabla \cdot \nabla \cdot (U^TU)

 Register Blogs Members List Search Today's Posts Mark Forums Read

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 Star-CD, 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::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
Attached Images
 divdivCellC2.jpg (73.8 KB, 35 views)

 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

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.
Attached Images
 divCellC2_X.jpg (53.9 KB, 19 views)

 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 non-uniform 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 rms-value 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 Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post hsieh Open Source Meshers: Gmsh, Netgen, CGNS, ... 32 September 13, 2011 05:50 liu OpenFOAM Running, Solving & CFD 6 December 30, 2005 18:27 kris CD-adapco 2 August 3, 2005 00:32

All times are GMT -4. The time now is 23:23.