Creating a solver for the nondimensionalized Navier-Stokes equation
 

Hi everyone,

I have recently started using OpenFoam.
I want to resolve the nondimensionalized Navier-Stakes equation (in which appears the Reynolds number).

My equation is ( you can also find it attached ) : Re . (U . div) U + grad(meanP) - laplacian(U) = - grad(fluctP)


http://www.cfd-online.com/Forums/dat...AASUVORK5CYII= meanP is a unit vector which is constant and P=meanP+fluctP
Re is the reynolds number



I want to resolve my equation for U and fluctP


For this, I have modified the IcoFoam solver as follows, but it seems that there is an error ( may be with the '' meanP '' )



fvVectorMatrix UEqn
(
fvm::ddt(U)
+ Re*(fvm::div(phi, U))
- fvm::laplacian(U)
+ magSqr(meanP)
);

solve(UEqn == -fvc::grad(fluctP));


I will be thankfull If someone can help me to create this solver or share a web site which treated the implementation of the nondimensionalized NS equation.


Best regards

Mehrez

Hi

Can someone help me please.

Thank you

please find attached the OpenFoam files of the solver (I have proceeded by modifying the IcoFoam solver to myIcoFoamB solver).

Hi,

there are many errors:

First: laplacian(U) is not a class. If you use laplacian you have to a coefficent like:

You set Re as reynoldsnumber. But the solver do not know what "Re" is. Its a none declared variable. The same like magSqr(meanP)

OF do not know what meanP is.
So you have to define

Re
meanP

and the other variables you implemented!

The following example is working:

Hi Mehrez,

First it would be better if you modify the Make/files as in the user guide chapter 3 to be:
When I try to compile your solver (initially) two errors occur from different places. This points out that you have made a punch of modifications then decided to compile ... Not a good idea .. It is always a good practice to compile regularly after small code changes to keep track of code bugs early on. It is a little bit tedious but very efficient. Back to the errors:

Then the compiler gives you some ideas about what could be the problem. Anyway the code that caused this error is:
Essentially this is not how you declare a volVectorField. The compiler is complaining that there is no lookup function for the volVectorField and if there were one (which isn't the case) you are doing it wrong and hence the error. volVectorField are defined like the p if you'll read as user input but if you need to calculate:
The second error:
Comes from the eq:
Maybe it will be solved if you fix the first one or not ...

Best regards,
Hisham

Hi Hisham, the laplacian class is wrong like I posted befor.
I just had a very short look into the code and dont see the declaration of the other variables :)

Well nice replay!

Hi guys

Thank you for your help.

I think that it is better if you take a look to my equation (which I have attached in my first message). You can see that I don't have : nu . laplacian (U) but I just have : laplacian (U)

Concerning the declaration of the variables : I have declared " Re " and " gradP " in the file createFields.H

dimensionedScalar Re
(
transportProperties.lookup("Re")
);

volVectorField gradP
(
transportProperties.lookup("gradP")
);

Oh sorry,

the laplacian(U) class existis! So its working :)

Hi Hisham

Thank you for your answer.

I think that I don't have to declare " gradP " like this :

volVectorField gradP ( IOobject ( "gradP", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), fvc::grad(p) );

As I said before, the '' gradP '' is a unit vector which is constant and which I put it like an input (so this is why OF don't have to compute it).
I have just to declare it like a vector which I will specify the value in constantProperties file.

Hmmm ... you calculate your gradP field with grad(p) or?
But how whould you get the information of the "vector" couse the pressure is just a scalar without a direction.

Hi Tobi

As you can see in my equation gradP is the gradient of the nondimensionalized mean pressure. I have it like an input (I will study my flow in function of given values of Re and gradP)

I think that the problem is how to define this vector (gradP) in the equation

I think that it is more clear like this

Equation : Re . (U . div) U + gradP - laplacian(U) = - grad(fluctP)

OF syntax :

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ Re*(fvm::div(phi, U))
- fvm::laplacian(U)
+ magSqr(gradP)
);

solve(UEqn == -fvc::grad(fluctP));


Input :
Re : scalar
gradP : vector

output (after computation)
U : vector field
fluctP : vector field

If gradP is a constant vector, you should not store it in a volVectorField, but in a dimensionedVector from the transportproperties.

Two more things:
- Are you sure Reynolds should appear like you wrote it. As far as I know, in the dimensionless NS equations, you have 1/Re in the diffusive terms.
- You're now adding magSqr(gradP), which is a scalar, to a vector equation. This can not be correct.

Dear Bernhard,

I'm sure that I have a right equation (attached here)

Thank you very much for your help. Actually I can compile my solver with your corrections.

fvVectorMatrix UEqn
(
fvm::ddt(U)
Re*(fvm::div(phi, U))
- fvm::laplacian(U)
+ gradP
);

and :

dimensionedVector gradP
(
transportProperties.lookup("gradP")
);


