[nandiganavishal]

Hi All,

I would like to solve the following equations using relaxation techniqes. As direct solvers are not available in OpenFoam.

ddt(C1) = K1 * laplcian(C1) + K2 * laplacian (phi)
ddt(C2) = K1 * laplcian(C2) + K2 * laplacian (phi)
K3 * laplacian (phi) = K4 * f(C1) + K5 * f(C2)

where k1,k2,k3,k4 and k5 are some constant values.. f(C1) and f(C2) means the equation is a function of C1 and C2.

I would like to know which tutorial i have to follow to solve these equations. I am fairly new to OpenFoam. K

Kindly do the needful help.

Thanking You

Vishal

[henrik]

I have no experience with PNP, but this might help.

http://groups.google.com/group/sci.math.num-analysis/browse_thread/thread/ec5e96 b20472d308?pli=1

To try relaxation methods is easy with OpenFOAM - Discretise each equation, solve one after another and hope for the best. This is your first step and should yield results if the the time step and/or under-relaxation is small enough.

Here are a few (uneducated) ideas to improve efficiency, but do step 1 first!:

a) Substitute eqn 3 into eqn 1 and eqn 2 to yield 1a and 2a. See whether this is better.
b) Substitution of eqn 2a into 1a is not possible, but you might want to create the matrix and substitute C2=H/A into 1a.

Good luck,

Henrik

[henrik]

There are two forms of relaxation. One is operating and on the field (Phi.relax()) to be used after solving and the other is operating on the matrix (PhiEqn.relax()). The latter is the one you want. What you coded has no effect, it only changes the starting point for the solver.

Also make sure to to have a subDict

relaxationFactors
{
Phi 0.7;
}

in system/fvSolution.

The underrelaxation factor is 1.0 if there is no entry. Make sure underrelaxation is working by rerunning the code (after changing relaxation) and see whether the solver output changes.

Henrik

[nandiganavishal]

Hi Henrik,

Thanks for the response... So Now the code should be something like this right....

fvScalarMatrix PhiEqn
(
fvm::laplacian(Phi) == -alpha * Z1 * C1 - alpha * Z2 * C2
);
PhiEqn.relax();

I would like to know does this relax the equation or the variable Phi.. I meant does this do something like

Phi = Phi * (1-relax) + Phi_old (relax)

I couldn't exactly comprehend what you replied last time..
Can you throw some more light on this,,,

Thanks a lot..

Regards

Vishal

[henrik]

Hi Vishal,

@code: Yes, if you don't forget PhiEqn.solve() later.

@equation: This is exactly what Phi.relax() does - Explicit relaxation.

In implicit relaxation (PhiEqn.relax()) the matrix is modified before solving. The solution is still approached taking it step by step, similar to explicit relaxation, but note that even if you were relaxing only one cell, the effect would be felt throughout the whole domain, in contrast to explicit relaxation. Furthermore, the diagonal dominance of the matrix is increased which is crucial for steady-state calculations where the matrix of a standard transport equation is diagonally equal at best.

@exactly comprehend: Please be more specific and will try and help.

[nandiganavishal]

Hi Henrik,

Thanks for the response.. Now I understood....

I just have a small query regarding implementing some boundary condition.. I would like to incorporate a fixed gradient b.c which is a function of gradient of another variable.. for example like I would like to incorporate

grad(C1) = Z1*C1*grad(Phi)where Z1 - constant

How this type of B.C should be incorporated...
I tried this but I don't know why it is not working...

The code I worte goes like this...

const polyPatchList& patches = mesh.boundaryMesh();

forAll (patches, patchI){

if(C1.boundaryField()[patchI].type()=="fixedGradient"){
fixedGradientFvPatchScalarField& nameOfYourChoice =
refCast<fixedgradientfvpatchscalarfield>(C1.bounda ryField()[patchI]);
scalarField& otherNameOfYourChoice = nameOfYourChoice.gradient();
otherNameOfYourChoice = -C1*Z1*fvc::snGrad(Phi);
}

}
Kindly let me know what is the mistake I am making... As the error I get is something like this

C1 Residual C2 Residual Phi Residual
Begin to solve C11#0 Foam::error::printStack(Foam:stream&) in
"/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in
"/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted: [0xb7f12420]
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&,
Foam::UList<double> const&) in
"/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::operator/(Foam::UList<double> const&, Foam::UList<double>
const&) in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libOpenFOAM.so"
#5 Foam::fixedGradientFvPatchField<double>::evaluate( Foam::Pstream::commsTypes)
in "/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libfiniteVolume.so"
#6 Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh>::GeometricBoundaryField::evaluate() in
"/home/vishal/OpenFOAM/vishal-1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam"
#7 Foam::fvMatrix<double>::solve(Foam::Istream&) in
"/home/vishal/OpenFOAM/OpenFOAM-1.5/lib/linuxGccDPOpt/libfiniteVolume.so"
#8 main in "/home/vishal/OpenFOAM/vishal-1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam"
#9 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#10 Foam::regIOobject::readIfModified() in
"/home/vishal/OpenFOAM/vishal-1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam"
Floating point exception

