|
[Sponsors] |
June 25, 2020, 05:56 |
New BC with differential equation
|
#1 |
Member
Join Date: Mar 2020
Posts: 66
Rep Power: 6 |
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 |
Member
Join Date: Mar 2020
Posts: 66
Rep Power: 6 |
can no one help me? I just need some ideas for the derivative
|
|
June 30, 2020, 08:08 |
|
#3 | ||
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
Quote:
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:
That depends on the solver. p_rgh = p - rho*gh |
|||
June 30, 2020, 09:16 |
|
#4 |
Member
Join Date: Mar 2020
Posts: 66
Rep Power: 6 |
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 | |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
Quote:
I think you are a more advanced user than I am. Hopefully, this post will help you with the derivative: Output of velocity gradients |
||
June 30, 2020, 11:11 |
|
#6 |
Member
Join Date: Mar 2020
Posts: 66
Rep Power: 6 |
Thank you for the link. What do you have to do? Maybe I can help you
|
|
July 1, 2020, 13:38 |
|
#7 | |
Senior Member
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9 |
Quote:
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]); } #}; |
||
July 1, 2020, 16:14 |
|
#8 |
Member
Join Date: Mar 2020
Posts: 66
Rep Power: 6 |
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 |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
||
July 2, 2020, 13:34 |
|
#11 |
Senior Member
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9 |
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 #} }; |
|
July 3, 2020, 09:58 |
|
#12 | |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
Quote:
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]; } 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 |
Senior Member
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9 |
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; 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++) Code:
for (Foam::label facei=0; facei<(field).size(); facei++) { field[facei] = (e - gTort_p[facei]) & patch().nf()()[facei]; } 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]; } I hope that helps you understand what is going on! |
|
July 6, 2020, 08:08 |
|
#14 | |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
Quote:
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 |
||
July 6, 2020, 11:29 |
|
#15 | ||
Senior Member
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9 |
Glad to know it helped to you! Concerning the other questions you have:
Quote:
Code:
typedef GeometricField<vector, fvPatchField, volMesh> volVectorField 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> { ... } } Quote:
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 | |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
Quote:
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 |
Senior Member
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9 |
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. |
|
July 7, 2020, 09:30 |
|
#18 |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
Quote:
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 |
Senior Member
Carlos Rubio Abujas
Join Date: Jan 2018
Location: Spain
Posts: 127
Rep Power: 9 |
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 |
Member
Rohit George Sebastian
Join Date: May 2017
Posts: 41
Rep Power: 8 |
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? |
|
|
|
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 |