[nikwin]

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

[eugene]

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).

[nikwin]

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.

[eugene]

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.

[nikwin]

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!

[eugene]

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.

[nikwin]

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?

[eugene]

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.

[nikwin]

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!

[eugene]

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.

[nikwin]

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?