[Philipadebayo]

Dear All,

I am currently working on modelling of melting and solidification of phase change materials. I have been able to develop a solver using piecewise enthalpy-temperature relationship with the help of this thread. Total liquid volume fraction vs Time

I was able to validate the model with some cases in the literature.

Now, I want to implement the experimental enthalpy-temperature relationship based for binary mixture as defined in the attached picture.

I implemented the equation in openfoam in the T.Eqn from the mentioned thread as shown below

{

volScalarField h2 = (cpL - cpS)*(Ta - Tm)*log(mag(h1)) + L/h1;
fvScalarMatrix TEqn
(
fvm::ddt(cp, T)

+fvm::div(phi*fvc::interpolate(cp), T)
- fvm::laplacian(DT/rho, T)
==- fvc::ddt(h2)- fvc::div(phi, h2)

);

TEqn.relax();

TEqn.solve();

rhok = 1.0 - beta*(T - TRef);

volScalarField Tstar = (Ta - Tm)*alpha + Tm;

DH = DH + omegaDH*cp*(T - Tstar);

}

and in the source.H file as below

volScalarField h1 =(T - Ta)/(Tm - Ta);

forAll(mesh.cells(),celli)
{
if (T[celli] < Tm.value() && Tm.value()==T[celli])
{
h1[celli]= (T[celli]-Ta.value())/ (Tm.value()-Ta.value());

}
else
{
h1[celli]=1;
}

if (DH[celli] > L.value())
{
DH[celli] = L.value();
}
if (DH[celli] < 0)
{
DH[celli] = 0;
}

};

alpha=DH/L;

A = -C*sqr(scalar(1)-alpha)/(rho*(pow(alpha,scalar(3))+C1));


I was able to compile and run the solver but the melting front progression is faster than the piecewise function for the same simulation parameter. I suspect that I may not be implementing the enthalpy-temperature relationship in the right way. I would appreciate any help and insight on how to better implement this code.

Thank you.