Wrong sign for solid-to-liquid mass transfer in Lee model

Summary

The Lee massTransferModel does not correctly handle the phase change from liquid to solid.

In OpenFOAM v2012, the behavior is correct with the bug fix (commit 7830cc50). However, commit 7399dbfe introduced another issue, which has led to problems in OpenFOAM v2206 and later.

An example case and a proposed fix are attached.

Steps to reproduce

Run the icoReactingMultiPhaseInterFoam solver with a case that includes solidification. An example case is attached.

solidMelting2D_heating_cooling.zip

Example case

Base case: icoReactingMultiPhaseInterFoam/solidMelting2D

Changes:

• Added liquid-to-solid massTransferModel settings in constant/phaseProperties.
• Modified T boundary condition at the left wall to trigger solidification after 50 s.
  • 0–50 s: heating
  • after 50 s: cooling

Without the fix: the run takes about 1100 s to complete, and no solidification is observed.
With the proposed fix: the run takes about 400 s, and solidification starts after 50 s from the left wall.

What is the current bug behavior?

Liquid-to-solid transition is not reproduced. After 50 s, solid should form from the left wall.

Attached File: Left, ofirinal code; Right, modified code.

solidMeltingCooling_vector_T_solid.avi

What is the expected correct behavior?

modifiedCode_vector_T_solid.avi

Environment information

• OpenFOAM version: v2206 to v2506
• Operating system: Ubuntu
• Compiler: binary package

Possible fixes

Lee.C.patch

Remove the negative sign at Line 120 of Lee.C (in KSp function, for negative C value).

Present:

-coeff*pos(Tactivate_ - refValue)

Fix:

coeff*pos(Tactivate_ - refValue)

Related discussion

https://www.cfd-online.com/Forums/openfoam-bugs/231579-bug-lee-c-icoreactingmultiphaseinterfoam.html

The recommendation in #1934 (closed) (closed) (by @m-tamura) is to include a negative sign in the KSu function when C < 0. This is correct.

For C<0, KSp does not need a negative sign and KSu needs a negative sign.

Summary of current code

MassTransferPhaseSystem::heatTransfer

The contribution from phase change is calculated with
eqn -= fvm::Sp(Sp*L.ref(), T) + Su*L.ref();.
In simplified form:

(Sp T + Su) L = C rho alpha (T + Ta) / Ta * L

Case C > 0

(Sp T + Su)
= C alpha rho T / Tactivate - C alpha rho
= C alpha rho T / Tactivate - C alpha rho Tactivate/Tactivate
= C alpha rho (T - Tactivate) / Tactivate

Case C < 0 (as in current code)

(Sp T + Su)
= - C alpha rho T / Tactivate - C alpha rho
= - C alpha rho T / Tactivate - C alpha rho Tactivate/Tactivate
= - C alpha rho (T + Tactivate) / Tactivate

T + Tactivate should be T - Tactivate.

Correction

For C < 0, the negative sign should be removed.

Case C > 0 (unchanged)

(Sp T + Su)
= C alpha rho T / Tactivate - C alpha rho
= C alpha rho T / Tactivate - C alpha rho Tactivate/Tactivate
= C alpha rho (T - Tactivate) / Tactivate

Case C < 0 (after fix)

(Sp T + Su)
= C alpha rho T / Tactivate - C alpha rho
= C alpha rho T / Tactivate - C alpha rho Tactivate/Tactivate
= C alpha rho (T - Tactivate) / Tactivate

Here, C is negative. Interpreting this as transferring the sign into the temperature term, we can write it using the absolute value of C (Cabs):

(Sp T + Su) = Cabs alpha rho (Tactivate - T) / Tactivate

This correctly evaluates the term when the field temperature is below the phase-change temperature.

Version comparison (standard code)

OF	KSp C>0	KSu C>0	KSp C<0	KSu C<0
2006 NG	coeff*pos(refValue - Tactivate_)	-coeff*pos(refValue - Tactivate_)	coeff*pos(Tactivate_ - refValue)	coeff*pos(Tactivate_ - refValue)
2012 OK	coeff*pos(refValue - Tactivate_)	-coeff*pos(refValue - Tactivate_)	coeff*pos(Tactivate_ - refValue)	-coeff*pos(Tactivate_ - refValue)
2206–2506 NG	coeff*pos(refValue - Tactivate_)	-coeff*pos(refValue - Tactivate_)	-coeff*pos(Tactivate_ - refValue)	-coeff*pos(Tactivate_ - refValue)
2506 fix	coeff*pos(refValue - Tactivate_)	-coeff*pos(refValue - Tactivate_)	coeff*pos(Tactivate_ - refValue)	-coeff*pos(Tactivate_ - refValue)