# New BC with differential equation

 Register Blogs Members List Search Today's Posts Mark Forums Read

 June 25, 2020, 05:56 New BC with differential equation #1 New Member   Join Date: Mar 2020 Posts: 28 Rep Power: 2 Hi everyone, I'm new to OpenFoam, I need your help. I need to make a new boundary condition and I am using codedFixedValue. Unfortunately I can't solve it. it is a differential equation like this: d(P − R1*Q)/dt +(1/R2*C) *(P(t=0) − R1*Q) =Q*C P and Q are respectively the pressure and the flow rate in the same outlet where I want to apply the equation. R1, C, R2 are constants. Would anyone be kind to help me? So far I have managed to implement the constants, the average flow rate in the outlet and the time read, but I have stuck on the derivative. Another question: the pressure I enter must already be divided by density or does OpenFoam do it?

 June 29, 2020, 09:28 #2 New Member   Join Date: Mar 2020 Posts: 28 Rep Power: 2 can no one help me? I just need some ideas for the derivative

June 30, 2020, 08:08
#3
New Member

Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
Quote:
 Originally Posted by Neb Hi everyone, I'm new to OpenFoam, I need your help. I need to make a new boundary condition and I am using codedFixedValue. Unfortunately I can't solve it. it is a differential equation like this: d(P − R1*Q)/dt +(1/R2*C) *(P(t=0) − R1*Q) =Q*C P and Q are respectively the pressure and the flow rate in the same outlet where I want to apply the equation. R1, C, R2 are constants. Would anyone be kind to help me? So far I have managed to implement the constants, the average flow rate in the outlet and the time read, but I have stuck on the derivative.

Hi, I am also trying to implement a custom boundary condition (for the first time), so we may be able to help each other. Can you show me the code you have so far?

Quote:
 Originally Posted by Neb Another question: the pressure I enter must already be divided by density or does OpenFoam do it?

That depends on the solver.

p_rgh = p - rho*gh

 June 30, 2020, 09:16 #4 New Member   Join Date: Mar 2020 Posts: 28 Rep Power: 2 Hi, the code I have so far is in this post. But now I need to be able to implement the derivative. verification of a new BC whit codedfixedValue

June 30, 2020, 10:59
#5
New Member

Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
Quote:
 Originally Posted by Neb Hi, the code I have so far is in this post. But now I need to be able to implement the derivative. verification of a new BC whit codedfixedValue

I think you are a more advanced user than I am.

 June 30, 2020, 11:11 #6 New Member   Join Date: Mar 2020 Posts: 28 Rep Power: 2 Thank you for the link. What do you have to do? Maybe I can help you

July 1, 2020, 13:38
#7
Member

Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
Quote:
 Originally Posted by Neb Hi everyone, I'm new to OpenFoam, I need your help. I need to make a new boundary condition and I am using codedFixedValue. Unfortunately I can't solve it. it is a differential equation like this: d(P − R1*Q)/dt +(1/R2*C) *(P(t=0) − R1*Q) =Q*C P and Q are respectively the pressure and the flow rate in the same outlet where I want to apply the equation. R1, C, R2 are constants.
Hello Neb!

I am aware that different temporal integration schemes are available inside OF libraries, but as far as I've checked they are designed to be applied onto complete fields instead of patch fields. Nevertheless as a straight solution you can try to implement a simple finite different approach using the values of previous time steps. On your equation:

You define .

And then just find .

Here is a code that does just that. That shall at least compile!

Code:
code
#{

const scalar dt = this ->db().time().deltaTValue(); //delta t reading for integration

// Mass flow field
const surfaceScalarField& phi =
db().lookupObject<surfaceScalarField>("phi");
const fvsPatchField<scalar>& phip =
patch().patchField<surfaceScalarField,scalar> (phi);
const fvsPatchField<scalar>& phip_old =
patch().patchField<surfaceScalarField,scalar> (phi.oldTime());

// Current patch p field
scalarField& p = *this;

// Constants
const scalar C   (4.82e0);
const scalar R1  (2.79245);
const scalar R2  (1.7);
const scalar p0 (1.07);

// Set the new values (overwrites the previoues p value)
forAll(p, facei) {
p[facei] =
dt*(phip_old[facei]*C - 1/R2/C*(p0-R1*phip_old[facei])) +
p[facei] - R1*(phip_old[facei] - phip[facei]);
}

