[nep101]

Hi everybody,

im new to Openfoam so i get stuck. I try solve this (dimensionless) Equ in a Simple Loop:

U-Ra*T*k=-grad(p)

were U is the velocity vectorfield, Ra is the Rayleighnumber, T is the temperature, k is a vector (Ra*T*k is a volScalarField) and p is the Pressure.

I want to creat a Matrix UEqn with
{
U-RaTk
}

But it dont work.

Code:
volScalarField V("V", U & mesh.C());
fvScalarMatrix UEqn
(
V-RaTk == fvc::grad(p) //Ratk is a VolScalarField
);



UEqn.solve();

Anybody got a hint?
Werner

[TomB]

Hi Werner,

I am working with Darcy's eqns in OpenFoam as well.

how did you define the k vector ? Can you multiply it by a scalar the way you have done ?

Is it okay to have a UEqn with the double equals in it. I think not.

I Look fwd to working with you !

Originally Posted by nep101

Hi everybody,

im new to Openfoam so i get stuck. I try solve this (dimensionless) Equ in a Simple Loop:

U-Ra*T*k=-grad(p)

were U is the velocity vectorfield, Ra is the Rayleighnumber, T is the temperature, k is a vector (Ra*T*k is a volScalarField) and p is the Pressure.

I want to creat a Matrix UEqn with
{
U-RaTk
}

But it dont work.

Code:
volScalarField V("V", U & mesh.C());
fvScalarMatrix UEqn
(
V-RaTk == fvc::grad(p) //Ratk is a VolScalarField
);



UEqn.solve();

Anybody got a hint?
Werner

[Cyp]

If you want to code a Darcy equation in OpenFOAM (or in whatever codes) it is very easy. In fact, it is much easier than a Navier-Stokes equation with PISO/SIMPLE algorithms...

The complete system is defined by :
1°) a continuity equation. Most of the time it is div (U) = 0
2°) a momentum equation, which is reduced to the Darcy's law in case of porous media and small Reynold's number. U = - K/mu*grad(P)

As in most of the (U,p) solvers, you must solve a pressure equation. It is derived by inserting the Darcy's law within the continuity equation :

div(U)=0=div(-K/mu*grad(p))=laplacian(-K/mu,p)

Consequently, in OpenFOAM you will code :

solve
(
        fvm::laplacian(K/mu,p)
);

U = -K/mu*grad(p);
U.correctBoundaryConditions();

In this example, the permeability K can be defined as :
- a dimensionedScalar (in case of isotropic homogeneous porous medium)
- a dimensionedTensor (in case of anisotropic homogeneous porous medium)
- a volScalarField (in case of isotropic heterogeneous porous medium or/and in case of velocy-dependance (Forchheimer correction for instance))
- a volTensorField (in case of anisotropic heterogeneous porous medium or/and in case of velocy-dependance)

If you have a source term in the continuity equation, for example :
div(U) = -q
Therefore the pressure equation becomes :
laplacian(-K/mu,p) = -q

I do not know your equation with your Ra*T*k term, but the method above is always available.

Regards,
Cyp

[Cyp]

This answer is the same for the two others topics you have opened...

[TomB]

Hi CyP

thanks very much for your guidance. I'll put some thought into it and reply again later.

The Ra*T*k term is there as we are looking at flow in porous media that is heated (from below) so the convection is the interesting thing.

Ra is rayliegh number
T is the temperature (im unsure if it will be Tref, or a 'local' Temp)
and k is just the unit vector in the vertical direction.

Cheers

[TomB]

Cyp,
I compiled a darcySimpleFoam based on your guidance.
i attached the solver code.

I based it on the simpleFoam as the name suggests, so it still has some of the simpleFoam headers etc.

I need to test it, see if it makes sense. say in a box of porous media with a pressure difference equal to rho*g*h from top to bottom,

From the way I have set up the pEqn, it doesnt seem to take any inputs, so how does it 'know' what the pressure difference is.. ?

Also, I understood from your post that I dont need a PISO/SIMPLE method to solver these equations.

What I do want though is a transient solver? Does this complicate things, or will the same code work?

Cheers

[Cyp]

Hi TomB/neph101,

Why are you persisting with PISO/SIMPLE loop ? If you want to base your darcy solver on an existing solver, you should start from laplacianFoam.

I suggest you to work in several step:
1°) Fisrt you code and compile a simple Darcy solver (as the one I wrote below)

1°bis) Add a simple scalar transport equation (see example in scalarTransportFoam)

