romant |
April 13, 2010 05:06 |
temperature equation for one phase in multiphase solver
Hej,
I have been working with the interPhaseChangeFoam solver for a while now, and have started writing a phase change model based on temperature. Right now the temperature is still given via the dictionary input.
I wanted to add the energy equation for the fluid phase to this solver. The governing equation is
http://img541.imageshack.us/img541/3011/energyeq.png
which I tried implementing like this
Code:
{
// const volScalarField rhoAlpha = cp * rho1 * alpha1; // for time derivative
volScalarField rhoAlpha
(
"rhoAlpha",
cp * rho1 * alpha1
);
surfaceScalarField phiAlpha
(
"phiAlpha",
phi * rho1 * cp * fvc::interpolate(alpha1) //check with rhoPhi from the alphaEqn.H file
// rhoPhi * cp //* fvc::interpolate(alpha1)
);
volScalarField lambdaFAlpha
(
"lambdaFAlpha",
lambdaF * alpha1
);
Pair<tmp<volScalarField> > mDotAlphal = twoPhaseProperties->mDotAlphal();
const volScalarField mDotAlphac = mDotAlphal[0]();
fvScalarMatrix TEqn
(
fvm::ddt(rhoAlpha,T)
+ fvm::div(phiAlpha, T)
- fvm::laplacian(lambdaFAlpha, T)
== mDotAlphac * ifg
);
TEqn.relax();
TEqn.solve();
Info << "Max(T): " << max(T).value() << " K "
<< "Min(T): " << min(T).value() << " K" << endl;
}
Unfortunately, I have had no success in actually using this within the solver, the results are completely bogus. One thing that is obvious is that there might be a problem with the alpha1 attachment to the equation.
My first thought was to maybe be able to extract the mesh that contains water alpha1>0, and impose boundary conditions onto the interface cells. Is there a way to do it, or does somebody else see what could be wrong with the implementation of this equation?
|