[Unkinunki]

Hello, everyone.

I am trying to validate the melting and solidification model by reproducing the experimental procedure found in (Jones et al., 2006). This model was also used by (Muhammad el al., 2015) and it is attached here. To simulate this experiment, I am using OF8 and the solver chtMultiRegion, in which the solid regions are the acrylic and the polycarbonate, whereas the pcm material is the fluid region.

As a first approach, I used constant thermal properties in all materials and the solution easily converged. Next, I used variable thermal properties for the fluid region, except the density. Again, the solutions converged well and some results are shown in the second figure.

My problem has started when I assumed a polynomial function or the Boussinesq approximation for the density. The values of p_rgh, p and Co explode and the solution displays horrible results in the fluid region after some iterations, despite not crashing. I suspect that this behavior occurs when the temperature begins to change in the fluid region, since it takes some time to diverge. Trying to solve this problem, I tried the following procedures:

• Refine Mesh
• Reduce time step to 0.00001 (Variables only take longer to diverge)
• Use an adjustable time step, which presented values in the order of 10^-15 (Co adjusts but solution stills messy, same as depicted in the figures)
• set dpdt off when I used the sensibleEnthalpy
• Limit pressure values by using a pMin = 99500 and pMax = 1e5;
• Relax p_rgh, U and e|h to 0.1.

None of these procedures have worked. So, I think that I might be setting the B.C in the fluid zone or using the pRef incorrectly. However, at this point, I am out of solutions in mind and if someone could help me, I would appreciate. I am attaching the temperature, velocity and p_rgh fields after the solution diverges.

Best regards,

#Code:

 
PIMPLE: Iteration 100

Solving for fluid region pcm
DILUPBiCGStab:  Solving for Ux, Initial residual = 0.025274, Final residual = 9.00165e-05, No Iterations 2
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.0200058, Final residual = 8.38601e-05, No Iterations 2
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.0911808, Final residual = 0.000685502, No Iterations 1
DILUPBiCGStab:  Solving for e, Initial residual = 1, Final residual = 0.00887018, No Iterations 2
Min/max T:296.15 343.15
DICPCG:  Solving for p_rgh, Initial residual = 0.43215, Final residual = 0.00425036, No Iterations 76
DICPCG:  Solving for p_rgh, Initial residual = 0.00382949, Final residual = 7.87853e-05, No Iterations 77
pressureControl: p min -7.88183e+28
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 11326.3, global = 0.0218294, cumulative = 0.0218294
DICPCG:  Solving for p_rgh, Initial residual = 0.108544, Final residual = 0.00107419, No Iterations 64
DICPCG:  Solving for p_rgh, Initial residual = 0.00104525, Final residual = 9.27925e-05, No Iterations 65
pressureControl: p min -7.88183e+28
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 11884.1, global = 0.0218294, cumulative = 0.0436589

Solving for solid region acrylic
DICPCG:  Solving for e, Initial residual = 0.0190569, Final residual = 1.83692e-06, No Iterations 1
DICPCG:  Solving for e, Initial residual = 1.8329e-06, Final residual = 3.44873e-10, No Iterations 1
Min/max T:296.15 328.937

Solving for solid region polycarbonate
DICPCG:  Solving for e, Initial residual = 0.00971627, Final residual = 1.11858e-06, No Iterations 1
DICPCG:  Solving for e, Initial residual = 1.11916e-06, Final residual = 2.52876e-10, No Iterations 1
Min/max T:296.15 343.15
PIMPLE: Not converged within 100 iterations
ExecutionTime = 26.66 s  ClockTime = 26 s

Region: pcm Courant Number mean: 1.53665e+10 max: 5.8385e+10
Region: acrylic Diffusion Number mean: 0.0100368 max: 11.2252
Region: polycarbonate Diffusion Number mean: 0.0110556 max: 0.0464376
Time = 0.07

Fluid Zone