To test this solver, I've taken the attached example "cavity" after putting U, p, Re, and gradP in a dimensionless form. I've done the same with the domain (I've commented the line " conversion to meters " ).

Executing "cavity" with "myIcoFoamB" solver, OpenFoam returns an error message on the dimensions:

-> FOAM FATAL ERROR:
incompatible dimensions for operation
[U [0 -1 0 0 0 0 0]] - [U [0 -2 0 0 0 0 0]]

I think it comes from the operators (grad gives : L ^ -1 and Laplacian gives : L ^ -2). Is there a way to put the components (x, y, z) of the domain in a dimensionless form in order to get dimensionless values ​​by applying one of the operators?


Thank you very much for your help.

Best regards.

Mehrez

Complicated one, you can choose to ignore the dimension-checking. I think you can do this in the controlDict you find in ..../OpenFOAM-2.1.x/etc/ .
Maybe you can put it in your own controlDict as well. I never tested this, but you can look here.

Hi,

you can put a dimensionSet on your p vector that the dimension check is working!

Tobi

Hi guys,

Again thanks for your precious help.

@ Bernhard : I have found the file " controlDict " , can you please precise me how to proceed.

@ Tobi : I've tried what you said but it doesn't work. it returns each time a different dimensions error ! I can't understand

Mehrez

Post your dimenionSet and the output error.

Re Re [ 0 0 0 0 0 0 0 ] 0.01;
gradP gradP [ 0 0 0 0 0 0 0 ] (0.9 0.2 0);
U : dimensions [0 0 0 0 0 0 0];
p : dimensions [0 0 0 0 0 0 0];

In the blockMeshDict, I have commented the following line :
//convertToMeters 0.1;

and this is what I get :






Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi


Starting time loop

Time = 0.005

Courant Number mean: 0 max: 0


--> FOAM FATAL ERROR:
incompatible dimensions for operation
[U[0 -1 0 0 0 0 0] ] - [U[0 -2 0 0 0 0 0] ]

From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
in file /opt/openfoam211/src/finiteVolume/lnInclude/fvMatrix.C at line 1316.

FOAM aborting

#0 Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#1 Foam::error::abort() in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#3
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#4
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#5 __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#6
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
Aborted (core dumped)
ubuntu@ubuntu-VirtualBox:~/OpenFOAM/mehrez-2.1.1/test/cavity$

What is that:

U : dimensions [0 0 0 0 0 0 0];
p : dimensions [0 0 0 0 0 0 0];


??? I dont understand the two lines.
Set

pGrad pGrad [0 -1 0 0 0 0 0]

dimensions [0 0 0 0 0 0 0];

This line represent the dimension of U and p. I have this line in the 2 files cavity/0/p and cavity/0/U

Dear Tobi,

Here is what I got by setting pGrad pGrad [0 -1 0 0 0 0 0]







Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi


Starting time loop

Time = 0.005

Courant Number mean: 0 max: 0


--> FOAM FATAL ERROR:
incompatible dimensions for operation
[U[0 -1 0 0 0 0 0] ] - [U[0 -2 0 0 0 0 0] ]

From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
in file /opt/openfoam211/src/finiteVolume/lnInclude/fvMatrix.C at line 1316.

FOAM aborting

#0 Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#1 Foam::error::abort() in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#3
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#4
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#5 __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#6
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
Aborted (core dumped)
ubuntu@ubuntu-VirtualBox:~/OpenFOAM/mehrez-2.1.1/test/cavity$

Even if I do this, it doesn't work :

0/p
dimensions [0 -1 0 0 0 0 0];

0/U
dimensions [0 0 0 0 0 0 0];

transportProperties
Re Re [ 0 -1 0 0 0 0 0 ] 0.01;
gradP gradP [ 0 -2 0 0 0 0 0 ] (0.9 0.2 0);

This is what I got :






Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi


Starting time loop

Time = 0.005

Courant Number mean: 0 max: 0
DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 6.17543e-06, No Iterations 17
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 3.68856e-06, No Iterations 14


--> FOAM FATAL ERROR:
LHS and RHS of + have different dimensions
dimensions : [0 2 0 0 0 0 0] + [0 4 -1 0 0 0 0]


From function operator+(const dimensionSet&, const dimensionSet&)
in file dimensionSet/dimensionSet.C at line 514.

FOAM aborting

#0 Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#1 Foam::error::abort() in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2 Foam::operator+(Foam::dimensionSet const&, Foam::dimensionSet const&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#3
at myIcoFoamB.C:0
#4
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
#5 __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#6
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/myIcoFoamB"
Aborted (core dumped)
ubuntu@ubuntu-VirtualBox:~/OpenFOAM/mehrez-2.1.1/test/cavity$

Any idea ? does someone know how to ignore the dimension-checking in the "controlDict" file ?

I'm stuck...

Oh its not my day today, sorry.
Now I understand what you are doing.

Hmmm I have to try it by my self but i have no time. It should be possible.

If you do
cd $FOAM_INST_DIR/OpenFOAM-2.1.0/etc

You can go to the line with dimensionSet (line 400 for me), and change the switch to 0. I assume this does not check the dimensions, but I am not sure about it.

If you don't have privileges to change that controlDict, you can set it in a .OpenFOAM/2.1.0/ folder in your home-directory.

Good luck, I hope this works.

Dear both,

It seems that it is working !

Thank you for your precious help.

I hope that this discussion will help other beginners.

Thank you.

Best regards.

Mehrez.

Dear Mehrez

It is nice that your solver is working now. However, I think you better track the dimensions of your equation ... It is a very nice feature to override. Especially, if you have problems with verification of your solver.

Regards,
Hisham

Dear Hisham,

If you look just a little above you will see that I tried to multiply the terms of my equation with coefficients of different dimensions to get a dimensionless form of my equation but it did not work and I do not know why ...

Otherwise I have another question, how to have the steady state? (is not more convenient if I modify another solver like SimpleFoam because removing the line ddt (U) on IcoFoam gives a false result)

thank you

Mehrez

You can edit the system/fvSchemes from:
to:
But that will just automatically ignore the ddt term. So it has the same shortcoming.

Regards,
Hisham

I just remembered. You can also try to base your solver on the pimpleFoam solver, which uses a merge of PISO and SIMPLE algorithms.

Hi Hisham,

I've found the following text on the OpenFoam tutorial :

" When specifying a time scheme it must be noted that an application designed for transient problems will not necessarily run as steady-state and visa versa. For example the solution will not converge if steadyState is specified when running icoFoam, the transient, laminar incompressible flow code; rather, simpleFoam should be used for steady-state, incompressible flow. "

So, to implement a solver that resolves the following equation (steady-state) :

fvVectorMatrix UEqn
(
// fvm::ddt(U)
Re*(fvm::div(phi, U))
- fvm::laplacian(U)
+ gradP
);

solve(UEqn == -fvc::grad(p));

I need to modify the simpleFoam solver rather than icoFoam.


Looking at the simpleFoam solver, I've found 2 files pEqn.H and UEqn.H but I don't know how to modify the equations... can you please tell me what I need to delete (like the turbulence terms).

Thank you.

Best regards.

Mehrez

OpenFOAM solvers use three major algorithms:
- SIMPLE: for steady state (e.g. simpleFoam)
- PISO: for transient problems (e.g. icoFoam and pisoFoam)
- PIMPLE: for steady state and transient because it is a merge of SIMPLE and PISO (e.g. pimpleFoam)

So as I mentioned earlier you need to base your solver on pimpleFoam. When you compare simpleFoam you need to compare it with pisoFoam because both has turbulence modelling while icoFoam has no turbulence modelling capabilities enabled.

Regards,
Hisham

Dear Hisham,

Thank you for your help.

In pimpleFoam there are 2 files to implement my equation : UEqn.H and pEqn.H , so it is different from icoFoam where I had one file to implement the equation.

Can I do some change on the icoFoam solver to be able to solve with pimple algorithm rather than piso ? If not, in which file do I have to put my equation in the pimpleFoam solver ?

Thank you so much for your help and excuse me for these questions.

best regards.

Mehrez

Hi,

did you see http://openfoam-extend.sourceforge.n...m/training.htm
take a look at the tutorials "New Developer"

http://switch.dl.sourceforge.net/pro...he_slides2.pdf
=>start @Slide10 to 19 Copy and rename scalarTransportFoam

Thank you elvis

Hi,
I am having the similar problem, please help me to deal with.
Thanks

How to turn off the dimension checking in OpenFOAM
 

Dear Bernhard,
I have copied the controlDict out into the case folder, and changed the dimensionSet from 1 to 0, since I have no permission to change the original local file. But still an error occurs:
I have set the dimensions in 0/p, 0/U and the transProperties to be [0 0 0 0 0 0 0].
It seems that the newly created solver does not listen to the modified controDict file.
I am really appreciated for any suggestions.
Thank you in advance.
Best Regards,
Bill

Why create a non-dimensional solver for NS? Use the existing ones!
 

Ok, let me leave you some comments about what you are trying to do:

(1) It seems to me that you are trying to solve the equations directly! But as you are currently implementing it will certainly diverge (even if you make it compile) because you are not giving an equation for p (the system is ill-posed), simply as that. Refer to a basic text of CFD to get acquainted with different solution algorithms. The book of Peyret is a nice introduction to the subject.

(2) You don't need this at all! If you choose the natural frequency of the system to be U/L then the dimensional and non-dimensional version of the NS eqns are equivalent, except that the viscous diffusivity multiplying the laplacian term becomes the inverse of the Reynolds number. In other words, when asked for nu in 'transportProperties' you insert 1/Re. No need to play with the dimensions either! PISO/SIMPLE doesn't care (should it?), but FOAM does dimension checking as a convenience for the user and the developer.

TL;DR: Don't bother implementing a new solver, just use ico-/simple-/pisoFoam considering nu equal to the inverse of Reynolds number, unless the natural freq. of the system is given at the boundaries.

the detailed descripition about the problem
 

Dear Santiago,
Thank you for the comments, I have described the problem detailly here, with the related equations and the decomposed pressure. https://www.cfd-online.com/Forums/op...tml#post628609

This thread may help to explain my problem.
Thank you for the help.
Cheers,
Bill:)