2°) Then, think to the way to add the temperature-dependent term in your solver. if U = Ra*k*T -grad(p) therefore, div(U) = Ra*k&grad(T) - laplacian(p) = 0 ...

3°) Consider transient effect by starting from a compressible continuity equation. Most of the time we transform the density of the time-derivative to a pressure via the perfect-gas law. There is probably other way. From example, you can solve an steadyState darcy at each time step while you temperature is solved in transient manner..

Regards,
Cyp

[nep101]

Hi Cyp, TomB!

First thanks for your help and tips.
But i dont get it.

My Equations for the thermal convection in a porous cylinder (dimensionless) is this:
ddt(T)+U&grad(T)-laplace(T)=0
Make it a little easier for testing:
U&grad(T)-laplace(T)=0,
div(U)=0 and
U+grad(p)-Ra*T*k=0

The Way i want to solve it is this:
1) solve the Tempeqn.
2) solve the continuity equation (this is what not work)
3) solve U

ad 1:
this is easy (i think)

solve
(
(U&fvc::grad(T))
- fvm::laplacian(T)
);

ad 2:
with div(U)=0=div(Ra*k*T)-laplacian(p)
or like Cyp already wrote
div(U) = Ra*k&grad(T) - laplacian(p) = 0

fvScalarMatrix pEqn
(
fvm::laplacian(p) == fvc::div(Ra*k*T)
);
pEqn.solve();

ad 3:

U=Ra*k*T-fvc::grad(p);

But something with the continuity equation went wrong. I get a pressure arround e+16 . And Openfoam need 1000 interations in this step.

Cheers

[TomB]

Hi Nep,

Why solve the Temp Eqn first? From what I understand, the scheme would be to solve first the pEqn, then UEqn and then TEqn.

I have been working on the scheme today, guided by Cyp's advice.

I have attached my solver.

I have issues with it,
1) I need to add the p_rgh (by that I mean pressure should increase with depth)
2) I need to add a TEqn
3) The UEqn needs the Ra*T term, which I still have to make a dictionary for and implement in the UEqn

I also want to do this all non-dimensionally, which is adding another level of confusion.

Cheers

Tom

[TomB]

Hi Cyp,

Thanks for your input. I was confused by the using the piso/simple loop. As you guessed, it was just a convenience to base my solver on an existing one. I'll re-base it on laplacianFoam and then follow your advice.

In regards to the temperature, would you recommend taking guidance from bouyantBoussinesq(Simple or Pimple)Foam ?

I do want to include the boussinesq approx, and the next thing I need to figure out is how to include the change in pressure with depth.

Cheers

Tom

[Cyp]

Your problem is very easy to code ! much easier than Navier-Stokes equation. In fact, coding a Darcy's law can help you to understand the (U,p) solver for Navier-Stokes equations...

For the moment, we forget the temperature scalar transport equation. So your system is defined by
1°) a continuity equation
2°) a momentum equation

You have to notice that the unknow variable in the momentum equation is the velocity. At this point, there is no differential equation that govern the pressure field. However, the knowledge of this pressure field is very important since it appears as the main source term in the momentum equation.

You must also notice that the continuity equation is never directly solve in a (U,p) solver. In fact, we use this equation as a constraint to built the equation that govern p. The pressure equation laplacian(-K/mu,p) will always satisfy the continuity equation, even if we do not directly solve div(U) = 0

The first step of your development is to code (and understand how to code) a simple Darcy's law ! Then you could the pressure effect.

Another remark : the divergence operator in OpenFOAM require a velocity flux (something like div(phi) where phi = linearInterpolate(U) & mesh.Sf() )

I can help you, but you must proceed step by step.

@++
Cyp

[nep101]

Hi Tom,

i think it is not so important in wich order you solve the equatios. My Problem is that something is going wrong with my continuity equation. After every interation my pressure is increases until Openfoam interrupts (at a point the pressure goes inf).

solve(fvm::laplacian(p) - fvc::div(Ra*k*T));

My Equat. are also nondimensional. What you must do is set dimensionSet to 0 in opt/openfoamXXX/etc/controlDict. Or you will get problems with the Units.

Adding rgh to the pressure is not so difficult i think look at the creatFileds in boussinesqSimpleFoam:

Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());

volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rhok*gh
);

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PIMPLE"),
pRefCell,
pRefValue
);

[nep101]

Hi Cyp,

Thanks for your input. But I really dont now how to code it.

