
[Sponsors] 
November 25, 2008, 16:17 
Hi All,
I would like to sol

#1 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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 

November 29, 2008, 15:59 
I have no experience with PNP,

#2 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
I have no experience with PNP, but this might help.
http://groups.google.com/group/sci.math.numanalysis/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 underrelaxation 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 

February 2, 2009, 08:03 
There are two forms of relaxat

#3 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
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 

February 3, 2009, 00:22 
Hi Henrik,
Thanks for the r

#4 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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 * (1relax) + 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 

February 3, 2009, 12:56 
Hi Vishal,
@code: Yes, if y

#5 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
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 steadystate 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. 

February 4, 2009, 17:19 
Hi Henrik,
Thanks for the r

#6 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #1 Foam::sigFpe::sigFpeHandler(int) in "/home/vishal/OpenFOAM/OpenFOAM1.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/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #4 Foam::operator/(Foam::UList<double> const&, Foam::UList<double> const&) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libOpenFOAM.so" #5 Foam::fixedGradientFvPatchField<double>::evaluate( Foam::Pstream::commsTypes) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libfiniteVolume.so" #6 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::evaluate() in "/home/vishal/OpenFOAM/vishal1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" #7 Foam::fvMatrix<double>::solve(Foam::Istream&) in "/home/vishal/OpenFOAM/OpenFOAM1.5/lib/linuxGccDPOpt/libfiniteVolume.so" #8 main in "/home/vishal/OpenFOAM/vishal1.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/vishal1.5/applications/bin/linuxGccDPOpt/electroosmoticF oam" Floating point exception Kindly do the needful Thanks Regards Vishal 

February 6, 2009, 02:38 
Dear Henrik,
Thanks for the

#7 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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 

February 6, 2009, 05:10 
Vishal,
I think you should

#8 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
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.cfdonline.com/cgibin/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 

February 6, 2009, 14:28 
Spend an "s" on "relaxationFac

#9 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
Spend an "s" on "relaxationFactor" will make your day.
Henrik 

February 6, 2009, 15:10 
Dear Henrik,
Thanks a lot f

#10 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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 

February 6, 2009, 17:19 
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 
Yes, now you are at right at the beginning where the becomes interesting:
a) Why do you have ddtterms if you are seeking a steadystate 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 

February 6, 2009, 19:53 
Dear Henrik,
Thanks for the

#12 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
Dear Henrik,
Thanks for the reply. "a) Why do you have ddtterms if you are seeking a steadystate 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 

February 7, 2009, 06:37 
a) Okay
b) Okay, I understand

#13 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
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 

February 7, 2009, 16:07 
Vishal,
@I) I tried to sugg

#14 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9 
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 

February 7, 2009, 16:33 
Dear Henrik,
The equations

#15 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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 

February 8, 2009, 06:22 
Shouldn't this be written as:

#16 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
Shouldn't this be written as:
Dragos 

February 8, 2009, 11:26 
Your entry div(((Z1*C1)*fvc::

#17 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
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 

February 8, 2009, 11:51 
It seems there is an extra spa

#18 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
It seems there is an extra space in div(((Z1*C1)*grad(Ph i)))


February 8, 2009, 12:05 
Dear Dragos,
No I just chec

#19 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 9 
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/OpenFOAM1.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 

February 9, 2009, 02:30 
Hello guys,
Can anyone tell m

#20 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
Hello guys,
Can anyone tell me why the following expression gives me a compiler syntax error?
Dragos 

Thread Tools  
Display Modes  


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 