Kindly do the needful

Thanks

Regards

Vishal

[nandiganavishal]

Dear Henrik,

Thanks for the response...

Here is how I have coded my equations

fvScalarMatrix C1Eqn
(
fvm::ddt(C1)

==
fvm::laplacian(D1,C1) + (D1 * Z1 * fvc::laplacian(Phi) * fvm::Sp(1.0,C1))+ (D1 * Z1 * (fvc::grad(Phi) & fvc::grad(C1)))
);
C1Eqn.relax();


eqnResidualC1 = C1Eqn.solve().initialResidual();
maxResidual = max(eqnResidualC1, maxResidual);

Similarly for C2

fvScalarMatrix C2Eqn
(
fvm::ddt(C2)

==
fvm::laplacian(D2,C2) + (D2 * Z2 * fvc::laplacian(Phi) * fvm::Sp(1.0, C2))+ (D2 * Z2 * (fvc::grad(Phi) & fvc::grad(C2)))

);
C2Eqn.relax();

eqnResidualC2 = C2Eqn.solve().initialResidual();
maxResidual = max(eqnResidualC2, maxResidual);


and for Phi

fvScalarMatrix PhiEqn
(
fvm::laplacian(Phi) == -alpha * Z1 * C1 - alpha * Z2 * C2
);
PhiEqn.relax();

eqnResidualPhi = PhiEqn.solve().initialResidual();
maxResidual = max(eqnResidualPhi, maxResidual);

Buy by using these equations, I don't see any relaxation being carried out even if I change my relaxation factor I don't see any change, further I wanted to solve steady state solution, but I get bizzare results with each iteration... Kindly let me know why is the relaxation not working.. I followed the same way given in the tutorial "BuoyantSimpleRadiation".

Kindly give your suggestions

Thanks a lot

Regards

Vishal

Thanks

Regards

Vishal

[henrik]

Vishal,

I think you should have opened one thread for the complete problem. It also makes it easier for you to find the answers ;-). For interested third parties: BCs are dealt with in this thread

http://www.cfd-online.com/cgi-bin/Op...c=1&post=22484

@code: Looks fine to me. Please post the output of the first two iterations for a run with and a run without underrelaxation as well as the relevant section in system/fvSchemes.

Henrik

[henrik]

Spend an "s" on "relaxationFactor" will make your day.

Henrik

[nandiganavishal]

Dear Henrik,

Thanks a lot for the suggestion, now I see the change..

But still the solution is blowing up with each iteration. I see that the residual is changing but it is not helping in giving me a converged solution. I am solving a steady state problem but still for every time step the error is picking up and I get NAN at the end... Kindly let me know where I am making an error...

Regards

Vishal

[henrik]

Yes, now you are at right at the beginning where the becomes interesting:

a) Why do you have ddt-terms if you are seeking a steady-state solution?

b) Change all relaxation factors to 0.001 and have a look how the fields are evolving. Where is it blowing up? What's causing the trouble?

c) fvc::grad(Phi) & fvc::grad(C2) is nasty. Bring this into a conservative form by using div(s*v) = s*div(v) + grad(s) . v

d) Just cosmetics, but it hurts my eyes:

(D2 * Z2 * fvc::laplacian(Phi) * fvm::Sp(1.0, C2))

is equal to

fvm::Sp(D2 * Z2 * fvc::laplacian(Phi), C2)

Henrik

[nandiganavishal]

Dear Henrik,

Thanks for the reply.

"a) Why do you have ddt-terms if you are seeking a steady-state solution? "

Though I used ddt terms, I mentioned steadyState in my fvschemes, it should have the same effect right.


"b) Change all relaxation factors to 0.001 and have a look how the fields are evolving. Where is it blowing up? What's causing the trouble? "

The reason for the blowing up of the solution is because of the changes in C1 and C2 with each iteration... The C1 and C2 initially are small (inlet and outlet = 0.1 mM).. and with each iteration they are increasing.. Initially they are around 0.5 mM for C1 and C2 is reaching negative values around -0.4 mM (I don't want this to happen), then after few iterations they are just blowing out reaching 10000 mM and all...

"c) fvc::grad(Phi) & fvc::grad(C2) is nasty. Bring this into a conservative form by using div(s*v) = s*div(v) + grad(s) . v "

I am not able to understand how exactly I have to write instead of fvc::grad(Phi) & fvc::grad(C2) in OpenFoam

"d)(D2 * Z2 * fvc::laplacian(Phi) * fvm::Sp(1.0, C2))

is equal to

fvm::Sp(D2 * Z2 * fvc::laplacian(Phi), C2) "

I did change this...

Basically I have my equations as


tau1 = -D1*grad(C1) - D1*Z1*C1*grad(Phi)
tau 2 = -D2*grad(C2) - D2*Z2*C2*grad(Phi)


where C1 C2 and Phi are concentrations and potential respectively

D1,D2,Z1,Z2 are constants
This is the set of coupled equations I am solving...

dC1/dt = - div (tau1)
dC2/dt = - div (tau2)

For Steady State, the left hand side is 0.

laplacian (phi) = - alpha * Z1*C1 - alpha*Z2 * C2

where alpha is a constant

I have expanded the div(tau1) and similarly div(tau2) and wrote everything as laplacian and gradient terms.. So I don't have any divergence terms as such hence even if I get negative values, I am not able to rectify them.

Kindly give your suggestions as to how I should go about solving these issues and whether my equations are correct for solving this kind of problem.

Thanks

Regards

Vishal

[henrik]

a) Okay
b) Okay, I understand that the solution has bounds (C1 > 0 C2 > 0) which are violated
c) One hint: s = Phi and v = grad(C2).

