furstj/myFoam

How to access density variable in codedFixedValue boundary condition?

Opened this issue · 3 comments

Hello professor furstj,
Recently, I successfully simulated a 2-D high subsonic compressor cascade using your myLusgsFoam solver with uniform total pressure and temperature inlet boundary condition. The "subsonicInletTotal" b.c. is used for inlet velocity boundary.
Nowadays I want to realize a non-uniform inlet total pressure profile along one direction using "codedFixedValue" boundary, where the static pressure is calculated by subtracting the dynamic pressure head (0.5rhomag(U)^2) from the specified total pressure. But I find that there are two forms of density can be access, "rho" and "rho_0".
The question is when the "rho_0" is used, the case can be successfully continued. But when the "rho" is used, the calculation is always divergent. I have no idea for the reason, also I can't find where the "rho_0" is defined. Therefore I hope to get your help. Meanwhile, I want to know whether there is a more reasonable way to define the total inlet pressure profile?
Thanks for professor's help!

Through residuals file "log", I found rho_ 0 may be the initial value in each iteration. But I still have no idea the reason for the divergence using "rho" in the "codedFixedValue" boundary condition, whereas "rho_0" does not.

Hi,

the rho_0 is the initial condition. The reason why rho doesn't work is in the way how the boundary conditions are implemented. Look at myFoam/mySolvers/myDbns/include/updateFields.H. There is (in the following order):

  1. U.correctBoundaryConditions(), then
  2. h.correctBoundaryConditions(); then
  3. p.correctBoundaryConditions(); and finally
  4. rho.boundaryFieldRef() = thermo.psi().boundaryField() * p.boundaryField();

It means that the boundary value of rho is calculated after the boundary value for pressure!

The quantities are tightly coupled in the case of compressible flows and it would be therefore better to calculate all boundary values at once in coupled way. This coupled approach is unfortunately quite difficult to implement for general boundary conditions in OpenFOAM.

BTW, the use of rho_0 instead of rho seems to be simple and working solution at least in the case of steady state problems.

On the other hand, non-uniform total pressure can be achieved simply by using the "standard" totalPressure boundary condition for pressure by setting the value "p0" to your data, i.e. not as "p0 uniform 101.325.0;" but prescribing the desired value for each face. You can use e.g. setExprBoundaryField utility (present e.g. in in OpenFOAM v2112)...

Thanks Professor Furstj.