Ok forget about the Temp so i got the 2 Equations:
div(U)=0 and U+grad(p)-Ra*k*T=0.
Ra is the Rayleighnumber and k is a Vector. I start with U=0 and heat a side of my object, so I can say my U is known and =0. So I must solve the pressure and satisfy the continuity equation.

I cannot solve div(U) the way i tried before.
Maybe with Rhie-Chow interpolation?
Like:
Ueqn = U - Ra*k*T
AU = UEqn().A();
U = UEqn().H()/AU
create phi
solve(fvm::laplacian(1.0/AU, p) == fvc::div(phi))
U=-grad(p)+Ra*k*T

What do you think Cyp or is this totally wrong????

[Cyp]

Forget the Temperature, so your velocity is U = -grad p (according to your dimensionless analysis).

Forget also the predictor/corrector algorithm. Just think about mathematic with pen and paper.

The first step is to code :
div(U) = 0
U = -grad(p)

I mentioned all the stages in my previous posts. At this point you may also forget the "phi" part.

let me know of your progression.

@++
Cyp

[TomB]

Hi Nep
Thanks.
Can you tell me how you set up the k vector?
Cheers

Tom

[nep101]

Hi Cyp,

solve
(
fvm::laplacian(p)
);

U=-fvc::grad(p);

I dont get it.
You said U=-grad(p) so with div(U)=0 you get laplacian(p)=0 and U=-grad(p).
But in my case i got from div(U)=0 laplacian(p)=div(Ra*k*T)=Ra*d(T)/dz.

Hi Tom,

I defined k as a volVectorField [0 0 1]

Cheers
Werner

[Cyp]

Indeed in your case you have another term in your equation. However it is easy to add it when you have a pressure/velocity solver that works fine !

Now that you have your velocity/pressure solver (if we neglect the temperature term) you can add the advection-diffusion temperature equation :

solve
(
    fvm::ddt(T) + fvm::div(phi,T) - fvm::laplacian(DT,T)
);

in OpenFOAM, "phi" is the velocity flux, defined as

phi = linearInterpolate(U) & mesh.Sf()

You must add this line just after the calculation of U

When you will have check the validity of your solver, you can add the temperature correction in the calculation of U and p.

Ra*k&fvc::grad(T)

Don't forget that you always have to work step by step when you code a solver. It is too difficult (and source of error) to program all your problem in one time...

@++
Cyp

[nep101]

Hi Cyp,

I realy thank you for your help, but I think I give up. It is allwasy the same and I dont know were the problem is.

My Code (another attempt):

volVectorField kT("kt",Ra*k*T);
volScalarField divT = fvc::grad(kT)().component(vector::Z);
solve
(
fvm::laplacian(p)- divT
);
U=-fvc::grad(p)+Ra*k*T;
phi = linearInterpolate(U) & mesh.Sf();
solve
(
fvm::ddt(T) + fvm::div(phi,T) - fvm::laplacian(T)
);

But it is alltime the same. After one or two interations the pressure explode (e+12) and thats it.

I cant compute the pressure right. So once again Cyp thanks for your time and your help.

Cheers
Werner

[Cyp]

Originally Posted by nep101

Hi Cyp,

I realy thank you for your help, but I think I give up. It is allwasy the same and I dont know were the problem is.

My Code (another attempt):


But it is alltime the same. After one or two interations the pressure explode (e+12) and thats it.

I cant compute the pressure right. So once again Cyp thanks for your time and your help.

Cheers
Werner

What are the boundary condition you applied ?

And I am not sure about the way you code divT.. I think you should use the snippet I previously proposed.

And you should also use U.correctBoundaryConditions() after the calculation of U since U is evaluated from another field

[nep101]

Hi Cyp,
i use a simple cube for testing.

T:

internalField uniform 450;

boundaryField
{

seite
{
type zeroGradient;
}
fab
{
type zeroGradient;
}


unten
{
type fixedValue;
value uniform 500;
}
oben
{
type fixedValue;
value uniform 450;
}

p:

internalField uniform 0;

boundaryField
{

seite
{
type zeroGradient;
}
fab
{
type zeroGradient;
}


unten
{
type zeroGradient;
}
oben
{
type zeroGradient;
}
}

and U:

internalField uniform (0 0 0);

boundaryField
{

seite
{
type fixedValue;
value uniform (0 0 0);
}
fab
{
type fixedValue;
value uniform (0 0 0);
}


unten
{
type fixedValue;
value uniform (0 0 0);
}
oben
{
type fixedValue;
value uniform (0 0 0);
}
}

with oben is the top and unten is the bottom, fab and seite are the side of the cube.

Cheers
Werner