One more thing: Do you understand what the conservative form of a transport equation is? If not, look it up. Any textbook on FVM will do. There is no chance you will solve this problem if you do not understand the basics - especially in the light of b)

d) Try suggestion a) from my first quote. This gives you two conservative equations

[henrik]

Vishal,

@I) I tried to suggest the following:

Replace

0 = K1 * laplacian(C1) + K2 * laplacian (phi)
0 = K1 * laplacian(C2) + K2 * laplacian (phi)
K3 * laplacian (phi) = beta(C1, C2)
beta (C1,C2) = K4 * f(C1) + K5 * f(C2)

by

0 = K1 * laplacian(C1) + K2/K3*beta(C1,C2)
0 = K1 * laplacian(C2) + K2/K3*beta(C1,C2)

This should be easier to solve because you have fewer equations and they are conservative. But you need to be careful when discretising the source to ensure boundedness at zero.

What do f(C1) and f(C2) look like?

Henrik

[nandiganavishal]

Dear Henrik,

The equations which I have quoted in the beginning of the forum is not the equations I am looking to solve... The exact equations which I would like to solve are

div {(D1*grad(C1) + D1*Z1*C1*grad(Phi))} = 0 -(1)
div {(D2*grad(C2) + D2*Z2*C2*grad(Phi))} = 0 -(2)
laplacian (phi) = - alpha * Z1*C1 - alpha*Z2 * C2 - (3)

where C1 C2 and Phi are concentrations and potential respectively and D1,D2,Z1,Z2, alpha are constants

Sorry for the confusion. I am trying to solve these equations,

Initially I have tried to expand each term on the LHS of equation 1/2 and wrote them as laplacian terms.. When I solve them like that then after every iteration I see a change in my C1 and C2 and eventually blowing out my solution...

Hope I have made myself clear now.. I am again sorry for the confusions..

Thanks

Regards

Vishal

[dmoroian]

Shouldn't this be written as:

fvScalarMatrix Eqn1(
fvm::div(D1*fvc::grad(C1)+D1*Z1*C1*fvc::grad(Phi))
);
fvScalarMatrix Eqn2(
fvm::div(D2*fvc::grad(C2)+D2*Z2*C2*fvc::grad(Phi))
);
fvScalarMatrix Eqn3(
fvm::laplacian(Phi)+alpha*Z1*C1+alpha*Z2*C2
);

Eqn1.relax();
Eqn2.relax();
Eqn3.relax();
Eqn1.solve();
Eqn2.solve();
Eqn3.solve();

Dragos

[dmoroian]

Your entry div(((Z1*C1)*fvc::grad(Phi))) Gauss upwind; in the dictionary is wrong. Write instead div(((Z1*C1)*grad(Phi))) Gauss upwind;

On the other hand, I don't know why it didn't work in the previous formulation.

Dragos

[dmoroian]

It seems there is an extra space in div(((Z1*C1)*grad(Ph i)))

[nandiganavishal]

Dear Dragos,

No I just checked it again... There is no space

fvm::laplacian(C1) + fvc::div(((Z1*C1)*fvc::grad(Phi)))

divSchemes
{
default none;
div(((Z1*C1)*grad(Phi))) Gauss upwind;
div(((Z2*C2)*grad(Phi))) Gauss upwind;
}

Calculating concentration distribution

Time = 0.2

C1 Residual C2 Residual Phi Residual


attempt to read beyond EOF

file: /home/vishal/OpenFOAM/OpenFOAM-1.5/tutorials/coupledFoam_2dSteady_majumdar_geome try/coupledFoam2dSteady_majumdar_geometry/system/fvSchemes::div(((Z1*C1)*grad(Ph i))) at line 31.

From function ITstream::read(token& t)
in file db/IOstreams/Tstreams/ITread.C at line 64.

FOAM exiting

But it is displayed like this... when I post it..

Regards

Vishal

[dmoroian]

Hello guys,
Can anyone tell me why the following expression gives me a compiler syntax error?

fvScalarMatrix Eqn1(
fvm::div(D1*fvc::grad(C1)+D1*Z1*C1*fvc::grad(Phi))
);

Dragos