laplacianfoam with solidification
 

Hi,
I'm trying to add solidification to laplacianfoam (I'm not interested in convective flux, but only in conduction with heat source).
My problem is that the temperature I have calculated in a simple test doesn't seem correct (a simple square cavity with uniform temperature inside, 700°C, and fixed temperature to the rigth wall, 500°C).
The temperature decreases too fast!

Modification to laplacianfoam:

Info<< "\nCalculating temperature distribution\n" << endl;

for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readSIMPLEControls.H"

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
solve
(
rho * cp * fvm::ddt(T)
+ Lh * rho * fvc::ddt(alpha1)
- fvm::laplacian(kT, T)
);
alpha1 = max(min((T-Ts)/(Tl-Ts),scalar(1)),scalar(0));
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Min(alpha2) = " << min(alpha2).value()
<< endl;
alpha2 = scalar(1) - alpha1;
rho = alpha1*rho1 + alpha2*rho2;
cp = alpha1*cp1 + alpha2*cp2;
kT = alpha1*kT1 + alpha2*kT2;
}

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return(0);

Modification to Createfield:

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);


dimensionedScalar rho1
(
transportProperties.lookup("rho1")
);

dimensionedScalar rho2
(
transportProperties.lookup("rho2")
);

dimensionedScalar cp1
(
transportProperties.lookup("cp1")
);

dimensionedScalar cp2
(
transportProperties.lookup("cp2")
);

dimensionedScalar kT1
(
transportProperties.lookup("kT1")
);

dimensionedScalar kT2
(
transportProperties.lookup("kT2")
);

dimensionedScalar Tl
(
transportProperties.lookup("Tl")
);

dimensionedScalar Ts
(
transportProperties.lookup("Ts")
);

dimensionedScalar Lh
(
transportProperties.lookup("Lh")
);

Info<< "Reading field T\n" << endl;

volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

Info<< "Reading field alpha1\n" << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh
);

Info<< "Calculating field alpha2\n" << endl;
volScalarField alpha2("alpha2", scalar(1) - alpha1);

Info<< "Reading field rho\n" << endl;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh
);

Info<< "Reading field cp\n" << endl;
volScalarField cp
(
IOobject
(
"cp",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh
);

Info<< "Reading field kT\n" << endl;
volScalarField kT
(
IOobject
(
"kT",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh
);

Can someone tell me where am I wrong?

Dario

Hi,
the problem is in the source term that is always 0.
How can I use the old value of alpha1, to compute:

Q = Lh*rho*(alpha1-alpha1OLD)?

Dario