CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

New BC with differential equation

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

Like Tree6Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 25, 2020, 05:56
Default New BC with differential equation
  #1
Neb
New Member
 
Join Date: Mar 2020
Posts: 28
Rep Power: 2
Neb is on a distinguished road
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 is offline   Reply With Quote

Old   June 29, 2020, 09:28
Default
  #2
Neb
New Member
 
Join Date: Mar 2020
Posts: 28
Rep Power: 2
Neb is on a distinguished road
can no one help me? I just need some ideas for the derivative
Neb is offline   Reply With Quote

Old   June 30, 2020, 08:08
Default
  #3
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by Neb View Post
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 View Post
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
RGS is offline   Reply With Quote

Old   June 30, 2020, 09:16
Default
  #4
Neb
New Member
 
Join Date: Mar 2020
Posts: 28
Rep Power: 2
Neb is on a distinguished road
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
Neb is offline   Reply With Quote

Old   June 30, 2020, 10:59
Default
  #5
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by Neb View Post
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.


Hopefully, this post will help you with the derivative:

Output of velocity gradients
RGS is offline   Reply With Quote

Old   June 30, 2020, 11:11
Default
  #6
Neb
New Member
 
Join Date: Mar 2020
Posts: 28
Rep Power: 2
Neb is on a distinguished road
Thank you for the link. What do you have to do? Maybe I can help you
Neb is offline   Reply With Quote

Old   July 1, 2020, 13:38
Default
  #7
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
crubio.abujas is on a distinguished road
Quote:
Originally Posted by Neb View Post
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!
Gerry Kan likes this.
crubio.abujas is offline   Reply With Quote

Old   July 1, 2020, 16:14
Default
  #8
Neb
New Member
 
Join Date: Mar 2020
Posts: 28
Rep Power: 2
Neb is on a distinguished road
Yes! It was just what I needed! Thank you
Neb is offline   Reply With Quote

Old   July 2, 2020, 02:31
Default
  #9
New Member
 
Join Date: Jun 2020
Posts: 1
Rep Power: 0
Mihan is on a distinguished road
Thank you very much for your helps and guides.
Mihan is offline   Reply With Quote

Old   July 2, 2020, 09:04
Default
  #10
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by Neb View Post
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
RGS is offline   Reply With Quote

Old   July 2, 2020, 13:34
Default
  #11
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
crubio.abujas is on a distinguished road
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.
Gerry Kan and RGS like this.
crubio.abujas is offline   Reply With Quote

Old   July 3, 2020, 09:58
Default
  #12
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by crubio.abujas View Post
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.
RGS is offline   Reply With Quote

Old   July 3, 2020, 11:39
Default
  #13
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
crubio.abujas is on a distinguished road
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 View Post
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 View Post
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 View Post
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 likes this.
crubio.abujas is offline   Reply With Quote

Old   July 6, 2020, 08:08
Default
  #14
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by crubio.abujas View Post
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
RGS is offline   Reply With Quote

Old   July 6, 2020, 11:29
Default
  #15
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
crubio.abujas is on a distinguished road
Glad to know it helped to you! Concerning the other questions you have:

Quote:
Originally Posted by RGS View Post
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 View Post
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!
crubio.abujas is offline   Reply With Quote

Old   July 7, 2020, 04:46
Default
  #16
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by crubio.abujas View Post
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?
RGS is offline   Reply With Quote

Old   July 7, 2020, 06:17
Default
  #17
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
crubio.abujas is on a distinguished road
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!




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.
crubio.abujas is offline   Reply With Quote

Old   July 7, 2020, 09:30
Default
  #18
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
Quote:
Originally Posted by crubio.abujas View Post
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!




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.
RGS is offline   Reply With Quote

Old   July 7, 2020, 10:36
Default
  #19
Member
 
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 84
Rep Power: 5
crubio.abujas is on a distinguished road
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!
crubio.abujas is offline   Reply With Quote

Old   July 13, 2020, 08:48
Default
  #20
RGS
New Member
 
Rohit George Sebastian
Join Date: May 2017
Posts: 17
Rep Power: 5
RGS is on a distinguished road
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?
RGS is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Explicit Runge-Kutta of a Differential Equation HakikiCanakkaleli OpenFOAM Programming & Development 1 September 25, 2014 06:18
Basic question about writing a differential equation related to finite element RbBb Main CFD Forum 1 April 10, 2014 03:50
Non-linearity Pressure Equation -- PISO algorithm gdeneyer OpenFOAM Programming & Development 1 August 23, 2012 05:19
DEFINE_GRID_MOTION using differential equation Ramsey FLUENT 1 January 14, 2011 03:13
differential equation for turbulent viscosity mateus FLUENT 0 March 4, 2005 01:45


All times are GMT -4. The time now is 10:35.