[wildfire230]

Hi all,

I am basically working with a modified scalarTransportFoam that considers the motion of two different scalar fields, one is some dissolved solute, and the other represents the concentration of a particle suspension. There is a phenomenon known as diffusiophoresis in which particles spontaneously move in the presence of a solute gradient.

Long story short, I need to modify the advection-diffusion equation with an additional term: gammap * grad(Cions) / Cions, where gammap is a constant and Cions is my solute ion concentration, which is updated every time step. That update just uses basically exactly the same equation as scalarTransportFoam, and that is working correctly.

Now, in my equation to update the particle concentration, I define a new updated velocity and phi:

volVectorField U_New = U + gammap / Cions * fvc::grad(Cions);
surfaceScalarField phiNew = linearInterpolate(U_New) & mesh.Sf();

and then my CpartEqn (particle concentration update equation) looks like this:

while (simple.correctNonOrthogonal())
{
volVectorField U_New = U + gammap/Cions*fvc::grad(Cions);
surfaceScalarField phiNew = linearInterpolate(U_New) & mesh.Sf();

fvScalarMatrix CpartEqn
(
fvm::ddt(Cpart)
+fvm::div(phiNew,Cpart)
-fvm::laplacian(Dpart,Cpart)
);

CpartEqn.relax();
fvOptions.constrain(CpartEqn);
CpartEqn.solve();
fvOptions.correct(Cpart);
}

However, I am getting some strange results, and it looks like something funny is happening in the phiNew step of the calculation. For example, if I just forget about updating U and use

volVectorField U_New = U;
surfaceScalarField phiNew = linearInterpolate(U_New) & mesh.Sf();

Then my updated Cpart actually gives different results than if I just use phiNew = phi;

Does that make sense??? Why would phiNew = linearInterpolate(U_New) & mesh.Sf(); give different results than phiNew = phi???

I can't understand.

Please help!

Thanks