CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

How to solve using relaxation techniques

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

Like Tree5Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   November 25, 2008, 16:17
Default Hi All, I would like to sol
  #1
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   November 29, 2008, 15:59
Default I have no experience with PNP,
  #2
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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 is offline   Reply With Quote

Old   February 2, 2009, 08:03
Default There are two forms of relaxat
  #3
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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
henrik is offline   Reply With Quote

Old   February 3, 2009, 00:22
Default Hi Henrik, Thanks for the r
  #4
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   February 3, 2009, 12:56
Default Hi Vishal, @code: Yes, if y
  #5
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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.
henrik is offline   Reply With Quote

Old   February 4, 2009, 17:19
Default Hi Henrik, Thanks for the r
  #6
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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 is offline   Reply With Quote

Old   February 6, 2009, 02:38
Default Dear Henrik, Thanks for the
  #7
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   February 6, 2009, 05:10
Default Vishal, I think you should
  #8
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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 is offline   Reply With Quote

Old   February 6, 2009, 14:28
Default Spend an "s" on "relaxationFac
  #9
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Spend an "s" on "relaxationFactor" will make your day.

Henrik
henrik is offline   Reply With Quote

Old   February 6, 2009, 15:10
Default Dear Henrik, Thanks a lot f
  #10
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   February 6, 2009, 17:19
Default Yes, now you are at right at t
  #11
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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
mm.abdollahzadeh likes this.
henrik is offline   Reply With Quote

Old   February 6, 2009, 19:53
Default Dear Henrik, Thanks for the
  #12
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   February 7, 2009, 06:37
Default a) Okay b) Okay, I understand
  #13
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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 is offline   Reply With Quote

Old   February 7, 2009, 16:07
Default Vishal, @I) I tried to sugg
  #14
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
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
henrik is offline   Reply With Quote

Old   February 7, 2009, 16:33
Default Dear Henrik, The equations
  #15
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   February 8, 2009, 06:22
Default Shouldn't this be written as:
  #16
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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 is offline   Reply With Quote

Old   February 8, 2009, 11:26
Default Your entry div(((Z1*C1)*fvc::
  #17
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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 is offline   Reply With Quote

Old   February 8, 2009, 11:51
Default It seems there is an extra spa
  #18
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
It seems there is an extra space in div(((Z1*C1)*grad(Ph i)))
dmoroian is offline   Reply With Quote

Old   February 8, 2009, 12:05
Default Dear Dragos, No I just chec
  #19
Senior Member
 
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 206
Rep Power: 9
nandiganavishal is on a distinguished road
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
nandiganavishal is offline   Reply With Quote

Old   February 9, 2009, 02:30
Default Hello guys, Can anyone tell m
  #20
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
dmoroian is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
CFD simulation techniques samar FLUENT 1 November 9, 2007 05:15
UDF - Debugging techniques PK FLUENT 2 February 6, 2007 17:44
TUI- solve/set/relaxation method Ugur FLUENT 0 September 21, 2004 10:05
Stabilization techniques? Mihail Galabov Main CFD Forum 1 April 23, 2004 11:57
Meshing Techniques phil Main CFD Forum 3 November 26, 2002 08:04


All times are GMT -4. The time now is 03:41.