#};
Let's hope you can advance with this new approach!

 July 1, 2020, 16:14 #8 New Member   Join Date: Mar 2020 Posts: 28 Rep Power: 2 Yes! It was just what I needed! Thank you

 July 2, 2020, 02:31 #9 New Member   Join Date: Jun 2020 Posts: 1 Rep Power: 0 Thank you very much for your helps and guides.

July 2, 2020, 09:04
#10
New Member

Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
Quote:
 Originally Posted by Neb Thank you for the link. What do you have to do? Maybe I can help you

Thank you for offering to help.

I am trying to implement the following BC:

where,

e is a unit vector
u is the tortuosity
n is the surface normal vector

 July 2, 2020, 13:34 #11 Member   Carlos Rubio Abujas Join Date: Jan 2018 Location: Spain Posts: 84 Rep Power: 5 Hello Rohit, May I assume that the tortuosity field (u) is already defined in the entire domain as part of your problem? I'm guessing also that the e unit vector is constant and set on the boundary condition, Is that right? In such case you may want to adapt the following code to your needs: Code: inlet { type codedFixedValue; value uniform 0; name vectorProyection; code #{ // Get the reference of the mesh const fvMesh& mesh = patch().boundaryMesh().mesh(); // Get the reference of the tortuosity const volScalarField& tort = db().lookupObject("u"); // Create a dummy object to store the gradient. // You shall adapt the dimensionSet to match your units volVectorField gTort ( IOobject ( "grad(u)", db().time().timePath(), db(), IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedVector ( "gTort0", dimensionSet(0, -1, 0, 0, 0, 0, 0), pTraits::zero ) ); // Set the unit vector const vector e(0.5, 0.5, 0); // Calculate the gradient of the tortuosity gTort = fvc::grad(tort); // Get a reference of the gradient field in current patch const vectorField& gTort_p = gTort.boundaryField()[patch().index()]; // Set the equation on each face of the patch scalarField& field = *this; forAll(field, facei) { // The operator & if the inner product field[facei] = (e - gTort_p[facei]) & patch().nf()()[facei]; } #}; codeInclude #{ #include "fvCFD.H" #}; codeOptions #{ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude #} }; It calculate the gradient on the entire domain just for using it in the boundary, which is a bit of a resources waste, but it shall work. Gerry Kan and RGS like this. July 3, 2020, 09:58 #12 New Member Rohit George Sebastian Join Date: May 2017 Posts: 17 Rep Power: 5 Quote:  Originally Posted by crubio.abujas Hello Rohit, May I assume that the tortuosity field (u) is already defined in the entire domain as part of your problem? I'm guessing also that the e unit vector is constant and set on the boundary condition, Is that right? In such case you may want to adapt the following code to your needs: Code: inlet { type codedFixedValue; value uniform 0; name vectorProyection; code #{ // Get the reference of the mesh const fvMesh& mesh = patch().boundaryMesh().mesh(); // Get the reference of the tortuosity const volScalarField& tort = db().lookupObject("u"); // Create a dummy object to store the gradient. // You shall adapt the dimensionSet to match your units volVectorField gTort ( IOobject ( "grad(u)", db().time().timePath(), db(), IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedVector ( "gTort0", dimensionSet(0, -1, 0, 0, 0, 0, 0), pTraits::zero ) ); // Set the unit vector const vector e(0.5, 0.5, 0); // Calculate the gradient of the tortuosity gTort = fvc::grad(tort); // Get a reference of the gradient field in current patch const vectorField& gTort_p = gTort.boundaryField()[patch().index()]; // Set the equation on each face of the patch scalarField& field = *this; forAll(field, facei) { // The operator & if the inner product field[facei] = (e - gTort_p[facei]) & patch().nf()()[facei]; } #}; codeInclude #{ #include "fvCFD.H" #}; codeOptions #{ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude #} }; It calculate the gradient on the entire domain just for using it in the boundary, which is a bit of a resources waste, but it shall work.

Hello Carlos,

Thanks a lot for this code, it helps me a lot. I still have to do a few things before I can run a simulation, but I will let you know if everything works.

You were right in both your assumptions - the tortuosity field is already defined in the entire domain, and the unit vector is a constant that only appears in the boundary condition.

