CFD Online Discussion Forums - Custom boundary condition: unexpected behavior with chtMultiRegionFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Custom boundary condition: unexpected behavior with chtMultiRegionFoam (https://www.cfd-online.com/Forums/openfoam-programming-development/129096-custom-boundary-condition-unexpected-behavior-chtmultiregionfoam.html)

Custom boundary condition: unexpected behavior with chtMultiRegionFoam

Dear OpenFOAM users,
I recently wrote a simple custom boundary condition which computes the fluid velocity at the wall using a formula of the type:

v = f(rho) * normal

While my code yields the expected results using the rhoPimpleFoam solver, everything is messed up when using it with chtMultiRegionFoam, even though I use the exact same parameters for the BC. The near-wall velocity has an inverted orientation, and the upstream flow has a stagnation point, even if the wall velocity is oriented downstream. This is obviously a physical and numerical nonsense ...

Is there something in particular that should be taken care of when using chtMultiRegionFoam with custom BCs?

Parameters with rhoPimpleFoam:

Code:

squareWall

    {

        type            pyrolysisVelocity;

        refValue        0.001;

        massPR          0.001;

        rhoS            -1;

        rhoName         rho;

    }

Result: https://www.dropbox.com/s/n7cfws80w3...Foam_outer.png

With chtMultiRegionFoam:

Code:

fluid_to_solid

            {

                type pyrolysisVelocity;

                refValue 0.001;

                massPR 0.001;

                rhoS -1;

                rhoName rho;

            }

Result: https://www.dropbox.com/s/8t2re4piih...Foam_outer.png

A little progress update: I went through the code and did some additional debugging. The coordinates of velocity on the patch seem to be computed correctly: I checked the values in the output files.

The results with a surfaceNormalFixedValue BC are okay:
https://www.dropbox.com/s/zui2qaz9tz...FixedValue.png
I'm guessing something wrong happens when running the updateCoeffs() method, maybe I shouldn't use the fixedValueFvPatchVectorField::operator==() member function?

I ran an additional test using swak4foam:

Code:

fluid_to_solid

    {

        type groovyBC;

        valueExpression "- 0.001 * normal()";

    }

It gave nearly the same result as obtained with chtMultiRegionFoam:
https://www.dropbox.com/s/0yjr5vca1y..._swak4foam.png

I computed the mass flux using the patchIntegrate utility.

Solver: chtMultiRegionFoam (using BCs for isothermal flow)
Inlet mass flux: 0.000811 kg/s

1. BC: surfaceNormalFixedValue, -0.001
fluid_to_solid mass flux: 0.000463 kg/s
Outlet mass flux: 0.00127 kg/s

2. BC: groovyBC, "- 0.001 * normal()"
fluid_to_solid mass flux: <something>e-20 = 0 kg/s
Outlet mass flux: 0.000811 kg/s

Mass conversation seems to be okay. Something might be going wrong with the phi field, because the value of U at the boundary seems to be fine:

Code:

        value           nonuniform List<vector>

40

(

(0.001 0 0)

(0 -0.001 0)

(0.001 0 0)

(0.001 0 0)

(0.001 0 0)

(0.001 0 0)

(0.001 -0 0)

(0.001 -0 0)

(0.001 -0 0)

(0.001 -0 0)

(0.001 -0 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(0 -0.001 0)

(0 0.001 0)

(-0.001 -0 0)

(0 -0.001 0)

(-0.001 -0 0)

(-0.001 -0 0)

(-0.001 -0 0)

(-0.001 -0 0)

(-0.001 0 0)

(-0.001 0 0)

(-0.001 0 0)

(-0.001 0 0)

(0 0.001 0)

(-0.001 0 0)

)

Any idea of what that could be? I attached my case directory.