0/pcm/T

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0/heater";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [ 0 0 0 1 0 0 0 ];

internalField   uniform 296.15;

    boundaryField
    {
        top
        {
            type            zeroGradient;
        }

        front
        {
		type			wedge;
	}
		
	back
        {
		type			wedge;
	}
		
	pcm_to_acrylic
	{
		type		compressible::turbulentTemperatureCoupledBaffleMixed;
		
		value		uniform 296.15;
			
		kappaMethod solidThermo;
			
		Tnbr		T;
			
	}

	pcm_to_polycarbonate
	{
		type		compressible::turbulentTemperatureCoupledBaffleMixed;
		
		value		uniform 296.15;
			
		kappaMethod solidThermo;
			
		Tnbr		T;
			
	}

0/pcm/p

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0/heater";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 1e5;

boundaryField
{

   top
        {
		type		calculated;
		value		$internalField;
        }

        front
        {
		type			wedge;
	}
		
	back
        {
		type			wedge;
	}

	pcm_to_acrylic
	{
		type		calculated;
		value		$internalField;
	}

	pcm_to_polycarbonate
	{
		type		calculated;
		value		$internalField;
			
	}

0/pcm/p_rgh

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0/fluid";
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [ 1 -1 -2 0 0 0 0 ];

internalField   uniform 99500;

boundaryField
{

        top
        {
		type		fixedFluxPressure;
		value		$internalField;
        }

        front
        {
		type			wedge;
	}
		
	back
        {
		type			wedge;
	}

	pcm_to_acrylic
	{
		type		fixedFluxPressure;
		value		$internalField;
	}

	pcm_to_polycarbonate
	{
		type		fixedFluxPressure;
		value		$internalField;
	}
    

}

// ************************************************************************* //

0/pcm/U

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0/fluid";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [ 0 1 -1 0 0 0 0 ];

internalField   uniform (0 0 0);

boundaryField

    {
        top
        {
			type		noSlip;
        }

        front
        {
			type			wedge;
		}
		
	back
        {
		type			wedge;
	}

	pcm_to_acrylic
	{
		type		noSlip;
			
	}

	pcm_to_polycarbonate
	{
		type		noSlip;
			
	}
		
	
}

// ************************************************************************* //

constant/pcm/thermophysicalProperties

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant/steel";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleInternalEnergy;
}

dpdt off;

mixture
{
    // PCM

    specie
    {
        molWeight       50;
    }
    
    equationOfState
    {
	rhoCoeffs<8>         (2456.56 -5.22 0 0 0 0 0 0 );
    }
    
    thermodynamics
    {
    Hf              0;
	Sf              0;
        CpCoeffs<8>         (-3726.924 18.96 0 0 0 0 0 0 );
    }
    transport
    {
		muCoeffs<8> 		(0.11994 -0.0006529 90 0 0 0 0 0);
		kappaCoeffs<8> 		(3.726502 -1.108e-2 0 0 0 0 0 0);
    }
}

// ************************************************************************* //

checkMesh

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 8-9b73cf21a682
Exec   : checkMesh
Date   : Sep 17 2020
Time   : 20:33:23
Host   : "fabricio-A320M-S2H"
PID    : 25153
I/O    : uncollated
Case   : /home/fabricio/Desktop/OF_cases/Phase_Change/Validation
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           40401
    internal points:  0
    faces:            80100
    internal faces:   39700
    cells:            20000
    faces per cell:   5.99
    boundary patches: 6
    point zones:      0
    face zones:       0
    cell zones:       3

Overall number of cells of each type:
    hexahedra:     19800
    prisms:        200
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    left                0        0        ok (empty)                        
    top                 100      201      ok (non-closed singly connected)  
    bottom              100      201      ok (non-closed singly connected)  
    right               200      402      ok (non-closed singly connected)  
    back                20000    20301    ok (non-closed singly connected)  
    front               20000    20301    ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (0 0 -0.00151577) (0.0347169 0.0658 0.00151577)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Wedge back with angle 2.5 degrees
    Wedge front with angle 2.5 degrees
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (1.71133e-14 5.13242e-18 -2.47923e-16) OK.
    Max cell openness = 3.18472e-16 OK.
    Max aspect ratio = 1.89533 OK.
    Minimum face area = 5.2623e-09. Maximum face area = 1.0472e-06.  Face area magnitudes OK.
    Min volume = 1.7313e-12. Max volume = 3.44528e-10.  Total volume = 3.46259e-06.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.330796 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

controlDict

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     chtMultiRegionFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         4000;

deltaT          0.01;

writeControl    timeStep;

writeInterval   1;

purgeWrite      0;

writeFormat     binary;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  no;

maxCo 1;

maxDeltaT 0.01;

// ************************************************************************* //

system/pcm/fvSchemes

*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)     bounded Gauss linearUpwind grad(U);
    div(phi,e)      bounded Gauss linearUpwind default;
    div(phi,R)      bounded Gauss linearUpwind default;
    div(phi,K)      Gauss linear;
    div(phi,Ekp)    Gauss linear;
    div(R)          Gauss linear;
    div(phiv,p)     Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


