[sahm]

Greetigns fellow openfoamers

I have developed a solver for my project, and wanted to implement a boundary condition on it. Here's the problem:

I need to solve this equation:

mdot = ∫ K*P/(R*T*Mu) *Grad(P) dA

mdot is the boundary condition that user defines (constant mass flow rate) and I need to solve this equation for Grad(P). The problem is that K, T, P and Grad(P) are not uniform on the boundary, but boundary pressure is.
How can I implement an integral like that?

One method I can think of, is to discretize Grad(P) ((P_bc-P)/Delta) and assume P_bc is constant everywhere, then solve that integral for P_bc, and set is as boundary condition. In that case, I need to calculate the following equation:

-mdot + ∫ K*P^2/(R*T*Δ*Mu ) dA
P_bc = —————————————————————
∫ K*P/(R*T*Δ*Mu ) dA

How can I calculate such integrals on the boundary cells? I tried to find the integral function, but couldn't find a good answer.

[chriss85]

The integral basically comes down to a sum over the boundary faces. Boundary faces are accessed with mesh.boundary()[patchI], or for a field, field.boundaryField()[patchI]. You can get information, for example the face area by mesh.boundary()[patchI].magSf(). Then it comes down to a loop over all involved faces in which you add up your integrand.

[sahm]

Thanks Chriss85
Here are some questions:

How can I get patchID for the current patch? I don't want to use the findPatchID(" Patch Name "); since it will be working for 1 name only. I am defining this as a boundary condition and should work with any patch name

I tried using
const polyPatch& cPatch = this->patch();
and got this error:

MassFlowFvPatchScalarField.C:175:43: error: invalid initialization of reference of type ‘const Foam:: PolyPatch&’ from expression of type ‘const Foam::fvPatch’
const polyPatch& cPatch = this->patch();

I'm confused with types defined in OF.

[chriss85]

You should be able to search the patch index in the list of patches, i.e. mesh.boundary(). There is also some indexing function you can use, but I don't recall it from my head.
For your second problem, try this->patch().patch(). polyPatch is a different class than fvPatch. Look in the respective header files to see the members of the classes.

[jherb]

This is an example of a boundary condition creating a Field based on the values on the patch:
https://github.com/OpenFOAM/OpenFOAM...chField.C#L111

[sahm]

Joachim Herb and Chriss

Thanks a lot for your comments. I made my boundary condition based on the fixedMean boundary condition as Joachim suggested. But I have more questions and problems to solve. I would appreciate if anyone could guide me through these questions. I will also include information that others might have questions about:

I got access to my K, T & P Using lookupObject:
eg.

    const volScalarField& K = //permeability from the solver
    this->db().objectRegistry::lookupObject<volScalarField> ("K_e");

I got the vectors from face center to cell centeres on the boundary Copied from fixedDisplacementFvPatchVectorField.C ) and dA

 vectorField delta = this->patch().delta();
 scalarField deltaC = this->patch().deltaCoeffs();
 scalarField dA = this->patch().magSf();

and mag (delta) should be equal to 1/deltaC (which is used to calculated gradients).

I can get the values on the patch:

    scalarField PatchValues(this->patchInternalField());

which is equal to P_.boundaryField()[patchindex].patchInternalField() . I got the index value using scalar

patchindex = this->patch().index();

But this part is unnecessary.

I'm calculating the gradient on the patch as:

 scalarField snGradP = (PatchValues-P_) * deltaC;

which apparently is how openFOAM does it. I tried to use functions like fvc::grad and snGrad(), but they require extra work. Any help on using these would be appreciated.

I'm calculating mdot using this code:

    scalar m_dot_calc = gSum (dA*snGradP*K*P_/( mu*Rbar*T));

But the problem is that m_dot_calc is being reported as negative. The reason is some terms in that summation come out as negative. Upon further check by printing values to log file, I noticed that snGradP*P_ gives negative values. Since snGradP is always positive, it means the P_ values come out as negative. I like to know why? How can I show the values from P_ that are being used in that equation? As I mentioned above P_.boundaryField()[patchindex].patchInternalField() are the same as PatchValues. So what values are being used for snGradP and m_dot ?

This is also how I'm calculating the equation in my first post (I mean using gSum and the rest). Another way that I can implement my boundary condition is to calculate m_dot and compare it with the user defined value, and then adjust the pressure on the boundary regarding the comparison. But even m_dot isn't working.

[sahm]

I asked openFOAM to write P file after, every time step, and noticed that P becomes negative in some cells. I changed the boundary condition to fixedValue, and noticed that the same thing is happening. So I have to make sure that my Solver is calculating P correctly. Interestingly the solver gave quite good results for the cases that I have run before , but for initial time steps, the P is coming out as negative which might be due to the source term I have in the mass transfer equation. (That's why the mass flow rate was being reported as negative) I will try to figure things out with the solver, but I am still a little confused on what happens in the boundary condition. The programmers guide is Ultra Vague on this matter. I think we should write some more guides and tutorials on how to do things when it comes to programming boundary conditions or other parts. The tutorials that I have seen are quite basic comparing to what I have seen people are trying to do.

[sahm]

Hello fellow FOAMers.

Thanks for your helps so far, I have made some progress in my work and my research. I have some more questions that I will describe bellow.

I finally made my boundary condition as described bellow, and I could run a case to test it. It seems it's working fine to some extent, however I think it is calculating the mass flow rate wrong.

Here's how the boundary condition works:

I calculate the m_dot as:
m_dot = ∫ K*P/(R*T*Mu) *Grad(P) dA

which is translated as this code:

    surfaceScalarField PGrad_ = fvc::snGrad(P_);// snGrad() is face Normal Gradient 

    scalarField m_dot_terms = dA*PGrad_.boundaryField()[patchindex]*K*P_/( mu*Rbar*T);
    scalar m_dot_e = gSum (m_dot_terms);

Then the code compares this m_dot with the prescribed value (input of the boundary condition) and adjusts the pressure so that the mass flow rate stays constant.

The code reports this mass flow rate and the reported value is the same as the prescribed value which means the code is working fine. However, when I calculated the amount of mass in the reactor and got a derivative of that, it seems like the mass flow rate into the system is about 5 times lower than it should be. I suspect that the m_dot is being calculated wrongly. Can someone tell me if I am calculating the pressure gradient on the boundary and m_dot correctly?

Thanks a lot.