[zakifoam]

Hi all,

I am trying to create an openfoam module to simulate phase changes (i.e. melting & solidification) using enthalpy-porosity technique (following Brent et al, 1988), that allows the use of fixed mesh. Basically, the module originates from icoFoam, added with energy equation (in terms of temperature), buoyancy term in momentum equation, and extra source terms in both momentum & energy equation to account for the phase change.

My temperature equation looks like the following:

volScalarField dDHdt = fvc::ddt(DH);

{
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
==
dDHdt/Cp
);

TEqn.solve();
}

where DH is the latent heat content of the cell, evaluated based on cell temperature. Assuming isothermal phase change, DH is zero if cell temperature is lower than melting point, or equals to latent heat of fusion if cell temperature is higher than melting point.

After running the melting simulation, I found out that the energy source term dDHdt/Cp is always zero throughout the computation. This makes me think that this term is still incorrectly handled in T equation above. I suppose, instead of explicitly specifying dDHdt, I should've put the time derivative of DH on the l.h.s implicitly along with T (referring to relevant thread here). But in this way, I still couldn't get any effect from the source term.

Has anyone ever tried to work on this problem? Any share of ideas & comments is welcome.

Thanks

zaki

[hartinger]

Hi Zaki,

i assume you update DH before using fvc::ddt(DH), right? That happend to me once.

regards
markus

[zakifoam]

Hi Markus,

Thanks for your response. Yes, I updated DH before using fvc::ddt(DH). This is carried out along with update of A for source term A*U in momentum equation that will force velocity to zero in solid cells.

Any suggestion for remedy?

regards,

zaki

[hartinger]

not really. just the usual, trace back the calculation of DH and check what makes that zero or constant...