// ************************************************************************* //

system/pcm/fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{

   rho
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-6;
        relTol          0.1;
    }

    rhoFinal
    {
        $rho;
        tolerance       1e-6;
        relTol          0.1;
    }

   p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-4;
        relTol          0.01;
        maxIter         5000;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "(U|h|e).*"
    {
        solver           PBiCGStab;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.01;
	maxIter		5000;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors 100;
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
    
  pRefCell        0;
  pRefValue       0;


outerCorrectorResidualControl
    {
        p_rgh
        {
            relTol          0;
            tolerance       0.01;
        }
    }

}

relaxationFactors
{
    fields
    {
	"rho.*"		1;
    "p_rgh.*"       0.7;
    }

    equations
    {
        "U.*"       0.7;
        "(h|e)"     0.7;
    }
}




// ************************************************************************* //

References:

Jones, Benjamin J., et al. "Experimental and numerical study of melting in a cylinder." International Journal of Heat and Mass Transfer 49.15-16 (2006): 2724-2738.

Muhammad, M. D., O. Badr, and H. Yeung. "Validation of a CFD melting and solidification model for phase change in vertical cylinders." Numerical Heat Transfer, Part A: Applications 68.5 (2015): 501-511.

[jherb]

What I have done once, when the polynomial density function resulted in crashs, was to fit or better mis-use the janaf and sutherland methods for cp and mu (for the temperature/pressure range of interest) and used Boussinesq for the equation of state:

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       sutherland;
    thermo          janaf;
    equationOfState Boussinesq;
    specie          specie;
    energy          sensibleInternalEnergy;
}

(mis-use, because Ts of the sutherland fit was negative)


Another way, I got it to work, was to create a new equation of state mode for the liquid including the compressibility.



So something like:

makeThermos
(
    rhoThermo,
    heRhoThermo,
    pureMixture,
    polynomialTransport,
    sensibleInternalEnergy,
    hPolynomialThermo,
    icoPolCompress,
    specie
);

template<class Specie, int PolySize>
class icoPolCompress;

The icoPolCompress is more or less the icoPolynomial class of OpenFOAM. The main changes are:

            //- Is the equation of state is incompressible i.e. rho != f(p)
            static const bool incompressible = true;

and

template<class Specie, int PolySize>
inline Foam::scalar Foam::icoPolCompress<Specie, PolySize>::rho
(
    scalar p,
    scalar T
) const
{
    scalar rho0 = rhoCoeffs_.value(T);

    // calc compressibility

    return rho0*(1 + (p - 1e5)*compr); // apply it to density
}

I am not sure, why this works, but it does. Before I had the same problem as you report.

[Bloerb]

