implementing this equation
 

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 x-velocity and y-velocity 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!!

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

Dear Henrik,
Thank you so much for the advice,it's kind of you.
Sorry for the delay in replying you.

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.065e-3
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.065e-3*(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 == nCorr-1 && 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??

Maybe ....
 

Hi isabel, maybe you can try to add an artifical compression term into this equation as similar as gamma*(1-gamma)*U of the gammaEqn.H? I don't know. But, why you not try it?

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.

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 Diagonal-ILU or Multi-grid. What is the matter with them? You think.

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.

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...

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... :o