CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

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

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

Reply
 
LinkBack Thread Tools Display Modes
Old   August 3, 2009, 03:59
Default 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
nikwin is on a distinguished road
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::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
Attached Images
File Type: jpg divdivCellC2.jpg (73.8 KB, 35 views)
nikwin is offline   Reply With Quote

Old   August 4, 2009, 06:11
Default
  #2
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12
eugene is on a distinguished road
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).
eugene is offline   Reply With Quote

Old   August 4, 2009, 07:22
Default
  #3
Member
 
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8
nikwin is on a distinguished road
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.
Attached Images
File Type: jpg divCellC2_X.jpg (53.9 KB, 19 views)
nikwin is offline   Reply With Quote

Old   August 4, 2009, 08:22
Default
  #4
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12
eugene is on a distinguished road
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.
eugene is offline   Reply With Quote

Old   August 6, 2009, 07:35
Default
  #5
Member
 
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8
nikwin is on a distinguished road
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!
nikwin is offline   Reply With Quote

Old   August 6, 2009, 07:43
Default
  #6
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12
eugene is on a distinguished road
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.
eugene is offline   Reply With Quote

Old   August 10, 2009, 04:22
Default
  #7
Member
 
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8
nikwin is on a distinguished road
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?
nikwin is offline   Reply With Quote

Old   August 10, 2009, 08:38
Default
  #8
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12
eugene is on a distinguished road
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.
eugene is offline   Reply With Quote

Old   August 12, 2009, 02:54
Default
  #9
Member
 
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8
nikwin is on a distinguished road
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!
nikwin is offline   Reply With Quote

Old   August 12, 2009, 06:23
Default
  #10
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12
eugene is on a distinguished road
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.
eugene is offline   Reply With Quote

Old   August 14, 2009, 04:19
Default
  #11
Member
 
Niklas Winkler
Join Date: Mar 2009
Location: Stockholm, Stockholm, Sweden
Posts: 73
Rep Power: 8
nikwin is on a distinguished road
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?
nikwin is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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 CD-adapco 2 August 3, 2005 00:32


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