[fabian_roesler]

Dear Foamers

For a project I am trying to get a solid/liquid phase change simulation running in OF 1.5-dev. A paraffin with non isothermal phase change is molten in a 2D rectangular test case scenario. I started from interFoam solver by changing especially the gamma equation to my needs:

Gamma equation:

dimensionedScalar pi = mathematicalConstant:i;
dimensionedScalar alpha = (Tl-Ts)/(scalar(2)*tan(scalar(0.49)*pi));
gamma = atan(scalar(1)/alpha*(T-(Tl+Ts)/scalar(2)))/pi+scalar(0.5);

Changig gamma from 0 (solid) to 1 (liquid) depending on the temperatures
Tl (upper phase change temperature) Ts (lower phase change temperature)

Momentum equation:

fvVectorMatrix UEqn

(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
);

UEqn.relax();

if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}

No changes at all, however no surface tension left

Energy equation:

volScalarField rhoCp("rhoCp", rho*cp);
surfaceScalarField rhoPhiCp("rhoPhiCp", rhoPhi*fvc::interpolate(cp));

fvScalarMatrix hEqn
(
fvm::ddt(rhoCp, T)
+ fvm::div(rhoPhiCp, T)
- fvm::laplacian(lambda, T)
);

solve
(
hEqn ==
- hs*fvc::ddt(rho, gamma) //Source term from Stefan condition
);

The source term derived from the Stefan condition takes care of the temperature gradient at the interface induced by the phase change.

Solver:

Like in interFoam, the solver works as follows:

1. Solving gamma equation
2. Solving coupled pressure and velocity equations
3. NEW solving energy equation NEW

The solver is very unstable what, in my opinion is due to the couppling of gamma and energy. Please correct me when I am wrong. My questions are:

Is it really the coupling of gamma and energy, making the solver instable?
Is there a way to solve the coupled equations (relaxation)?

Thanks allot for your replies.

Best regards

Fabian