[Muerio]

Hello everyone,

I'm currently trying to implement a pair of coded BCs in a simpleFoam case.

The pair of BCs should represent an inlet and outlet with a adapting flow rate on the inner faces parallel to the XZ-Plane (see attachment).

As a first proof of concept I ran the case with flowRateInletVelocity BCs. For the upper face, I defined:

type flowRateInletVelocity;
rhoInlet 1071.2;
value uniform (0 0 0);
massFlowRate 0.0015;

The lower face is the same with a negative massFlowRate. The resulting flowfield can be seen in the attached image.

Now I want to vary the flow rate based on the pressure difference between the two BCs. For this I export the pressure difference using a fieldValueDelta function in the controlDict.

At first I tried to replace the value for the flow rate with a codestream. In the code block I read the last line from the exported fieldValueDelta as the pressure difference and multiplied it with the density resulting in delta_p. From this I calculated my new flow rate flow_rate and tried writing it to the BC using

scalarField mass_flow_rate(1, flow_rate);
mass_flow_rate.writeEntry("", os);

in the code block.
I noticed that the codestream is only evaluated at the first time step, so the flow_rate would not adapt to the change in pressure difference.

Now I tried using codedFixedValue BCs. Again I read the last line from the fieldValueDelta and store it in deltaP. From there I use my scaling factor hydResistance to calculate the volumetric flow rate for the next time step. The scaling vector was determined in the "proof of concept" run. Divided by the patch area, I get the average velocity for the patch which I then write to all faces on the patch.

// Calc volFlow from deltaP 
double hydResistance = 1.1e+08;     //in kg / (m^4 * s)
double volFlow = (deltaP * rho) / hydResistance; //Pressure normalized with rho!    // m^2 / s^2 * kg / m^3 /( kg / (m^4 * s)) = m^3 / s

// Avg. velocity on Patch equals volFlow / patchArea
double velocity = volFlow / area_1;                                                 // m^3 / (s * m^2) = m / s

forAll(Cf, faceI) // Write avg. velocity to Patch faces
{
    U [faceI] = vector(0, velocity, 0);
}

Running the case results in very large changes in the pressure difference and ultimately fails when the values get too large.

# Source1         : pressureDifferencePatch1.region1
# Source2         : pressureDifferencePatch1.region2
# Operation       : subtract
# Time            	p
1                 	1.2744304410e+01
2                 	1.0534881848e+00
3                 	-2.1438541297e+03
4                 	5.3147321704e+04
5                 	3.6129157963e+06
6                 	-2.1841268573e+10
7                 	-3.1164728773e+13

For all cases, the settings for the BCs for the pressure remain the same:
All BCs use zeroGradient except for the main outlet at the top of the geometry, which is a fixedValue of uniform 0.

So now I don't know where to go from here and have a few questions:

1. Is there a way to use codestream in a flowRateInletVelocity condition so it is evaluated in every iteration?
2. Is there a mistake in my codedFixedValue condition I haven't noticed?
3. What could be changed in the codedFixedValue condition to stop the pressure values from going crazy?
4. Is there another possibility to achieve my goal?

Thank you for your help!