View Issue Details

IDProjectCategoryView StatusLast Update
0001548OpenFOAM[All Projects] Bugpublic2015-06-18 15:14
Reporteruser528Assigned Tohenry 
Status resolvedResolutionfixed 
PlatformLinuxOSSUSE Linux Enterprise ServerOS Version11
Product Version 
Fixed in Version 
Summary0001548: rhoCentralFlow density buildup near outlet
DescriptionIn rhoCentralFlow, if I set at some boundary patch rather low fixedValue pressure and zeroGradient velocity, the density near that patch starts to accumulate instead of freely leaving the computational domain.
Steps To ReproduceUnzip the attached archive. It contains two identical (different are only the results) OpenFOAM cases. The "beforeCorrection" case has the results produced by 'standard' rhoCentralFoam, the "afterCorrection" has the results produced by corrected rhoCentralFoam (see below).

Each test has a nonuniform initial distribution of pressure, with two steps (jumps) inside the volume, and also it has a low fixedValue pressure at the "inlet" patch (sorry for the misleading patch name), thus creating one more jump similar to the internal ones.

Run rhoCentralFoam on the "beforeCorrection" test case. Observe that in the cell layer adjacent to the "inlet" patch, the density becomes rather high, while I expect the behavior near the "inlet" patch be similar to the behavior near the internal jumps.
Additional InformationI am not sure that this is indeed a bug in OpenFOAM, maybe I am missing some special boundary condition class to be used in such a situation, but nevertheless my choice of BCs seems rather sensible, and the results do not seem so.

My understanding of the reasons is as follows. The key idea of rhoCentralFlow algorithm is to have on each face two pairs of values for density (and other variables): one value is obtained from owner cell, and the other from neighbor cell. This works perfectly for internal faces, but not for the boundary, because the interpolate() function always uses the boundary field values of density to set the boundary field values for both rho_pos and rho_neg. Therefore, rho_pos is always equal to rho_neg on the boundary, which produces wrong flux across the boundary, etc. I can suppose that it has more sense to take the boundaryInternalField of density for rho_pos (similarly for other interpolate variables).

I have done corresponding corrections, see the diff in the attached archive. I am not completely sure that this is physically correct, but there is no density buildup, and the overall density and velocity profiles near the "inlet" patch look very similar to the profiles near the internal jumps.

Note also that rhoPimpleFoam (after the required corrections to fvSchemes and fvSolution) produces results more-or-less in agreement with the corrected version, though with greater oscillations.
TagsNo tags attached.



2015-02-26 11:20 (631,150 bytes)


2015-02-26 11:39

manager   ~0003918

Thanks for the cases, analysis and patch. I agree with your conclusions: the boundary conditions for the positive and negative fluxes should indeed be evaluated differently. Ideally this would be achieved by implementing a special type which includes both the positive and negative fluxes and implements appropriate boundary treatment. It is quite a lot of coding for a modest gain in code organisation but would probably be worth it in the longer term as the "central" codes would be simpler and easier to maintain. In the short-term a correction like the one you propose would be OK, I will put in a slightly more efficient version later today.


2015-02-26 17:12

manager   ~0003919

Resolved by commit 11ddd071091adc12602d3e56851a1e249ace7310

I have wrapped the interpolation and BC correction into a single templated function with some argument juggling to simplify use.

Let me know if this resolves the issues you have and does not introduce any new ones.


2015-02-27 09:28


Yes, your fix seems to work for me and I see no new issues.


2015-03-13 07:28


I found a funny bug in the fix. In case the mesh has no internal faces (for example, on a one cell mesh), the check
if (dir[0] > 0)
obviously fails.

I have replaces the code with
    forAll(sf.boundaryField(), patchi)
        if ((!sf.boundaryField()[patchi].coupled())
            && (dir.boundaryField()[patchi].size() > 0) && (dir.boundaryField()[patchi][0]>0))
            sf.boundaryField()[patchi] =
which should work even in such a case.

(Definitely, it is somewhat strange to run rhoCentralFoam on a one cell mesh, but I have some more advanced solver based on rhoCentralFoam code. It includes, in particular, chemistry, and I was in fact running the solver a one cell mesh to check for correctness of chemistry rates.)


2015-03-13 11:55

manager   ~0004107

Resolved by commit 13bd74c439ce17a97968a2a2437e0db20e47be43


2015-06-18 15:14

manager   ~0004966

I have reverted this correction to the BCs as it adversely affects fixed-value BCs and is formulated to fix an issue with an unphysical case. Further analysis of the handling of fixed pressure outlet conditions as the Mach number approaches 1 is required.

Issue History

Date Modified Username Field Change
2015-02-26 11:20 user528 New Issue
2015-02-26 11:20 user528 File Added:
2015-02-26 11:39 henry Note Added: 0003918
2015-02-26 17:12 henry Note Added: 0003919
2015-02-27 09:28 user528 Note Added: 0003921
2015-02-27 09:31 henry Status new => resolved
2015-02-27 09:31 henry Resolution open => fixed
2015-02-27 09:31 henry Assigned To => henry
2015-03-13 07:28 user528 Note Added: 0004099
2015-03-13 07:28 user528 Status resolved => feedback
2015-03-13 07:28 user528 Resolution fixed => reopened
2015-03-13 11:55 henry Note Added: 0004107
2015-03-13 11:55 henry Status feedback => resolved
2015-03-13 11:55 henry Resolution reopened => fixed
2015-03-24 00:17 liuhuafei Issue cloned: 0001619
2015-06-18 15:14 henry Note Added: 0004966