[MamoruTamura]

Hello everyone,

I am trying a simulation of melting and solidification of matter, using icoReactingMultiphaseInterFoam in OpenFOAM v2006.
However, I found a phenomenon which seems to be a bug.
Since I'm a beginner of OpenFOAM, please let me consult about this.
To come to the point, I guessed there is a quite simple bug of a minus sign in Lee.C.
My simulation files are attached.
solidMelting2D_heating_cooling.zip

To try a simulation, I refered a tutorial of solidMelting2D, and changed its parameters, as follows.

In system/controlDict, endTime and writeInterval were changed,

endTime         200;
writeInterval   2;

In 0/T, left and right boundaryField were changed,

internalField   uniform 301.3;

boundaryField
{
    left
    {
    	type externalWallHeatFluxTemperature;
        mode power;
        Q table
        (
            (0      1000)
            (50     1000)
            (50.01  -1000)
            (100    -1000)
        );
        kappaMethod fluidThermo;
        value $internalField;
    }

    right
    {
        type           fixedValue;
        value          $internalField;
    }

In constant/phaseProperties, massTransferModel for a phase change from liquid to solid was add,

massTransferModel
(
    (solid to liquid)
    {
        type            Lee;
        C               40;
        Tactivate       302.78;
    }
    (liquid to solid)
    {
        type            Lee;
        C               -40;
        Tactivate       302.78;
    }
);

I changed the left boundary condition of T, to heat their with 1000W for the first 50 seconds, and after that, cool their with -1000W.
Thus, we can expect the melt and solidification of matter.

This is a result of simuation, alpha.solid is shown.
by default Lee.gif

First, as expected, the matter showed a melt.
However, in the middle of simulation (after t = 36s), the simulation showed a strange behavior.

I think there is a bug in the KSu function of Lee.C, because of the following reason.
For a different purpose, I investigated the latent heat calculation in icoReactingMultiphaseInterFoam.
It is calculated in the heatTransfer function of MassTransferPhaseSystem.C, whose line 294 is,

                eqn -= fvm::Sp(Sp*L.ref(), T) + Su*L.ref();

these Sp and Su would be from Lee class.

In Lee.C, KSp function calculates,

         volScalarField from
         (
             min(max(this->pair().from(), scalar(0)), scalar(1))
         );
  
         const volScalarField coeff
         (
             C_*from*this->pair().from().rho()*pos(from - alphaMin_)
             /Tactivate_
         );
  
         if (sign(C_.value()) > 0)
         {
             return
             (
                 coeff*pos(refValue - Tactivate_)
             );
         }
         else
         {
             return
             (
                 coeff*pos(Tactivate_ - refValue)
             );
         }

and KSu function calculates,

         volScalarField from
         (
             min(max(this->pair().from(), scalar(0)), scalar(1))
         );
  
         const volScalarField coeff
         (
             C_*from*this->pair().from().rho()*pos(from - alphaMin_)
         );
  
         if (sign(C_.value()) > 0)
         {
             return
             (
                 -coeff*pos(refValue - Tactivate_)
             );
         }
         else
         {
             return
             (
                 coeff*pos(Tactivate_ - refValue)
             );
         }

For (phase change from solid to liquid), these functions are,

for

Thus, the above mentioned MassTransferPhaseSystem.C line 294 would be,

eqn -=

which describes the energy source due to the mass change of phase change, and the mass change corresponds to the Lee model explained in,
https://www.openfoam.com/releases/op...nd-physics.php

On the other hand, for (phase change from liquid to solid), from above Lee.C code, Sp and Su are calculated as,

for

There is not a minus sign in Su calculation, thus, the sign of Sp and Su are identical; MassTransferPhaseSystem.C line 294 would be,

eqn -=

which is different with the Lee model shown in above link.
I guessed a minus sign is required in Su also for , like as,

and the case of KSu function should be,

         else
         {
             return
             (
                 -coeff*pos(Tactivate_ - refValue)
             );
         }

After this modification, I compiled icoReactingMultiphaseInterFoam, and retried the simulation.
This is a result.
by modified Lee.gif
As expected initially, the matter showed the melt and solidification.
I think the simulation was successful.

Let me hear what you think.

Best regards,
Mamoru

[MamoruTamura]

For this topic, I reported on

https://develop.openfoam.com/Develop.../-/issues/1934

It was solved.

Thanks!