add source term in energy equation
 

Hi, Foamers,
I added a source term in temperature equation, but the strange thing is temperature becomes unbounded hence unphysical. I have no any idea. The following are relevant codes.
The basic idea is : I added a source term in TEqn, and the source term is also dependent on T. When I deal with ddt(alphaTemp), I already stored .oldTime() each time.
I cannot see where the problem would be.

-----------------
fvScalarMatrix TEqn
(
fvm::ddt(cp, T)
+ fvm::div(phi*fvc::interpolate(cp), T)
- fvm::laplacian(lambda/rho, T)
==
ST
);
-------------------
ST = (
- hs*fvc::ddt(alphaTemp)
- hs*fvc::div(phi,alphaTemp)
);
-------------------
forAll(mesh.cells(), celli)
{
if (T[celli] < Ts.value())
{
alphaTemp[celli] = 0;
}
else if (T[celli] > Tl.value())
{
alphaTemp[celli] = 1;
}
else
{
alphaTemp[celli] = (T[celli]-Ts.value())/(Tl.value()-Ts.value());
}
};
----------------------
Any response is welcome and appreciated. :)

Any idea? It takes me quite a long time to figure out. :(
Many thanks to any comments. :)

I haven't looked too deeply at your code, and I don't know if there's some fundamental problem... like, say the source term has the wrong sign. I'll assume you know what you are doing, and everything is okay.

What I was going to mention is that numerical models sometimes behave very badly around switching functions, especially if you are taking derivatives of them. So your if/else structure may be causing the problems. You could try making it a smooth function without any conditionals... it will also make your code faster. In fact, you could probably make it a field operation instead of a cell-by-cell loop. Google logistic functions. Not sure if this will help, just an idea.

Hello David,
Thanks for your comments. I tried: alphaTemp = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5) which is a field operation, but still didn't figure out my problem. Also, the sign of source term has no problem.

I believe this is stability problem about solving TEqn with source term.
----- part of fvSolution -------
T
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0.1;
}

TFinal
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-9;
relTol 0;
}

PIMPLE
{
momentumPredictor no;
nOuterCorrectors 3;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
----------------------------------------------
My problem is a 2D cavity melting, initial temperature of the whole domain is 301.45, BC is left side 311.15 fixed, right side 301.45 fixed and up and bottom insulated.
---------------- part of log (temperature) --------------
Time = 0.05
uniform 301.45;
PIMPLE: iteration 2
(306.456 303.191 302.055 301.66 301.523 301.475 301.459 301.453 301.451 301.45)
PIMPLE: iteration 3
(218.682 189.548 262.537 287.918 296.744 299.814 300.881 301.253 301.384 301.434);

Time = 0.109375
(306.456 303.191 302.055 301.66 301.523 301.475 301.459 301.453 301.451 301.45);
PIMPLE: iteration 2
(308.363 304.81 302.974 302.112 301.73 301.565 301.497 301.469 301.457 301.452);
PIMPLE: iteration 3
(301.785 278.104 243.453 279.972 293.494 298.502 300.358 301.048 301.307 301.415)
--------------------------------------------
I also tried to set pimple iterations even as 25, still didn't work...
Is this stability problem, if yes, how can I handle it? Many thanks in advance.