Code:
    // Set the equation on each face of the patch
scalarField& field = *this;
forAll(field, facei)
{
// The operator & if the inner product
field[facei] = (e - gTort_p[facei]) & patch().nf()()[facei];
}
is facei a standard template variable in openfoam?
what does the ()[facei] do?
do we not need to specify somewhere that the expression should evaluate to zero?

Also, how do you know which libraries to include in codeOptions?

I hope these questions are not stupid or too basic.

July 3, 2020, 11:39
#13
Member

Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
No stupid questions at all! It can be quite frustrating to understand these kind of things (I've been there), even now I'm just trying to make myself comfortable with how the code works, so maybe my explanation is not quite exact. If any of you can improve this explanation do not hesitate in letting me know!

Code:
  scalarField& field = *this;
First you have to understand that is polymorphism. The code here executed lives inside an object of type fvsPatchField , but by calling it this way you tell the code to treat it as a Field object. It can behave as such because the class has inhered from Field.

Quote:
 Originally Posted by RGS is facei a standard template variable in openfoam?
No, OF has no builtin variable called facei. You may be wondering then, how in earth does the code know what is facei if it has never been defined?. The answer... it has been defined in the forAll statement.

Let me explain this a little more. forAll is not a function nor a method. It's a macro. Typically macros are writer in upper letters to avoid confusion, but smart developers of OF consider it is more readable this way. The macro is defined as:
Code:
 #define forAll(list, i) \
for (Foam::label i=0; i<(list).size(); i++)
Before compilation it replaces the variable value of "i" to whatever text you provided, in this case facei. Macros are used to ease the writing process and improve the readability, but if you are not aware of them it can obscure the code. So what is really happening under the hood is:
Code:
for (Foam::label facei=0; facei<(field).size(); facei++)
{
field[facei] = (e - gTort_p[facei]) & patch().nf()()[facei];
}
Now the definition of facei shall be evident.

Quote:
 Originally Posted by RGS what does the ()[facei] do?
What you want is the normal vector of that face. Recall that your starting position (this) is a fvsPatchField, so you have to use the methods available by this class. First you call the patch() method, which returns you a fvPatch object. This object connect all the geometrical information of the patch, such as nf. But if you watch carefully the method nf() (here) it does not return a vector field, but a tmp<vectorField>. tmp objects are related with memory management and other issues, but the idea is that they are pointing to the relevant information, but they are NOT the relevant information by themselves. You need to do an extra step to get the normals, and that is to call the operator(). Then you have a list of all the faces and you select them individually by the operator[].

I've compressed all these step by chaining them together, but what is really going on can be written as follows:
Code:
// Get the temporal object nf
const tmp<vectorField> tmpNf = patch().nf();

// Get the vector field associated with the temporal object
const vectorField& nf_p = tmpNf();

for (Foam::label facei=0; facei<(field).size(); facei++)
{
field[facei] = (e - gTort_p[facei]) & nf_p[facei];
}
Quote:
 Originally Posted by RGS Also, how do you know which libraries to include in codeOptions?
To be honest, there are for sure quite more libraries attached than strictly required. I've just copied the ones appearing into simpleFoam. You can check it into the source code.

I hope that helps you understand what is going on!

July 6, 2020, 08:08
#14
New Member

Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
Quote:
 Originally Posted by crubio.abujas No stupid questions at all! It can be quite frustrating to understand these kind of things (I've been there), even now I'm just trying to make myself comfortable with how the code works, so maybe my explanation is not quite exact. If any of you can improve this explanation do not hesitate in letting me know! Code:  scalarField& field = *this; First you have to understand that is polymorphism. The code here executed lives inside an object of type fvsPatchField , but by calling it this way you tell the code to treat it as a Field object. It can behave as such because the class has inhered from Field. No, OF has no builtin variable called facei. You may be wondering then, how in earth does the code know what is facei if it has never been defined?. The answer... it has been defined in the forAll statement. Let me explain this a little more. forAll is not a function nor a method. It's a macro. Typically macros are writer in upper letters to avoid confusion, but smart developers of OF consider it is more readable this way. The macro is defined as: Code:  #define forAll(list, i) \ for (Foam::label i=0; i<(list).size(); i++) Before compilation it replaces the variable value of "i" to whatever text you provided, in this case facei. Macros are used to ease the writing process and improve the readability, but if you are not aware of them it can obscure the code. So what is really happening under the hood is: Code: for (Foam::label facei=0; facei<(field).size(); facei++) { field[facei] = (e - gTort_p[facei]) & patch().nf()()[facei]; } Now the definition of facei shall be evident. What you want is the normal vector of that face. Recall that your starting position (this) is a fvsPatchField, so you have to use the methods available by this class. First you call the patch() method, which returns you a fvPatch object. This object connect all the geometrical information of the patch, such as nf. But if you watch carefully the method nf() (here) it does not return a vector field, but a tmp. tmp objects are related with memory management and other issues, but the idea is that they are pointing to the relevant information, but they are NOT the relevant information by themselves. You need to do an extra step to get the normals, and that is to call the operator(). Then you have a list of all the faces and you select them individually by the operator[]. I've compressed all these step by chaining them together, but what is really going on can be written as follows: Code: // Get the temporal object nf const tmp tmpNf = patch().nf(); // Get the vector field associated with the temporal object const vectorField& nf_p = tmpNf(); for (Foam::label facei=0; facei<(field).size(); facei++) { field[facei] = (e - gTort_p[facei]) & nf_p[facei]; } To be honest, there are for sure quite more libraries attached than strictly required. I've just copied the ones appearing into simpleFoam. You can check it into the source code. I hope that helps you understand what is going on!

Thank you for your reply. It was just detailed enough for me to understand and to help me to navigate through the source code.

I have one last question. Should I not somehow specify that the expression has to evaluate to zero at the surfaces? Or is that what the code already means?

I also have a comment that might help other beginners like me. The operator () has been redefined in OpenFOAM. Not knowing this was confusing for me in the beginning.

July 6, 2020, 11:29
#15
Member

Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
Glad to know it helped to you! Concerning the other questions you have:

Quote:
 Originally Posted by RGS I have one last question. Should I not somehow specify that the expression has to evaluate to zero at the surfaces? Or is that what the code already means?
Well in fact when you created the tortuosity field you've used a volVectorField type. This type is defined as:
Code:
typedef GeometricField<vector, fvPatchField, volMesh> volVectorField
So the class is a GeometricField object, particularized with a set components in a template. The first one is the type of data contained in the field (vector in this case), the second (fvPatchField) is related with the boundary conditions and how to treat them, and the last one is the way it stores the information. volMesh stores the information in the cell centers, in contrast surfaceMesh stores them in the face centers and pointMesh does it in the vertices.

The fvPatchField is read from the boundary condition definition. When you define a boundary condition the file has this form:
Code:
dimensions      [0 2 -2 0 0 0 0];

// This is used to fill the internalField
internalField   uniform 0;

// This is used to fillup the fvPatchField object
boundaryField
{
<patchName>
{ ...  }
}
So there you go, when you get the volVectorField from the objectRegistry(which have used the BC definition file) it already knows all the information on how to treat the boundaries for the gradient calculation. If you want to set the boundaries to 0, you just have to define it the file 0/u.

Quote:
 Originally Posted by RGS I also have a comment that might help other beginners like me. The operator () has been redefined in OpenFOAM. Not knowing this was confusing for me in the beginning.
Yes, in C++ you can redefine how the objects are treated with special function called operators. Each operation on the object (+, -, =, ==, >, ...) can be redefined inside to suit your needs. To understand them you may have to look around on the definition on the class and the one it has inherited.

It is a good point to remark these doubts as for sure more people are struggling with the same issues. I hope my answers serve you, let me know if is it still unclear to you!

July 7, 2020, 04:46
#16
New Member

Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
Quote:
 Originally Posted by crubio.abujas Well in fact when you created the tortuosity field you've used a volVectorField type. This type is defined as: Code: typedef GeometricField volVectorField So the class is a GeometricField object, particularized with a set components in a template. The first one is the type of data contained in the field (vector in this case), the second (fvPatchField) is related with the boundary conditions and how to treat them, and the last one is the way it stores the information. volMesh stores the information in the cell centers, in contrast surfaceMesh stores them in the face centers and pointMesh does it in the vertices. The fvPatchField is read from the boundary condition definition. When you define a boundary condition the file has this form: Code: dimensions [0 2 -2 0 0 0 0]; // This is used to fill the internalField internalField uniform 0; // This is used to fillup the fvPatchField object boundaryField { { ... } } So there you go, when you get the volVectorField from the objectRegistry(which have used the BC definition file) it already knows all the information on how to treat the boundaries for the gradient calculation. If you want to set the boundaries to 0, you just have to define it the file 0/u.

Hello! Thank you for the explanation about GeometricField objects - it helps me understand better what OpenFOAM is doing.

I think you misunderstood what I am trying to implement. I don't want to set the value of the field at the boundaries to zero. I want the expression to evaluate to zero at the boundaries. So the tortuosity or u should not be zero at the boundaries. I would like OpenFOAM to assign a value to u such that the expression becomes zero at the boundaries. Is this even possible?

 July 7, 2020, 06:17 #17 Member   Carlos Rubio Abujas Join Date: Jan 2018 Location: Spain Posts: 84 Rep Power: 5 Ok, now I'm getting you. Sorry, I totally missed the point. In such a case what you want is to change u in the boundaries, so the field validate your condition. Right? I've been thinking about it and I think you need further information. Let me explain: In this case you need to work with a boundary condition in u (tortuosity field) which is the one related to the condition. You can set a gradient into this field, but your condition does not totally define the expected value, only part of its information. What you're asking to is that, when subtracted to the resultant vector is perpendicular to . But there are infinite number of vectors that can fulfill that request! You can add a corrector term to other expression to ensure that this condition is satisfied, but you need further definition on this gradient. Also you need to take into account that the relevant information to OpenFoam is the value of the gradient projected on the line connecting the face center and the cell center. This is the factor that will affect the rest of your field. If you figure out a way to define it, then it can be implemented. RGS likes this.

July 7, 2020, 09:30
#18
New Member

Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
Quote:
 Originally Posted by crubio.abujas Ok, now I'm getting you. Sorry, I totally missed the point. In such a case what you want is to change u in the boundaries, so the field validate your condition. Right? I've been thinking about it and I think you need further information. Let me explain: In this case you need to work with a boundary condition in u (tortuosity field) which is the one related to the condition. You can set a gradient into this field, but your condition does not totally define the expected value, only part of its information. What you're asking to is that, when subtracted to the resultant vector is perpendicular to . But there are infinite number of vectors that can fulfill that request! You can add a corrector term to other expression to ensure that this condition is satisfied, but you need further definition on this gradient. Also you need to take into account that the relevant information to OpenFoam is the value of the gradient projected on the line connecting the face center and the cell center. This is the factor that will affect the rest of your field. If you figure out a way to define it, then it can be implemented.

Yes, that's true! I didn't think about that. Thank you! We are still trying to understand the problem (the equations to be solved and the BCs), so it is very much possible that something is missing or even wrong.

Do you by any chance have an example of a BC in which the value inscribed on the field was somehow determined by OpenFOAM? I am very curious to see what the implementation would look like.

 July 7, 2020, 10:36 #19 Member   Carlos Rubio Abujas Join Date: Jan 2018 Location: Spain Posts: 84 Rep Power: 5 Thinking about it, the condition can be applied if the normal is pointing to the cell center. If you're working with simple meshes there's a good chance that this is the case. If you have a big incidence of non-orthogonality then this approach could lead to some numerical errors, but maybe these error can be assumed. You can try this approach as a first approximation and see is it fits your results. I've been recently working through the book "The Finite Volume Method in Computational Fluid Dynamics", there is a chapter which tell you how OpenFoam deals with the non-orthogonality and some strategies to reduce the errorr. The book is openly available, so feel free to take a look on it!

 July 13, 2020, 08:48 #20 New Member   Rohit George Sebastian Join Date: May 2017 Posts: 17 Rep Power: 5 Hello, I was finally able to run a simulation with the code that you sent me, and I got the following error: Code: --> FOAM FATAL ERROR: Attempt to cast type patch to type lduInterface From function To& Foam::refCast(From&) [with To = const Foam::lduInterface; From = const Foam::fvPatch] in file /home/ubuntu/OpenFOAM/OpenFOAM-7/src/OpenFOAM/lnInclude/typeInfo.H at line 114. FOAM aborting Do you know what might be the cause of the error?