CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   New BC with differential equation (https://www.cfd-online.com/Forums/openfoam-programming-development/228271-new-bc-differential-equation.html)

Neb June 25, 2020 06:56

New BC with differential equation
 
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?

Neb June 29, 2020 10:28

can no one help me? I just need some ideas for the derivative

RGS June 30, 2020 09:08

Quote:

Originally Posted by Neb (Post 775967)
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 (Post 775967)
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

Neb June 30, 2020 10:16

Hi, the code I have so far is in this post. But now I need to be able to implement the derivative.

https://www.cfd-online.com/Forums/op...tml#post773341

RGS June 30, 2020 11:59

Quote:

Originally Posted by Neb (Post 776441)
Hi, the code I have so far is in this post. But now I need to be able to implement the derivative.

https://www.cfd-online.com/Forums/op...tml#post773341


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


Hopefully, this post will help you with the derivative:

https://www.cfd-online.com/Forums/op...tml#post328925

Neb June 30, 2020 12:11

Thank you for the link. What do you have to do? Maybe I can help you

crubio.abujas July 1, 2020 14:38

Quote:

Originally Posted by Neb (Post 775967)
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:

\frac{d}{dt}\left(P-R_1\phi\right) + \frac{1}{R_2C}\left(P_0-R_1\phi\right)=\phi{}C
You define \frac{d}{dt}\left(P-R_1\phi\right) = \frac{1}{\Delta{}t}\left[P^t-P^(t-1)-R_1(\phi^t-\phi^{t-1})\right].

And then just find P^t = f(P^{t-1}, \phi^t, \phi^{t-1}).


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!

Neb July 1, 2020 17:14

Yes! It was just what I needed! Thank you

Mihan July 2, 2020 03:31

Thank you very much for your helps and guides.

RGS July 2, 2020 10:04

Quote:

Originally Posted by Neb (Post 776454)
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:


\left( \textbf{e} - \nabla \textbf{u} \right) . \textbf{n} = 0


where,

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

crubio.abujas July 2, 2020 14:34

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<volScalarField>("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<vector>::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.

RGS July 3, 2020 10:58

Quote:

Originally Posted by crubio.abujas (Post 776710)
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<volScalarField>("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<vector>::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.

Could you please help me understand your code a bit better? I didn't fully understand the lines


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.

crubio.abujas July 3, 2020 12:39

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 (Post 776801)
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 (Post 776801)
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 (Post 776801)
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!

RGS July 6, 2020 09:08

Quote:

Originally Posted by crubio.abujas (Post 776811)
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<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];
 }

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.

Link: https://cpp.openfoam.org/v7/classFoa...81a2dd93376d82

crubio.abujas July 6, 2020 12:29

Glad to know it helped to you! Concerning the other questions you have:

Quote:

Originally Posted by RGS (Post 776978)
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:
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 (Post 776978)
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!

RGS July 7, 2020 05:46

Quote:

Originally Posted by crubio.abujas (Post 777004)
Well in fact when you created the tortuosity field you've used a volVectorField type. This type is defined as:
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.


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?

crubio.abujas July 7, 2020 07:17

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 \bf{\nabla{}u} is that, when subtracted to \bf{e} the resultant vector is perpendicular to \bf{n}. But there are infinite number of vectors that can fulfill that request!


https://i.ibb.co/QFh601x/condition-visualized.png

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 July 7, 2020 10:30

Quote:

Originally Posted by crubio.abujas (Post 777069)
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 \bf{\nabla{}u} is that, when subtracted to \bf{e} the resultant vector is perpendicular to \bf{n}. But there are infinite number of vectors that can fulfill that request!


https://i.ibb.co/QFh601x/condition-visualized.png

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.

crubio.abujas July 7, 2020 11:36

Thinking about it, the condition \nabla{}u\cdot{}n = e\cdot{}n 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!

RGS July 13, 2020 09:48

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?


All times are GMT -4. The time now is 07:57.