
[Sponsors] 
July 10, 2009, 05:27 
implementing this equation

#1 
New Member
andrianomena
Join Date: Mar 2009
Location: madagascar
Posts: 23
Rep Power: 9 
Hi there, It's been a whlile I haven't logged in, Please,Can anybody show me how to implement this equation of species in a chemical model? d/dx(rho*Ux*xi)+d/dy(rho*Uy*xi) = Sxi d/dx and d/dy stand for the derivative with respect to x and y Ux,Uy are the xvelocity and yvelocity components xi is the mass fraction of a species Sxi is the source term that generates the species at this stage the pressure,the temperature,U and rho have just been solved Here is what I tried but it diverged: fvm::div(phiv,xi) = Sxi; phiv is surfaceScalarField ("fvc::interpolate(rho*U) & mesh.Sf()") when I run it,it crashes.(error... prinstack...sigpFe) So my question is,how to implement this in the code? Cause I'm quite sure that problem lies in the way I implement it but I just don't know how.I'm not very experienced in OpenFoam Thanks!!


July 12, 2009, 14:10 

#2 
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 277
Rep Power: 10 
Dear sambatra,
1) Make sure the flux is conservative 2) Make the source implicit 3) Underrelax Read a book on Finite Volumes if this does not ring a bell. Henrik 

July 14, 2009, 05:26 

#3 
New Member
andrianomena
Join Date: Mar 2009
Location: madagascar
Posts: 23
Rep Power: 9 
Dear Henrik,
Thank you so much for the advice,it's kind of you. Sorry for the delay in replying you. 

August 14, 2009, 14:28 

#4 
Senior Member
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 9 
Sorry for bother you.
I am studying boiling (growing of a vapor bubble in a liquid medium). I have to add a mass source to the continuity equation: div(u) = cte*grad(T)*grad(gamma) gamma is the liquid volume fraction, and cte=1.065e3 So I have programmed the pEqn.H as this: volScalarField rUA = 1.0/UEqn.A(); surfaceScalarField rUAf = fvc::interpolate(rUA); U = rUA*UEqn.H(); surfaceScalarField phiU ( "phiU", (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi) ); phi = phiU + ( fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma)  ghf*fvc::snGrad(rho) )*rUAf*mesh.magSf(); adjustPhi(phi, U, pd); volVectorField gradT = fvc::grad(T); volVectorField gradgamma = fvc::grad(gamma); sourcemasa = 1.065e3*(gradT & gradgamma); for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pdEqn ( fvm::laplacian(rUAf, pd)  fvc::div(phi)  sourcemasa ); pdEqn.setReference(pdRefCell, pdRefValue); if (corr == nCorr1 && nonOrth == nNonOrthCorr) { pdEqn.solve(mesh.solver(pd.name() + "Final")); } else { pdEqn.solve(mesh.solver(pd.name())); } if (nonOrth == nNonOrthCorr) { phi = pdEqn.flux(); } } U += rUA*fvc::reconstruct((phi  phiU)/rUAf); U.correctBoundaryConditions(); The solver runs Ok, but I have convergence problems. Does anybody know why?? 

August 16, 2009, 08:36 
Maybe ....

#5 
Senior Member
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 10 
Hi isabel, maybe you can try to add an artifical compression term into this equation as similar as gamma*(1gamma)*U of the gammaEqn.H? I don't know. But, why you not try it?


August 20, 2009, 06:48 

#6 
Senior Member
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 10 
Hi isabel, I just remind it. If you add a source term to pdEqn, you think, whether or not this equation can still keep symmetrical?? If not, you should choose PBiCG method to solve it (in system / fvSolution) but not PCG; if yes, you can still use PCG. I can not make sure, because I don't know how to analyze the coefficient matrix structure to the momentum, continuity and gamma fraction volume equations in OpenFOAM, f. e. their eigenvalues, diagnal predominance et. at. Do you have some direct methods? Thanks.


August 20, 2009, 10:05 

#7 
Senior Member
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 10 
Hi isabel, in fact, I ever tested in my own code (not OpenFOAM) that, whether the pressure equation is symmetrical or not, it can been solved very well by PBiCG method. So, please don't worry about anything and try it.
But my preconditioning method is just ILU(0) but not DiagonalILU or Multigrid. What is the matter with them? You think. 

August 26, 2009, 06:42 

#8 
Senior Member
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 9 
Thank you very much, sandy, but when I try to change PCG to PBiCG, I obtain this error:
Unknown symmetric matrix solver PBiCG Valid symmetric matrix solvers are : 4 ( ICCG smoothSolver PCG GAMG ) It means that my matrix is still symmetric despite I have added a source. 

August 26, 2009, 10:50 

#9 
Senior Member
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 10 
Hi isabel, I wake up again. Source term = Sp * psi +Sc, so Sp actually is acutally added into diagonal elements in the ldumatrix. So, the pressure equation is still symmetric. If you add a ddt term as like as the sonicTurbFoam, the pressure equation will loss symmetric, I think. Sorry...
Last edited by sandy; August 27, 2009 at 08:53. 

August 27, 2009, 08:55 

#10 
Senior Member
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 10 
Hi isabel, I wake up again. Source term = Sp * psi +Sc, so Sp actually is acutally added into diagonal elements in the ldumatrix. So, the pressure equation is still symmetric. If you add a ddt term as like as the sonicTurbFoam, the pressure equation will loss symmetric, I think. Sorry...


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Calculation of the Governing Equations  Mihail  CFX  7  September 7, 2014 06:27 
continuity equation  Rafal  Main CFD Forum  4  November 29, 2006 10:27 
Implementing additional equation?  Pratap  FLUENT  1  December 4, 2004 21:32 
implementing continuity equation at walls  Subhra Datta  Main CFD Forum  3  November 27, 2003 00:57 
Diffusion Equation  izardy amiruddin  Main CFD Forum  2  July 4, 2002 08:14 