Your settings seem well suited for the problem. Are you certain that is it not a problem of your polynomials getting out of bound? If for some reason the temperature drops or increases in such a fashion that one of your values becomes negative your simulation is bound to blow up. Might want to add an fvOption to limit the temperature between those ranges, for numerical reasons. From a quick glance your temperature range in which all of your polynomials are positive is quite small.
Why the slight difference in initial values between p and p_rgh?
maxCo=1 seems a bit high and you did not specify maxDi. Although this does not seem diffusion bound. Nevertheless maybe 0.5 for both, especially with Euler time stepping.
The rest is your divSchemes. Does it work with upwind schemes? You are using bounded with your schemes. This is usually only done with steadyState cases. This shouldn't be there, but might be ok due to the outerCorrectors....
Since your mesh seems perfectly orthogonal, and only slightly skewed you could even consider switching your laplacian/snGradSchemes to uncorrected or orthogonal. Leads to a slightly increased stability and should be fine on such a nice grid.

In fvSolution consider setting the relTol to zero for the Final parameters.

[peterhess]

What about the gravity?

There is no inlet or oulet bc in the velocity field.

That why I guess, you have free convection.

The prgh field shows no gravity effect.

The temperature shows no gravity effect also.

Post the gravity field please.

Regards

Peter

[Unkinunki]

Sorry for the late reply and thanks for helping me.

So, the first mistake that I noticed was that the Boussinesq approximation is already implemented in the source term related to the melting/solidification process.

solidificationMeltingSource.C

forAll(cells_, i)
    {
        const label celli = cells_[i];

        const scalar Vc = V[celli];
        const scalar alpha1c = alpha1_[celli];

        const scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
        const vector Sb = rhoRef_*g*beta_*deltaT_[i];

        Sp[celli] += Vc*S;
        Su[celli] += Vc*Sb;
    }

Thus, it is not necessary to apply again in the thermophysicalProperties. There, one can use a constant density.

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleInternalEnergy;
}

dpdt off;

mixture
{
    // PCM

    specie
    {
        molWeight       50;
    }
    
    equationOfState
    {

	rhoCoeffs<8> (769 0 0 0 0 0 0 0);

    }
    
    thermodynamics
    {
    	Hf              0;
	Sf              0;
        CpCoeffs<8>         (-3726.924 18.96 0 0 0 0 0 0 );
    }
    transport
    {
		muCoeffs<8> 		(0.11994 -0.0006529 9e-7 0 0 0 0 0);
		kappaCoeffs<8> 		(2.16893 -5.89361e-3 0 0 0 0 0 0);
    }
}

The remaining mistakes were related to some polynomials being out of bound, as stated by Bloerb. I am sharing some files and results just for info, thanks again!

fvOptions

options
{    
    solidificationMeltingSource1
    {
        type            solidificationMeltingSource;
        active          yes;

        selectionMode   all;
        
       Tsol            309.55;
        Tliq		309.65;
        L               248000;
       thermoMode      thermo;
        beta            8.161e-4;
        rhoRef          769;
        Cu              1.6e+06;
        q               1.0e-03;
    }

limitT
{
    type            limitTemperature;
    active          yes;

    selectionMode   all;
   min               296.15;
    max             343.15;
}

}


// ************************************************************************* //

fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         leastSquares;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,e)      Gauss linearUpwind default;
    div(phi,K)      Gauss linearUpwind default;
    div(phiv,p)     Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

[Rakesh_Kumar]

Hello Foamers,
I'm also working on a similar problem, but also involves air phase at top of liquid initially, and also densities different in liquid and solid of PCM, so that a shrinkage void forms after solidification.
I have tried to setup this using interFoam but not handling, also tried to couple solidification fvOption with interFoam but I ran into errors.



Anyone know how to setup this problem???
If anyone has any idea of, how to setup this problem, please explain.
Any hep regarding this problem is greatly appreciated.



Thanks in advance.




-Rakesh