[Unkinunki]

Hello, everyone.

I using OF v8 to simulate case with phase change, in which the thermal properties are constant during the solid and liquid phases. As I cannot describe the thermal properties by using a polynomial function, I had to change the chtMultiRegionFoam to read tabulated properties. The problem is related to a drastic property change in a small temperature interval (0.1 K), which leads to the following error:

Region: pcm Courant Number mean: 5.9305e-07 max: 1.08311e-05
Region: acrylic Diffusion Number mean: 0.0100013 max: 11.2252
Region: polycarbonate Diffusion Number mean: 0.0112191 max: 0.04648
Time = 35.27


Solving for fluid region pcm
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.000943529, Final residual = 2.39963e-09, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.00156482, Final residual = 1.52596e-08, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1.57125e-05, Final residual = 1.79032e-10, No Iterations 1
DILUPBiCG:  Solving for e, Initial residual = 0.000157825, Final residual = 6.02885e-10, No Iterations 2
         [15] iter          Test           e/h          Cv/p          Tnew
         [15] 0        309.65       47981.1          2400       298.897
         [15] 1       298.897       1437.98          1926       309.662
         [15] 2       309.662         48011          2400       298.897
         [15] 3       298.897       1437.98          1926       309.662
         [15] 4       309.662         48011          2400       298.897
         [15] 5       298.897       1437.98          1926       309.662
         [15] 6       309.662         48011          2400       298.897
         [15] 7       298.897       1437.98          1926       309.662
         [15] 8       309.662         48011          2400       298.897
         [15] 9       298.897       1437.98          1926       309.662
         [15] 10       309.662         48011          2400       298.897
         [15] 11       298.897       1437.98          1926       309.662
         [15] 12       309.662         48011          2400       298.897
         [15] 13       298.897       1437.98          1926       309.662
         [15] 14       309.662         48011          2400       298.897
         [15] 15       298.897       1437.98          1926       309.662
         [15] 16       309.662         48011          2400       298.897
         [15] 17       298.897       1437.98          1926       309.662
         [15] 18       309.662         48011          2400       298.897
         [15] 19       298.897       1437.98          1926       309.662
         [15] 20       309.662         48011          2400       298.897

Well, I tried to relax the energy equation and reduce the time-step to approximately 0.0002. However, none of these atempts worked.

I don't know if I need to implement a new thermophysical model or change the temperature conversion method. I would appreciate if anyone could help me with this.

Best regards.

thermophysicalProperties

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       tabulated;
    thermo          hTabulated;
    equationOfState icoTabulated;
    specie          specie;
    energy          sensibleInternalEnergy;
}

dpdt off;

mixture
{
    // PCM

    specie
    {
        molWeight       50;
    }
    
    equationOfState
    {

	rho ((296.15 769) (343.15 769));

    }
    
    thermodynamics
    {
    	Hf              0;
	Sf              0;
        Cp 	((296.15 1926) (309.55 1926) (309.65 2400) (343.15 2400)) ;
    }
    transport
    {
		mu	((296.15 0.0061) (300 0.0055) (310 0.0043) (320 0.00349) (330 0.00284)  (343.15 0.00224));
		kappa 	((296.15 0.423) (309.55 0.423) (309.65 0.146) (343.15 0.146));
    }
}

fvSchemes

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind grad (U) ;
    div(phi,h)      bounded Gauss upwind default;
    div(phi,R)      bounded Gauss upwind 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;
}

fvSolution

solvers
{
    rho
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0.001;
        maxIter     	5000;
    }

    rhoFinal
    {
        $rho;
        tolerance       1e-7;
        relTol          0;
    }
   p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0.001;
	   maxIter		5000;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

   "(U|h|e|k|epsilon|R)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-7;
        relTol          0.001;
        maxIter     5000;
    }

    "(U|h|e|k|epsilon|R)Final"
    {
        $U;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
    

pRefCell 0;
pRefValue 0;

outerCorrectorResidualControl
   
 {
        p_rgh
        {
            relTol          0;
            tolerance       0.1;
        }

        e
        {
            relTol          0;
            tolerance       0.1;
        }

    }

}

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

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

controlDict

application     chtMultiRegionFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         4550;

deltaT          0.01;

writeControl    timeStep;

writeInterval   5000;

purgeWrite      0;

writeFormat     binary;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  no;

maxCo 0.1;

maxDi 0.1;

maxDeltaT 0.01;

[mikulo]

Hello Fabricio,

I have seen that before, the reason for that is the Newton-Raphson method in the thermoI.H. Since the function is discontinuous, the solution won't converge using Newton-Raphson as by default from OpenFOAM. One solution somehow is to use the Bisection method. Try to change a custom file, say myThermoI.H. Please tell me then if it works.

Mike

[erinsam]

Hello Michael and Fabrício,

I have been trying to use tabulated thermophysical properties just like Fabrício in OFv2006's compressibleInterFoam, but I am getting an issue that "Maximum number of iterations exceeded". Does this have to do with the convergence of the Newton-Raphson method used in the thermo package to obtain temperature from he (enthalpy or internal energy) that is solved for in the energy equation? Is that what you are suggesting Michael?

Thanks for reading.

Erin

[mikulo]

Originally Posted by erinsam

Hello Michael and Fabrício,

I have been trying to use tabulated thermophysical properties just like Fabrício in OFv2006's compressibleInterFoam, but I am getting an issue that "Maximum number of iterations exceeded". Does this have to do with the convergence of the Newton-Raphson method used in the thermo package to obtain temperature from he (enthalpy or internal energy) that is solved for in the energy equation? Is that what you are suggesting Michael?

Thanks for reading.

Erin

Hi Erin,

Yes, that could be one of the reasons. However, try to refine your mesh too since refining the mesh will minimize the effect of discontinuities. Since you are using tabulated properties it means that the slope of the enthalpy curve with respect to temperature (Cp) is changing from time to time so Newton-Raphson could be detrimental.

Mike