[zakifoam]

Dear all,

I am interested in the VOF simulation of free surface flow driven by surface tension gradient due to surface temperature variation.

To do this, I have tried to add energy equation (temperature-based) in interFoam, in which the thermal conductivity is weighted by volume fraction, just like density and viscosity.

Additionally, at the interface, the smeared surface force has the normal and tangential components of the following form:



where



Here sigma_0 is surfac tension at reference temperature, sigma_T is dsigma/dT, i.e. the temperature gradient of surface tension, and T_0 is the reference temperature.

The first term of the normal force already exists in interFoam, so I need to implement the second term of the surface force and the tangential force.

I ended up with the following in UEqn.H :

surfaceScalarField muf = twoPhaseProperties.muf();

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
//- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);

if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma)
+ fvc::interpolate(dsigmadT*(T-T0)*interface.K())*fvc::snGrad(gamma)
+ fst
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}


and in pEqn.H :

phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma)
+ fvc::interpolate(dsigmadT*(T-T0)*interface.K())*fvc::snGrad(gamma)
+ fst
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();


Here, I have implemented the second term for the normal force (highlighted in blue), but not the whole tangential component of the surface force (in red).

This tangential part is a bit confusing, since in UEqn.H, everything in fvc::reconstruct brackets should be of the type surfaceScalarField. I tried to specify fst as follows:

volVectorField nhat = gradGamma/(mag(gradGamma) + deltaN);

volVectorField fst = dsigmadT*(fvc::grad(T) - nhat*(nhat & fvc::grad(T)))*mag(gradGamma);

Here, fst is a volVectorField, which can not be used in UEqn.H and pEqn.H above.

Could anyone help me with this problem?

Thanks a lot,


Zaki