[mm66]

Dear Foamers,

I am trying to simulate the melting process of a solid body (phase change material) in 3D using icoReactingMultiPhaseInterFoam (first question: is this the best solver for melting? ). I included the buoyant melted movement by Boussinesq as the equation of state for the liquid. The problem is that when I run the simulation, it goes fine until towards the end of melting when suddenly the Courant number goes higher and higher until the alpha value becomes negative... I tried this in 2D and it never happened. But in 3D despite using a relatively fine mesh this happens... I also tried with adjustableRunTime which went down to 10^-6... I also checked the mesh quality:

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

Create mesh for time = 0

Time = 0

Mesh stats
    points:           206283
    faces:            2290188
    internal faces:   2234676
    cells:            1131216
    faces per cell:   4
    boundary patches: 3
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     0
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    1131216
    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
    Symmetry            19091    9950     ok (non-closed singly connected)
    Tube                18819    9679     ok (non-closed singly connected)
    Wall                17602    9007     ok (non-closed singly connected)

Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
    No cellZones found.

Checking geometry...
    Overall domain bounding box (0 0 0) (0.175 0.175 0.5)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.30135e-17 -3.18322e-16 2.85069e-16) OK.
    Max cell openness = 2.33158e-16 OK.
    Max aspect ratio = 4.71177 OK.
    Minimum face area = 6.46877e-07. Maximum face area = 4.8044e-05.  Face area magnitudes OK.
    Min volume = 2.68775e-10. Max volume = 8.59991e-08.  Total volume = 0.01183.  Cell volumes OK.
    Mesh non-orthogonality Max: 52.2945 average: 14.4077
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.623961 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

Here is the phaseProperties file:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1812                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      phaseProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
type    massTransferMultiphaseSystem;

phases  (solid liquid);

liquid
{
    type            pureMovingPhaseModel;
}

solid
{
    type            pureStaticSolidPhaseModel;
}

interfacePorous
(
    (solid and liquid)
    {
        type            VollerPrakash;
        solidPhase      alpha.solid;
        Cu              1e5;
    }
);

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

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

And here is the controlDict:

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

application     icoReactingMultiphaseInterFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         4000;                                                                                                                                                                                           deltaT          0.01;                                                                                   
writeControl    runTime;

writeInterval   5;

purgeWrite      0;

writeFormat     binary;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable no;

adjustTimeStep  no;

//maxDeltaT       0.01;

//maxCo           1.;
//maxAlphaCo      1.;
//maxAlphaDdt     2.3;
//maxDi           25;

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

Any ideas? Is there a way let's say in fvOptions to set upper and lower limits for the alpha value between [0,1]?

Thanks,
MJ

PS: Now I am trying to further refine the mesh, hoping it is possible to do so...

[wolfindark]

same problem i have....

[TeresaT]

Hi,

I think this is the right kind of solver. I am working with icoReactingMultiphaseInterFoam as well but on a different topic. A Laser melts my material. You said "towards the end" and "same problem" so this happens when both of your solids are melted completely and there is only liquid left?

Regards,
Teresa

[wolfindark]

Originally Posted by TeresaT

Hi,

I think this is the right kind of solver. I am working with icoReactingMultiphaseInterFoam as well but on a different topic. A Laser melts my material. You said "towards the end" and "same problem" so this happens when both of your solids are melted completely and there is only liquid left?

Regards,
Teresa

I also try to simulate laser melting. but at the same time liquid metal evaporation. When I added liquid to gas evaporation, simulation diverges.

[TeresaT]

Liquid to gas is my next goal and I just startet the first simulation with it included. I first wanted to see that melting works fine.

I will see if my case will diverge too.

[mm66]

Thanks for replying
Mine is pure melting of a solid. It runs well let's say up to 60% liquid (this value changes based on mesh quality, etc.) and then Courant number begins increasing until it messes everything (alpha becomes negative and divergence happens...)
Have tried everything that might affect this (changing settings in fvSolutions, very fine mesh, etc.) still getting the same error...
Any idea/input is appreciated...

Regards,
MJ

[TeresaT]

You probably have checked the quality of your mesh and if it crashes when liquid flows to a particular area?

Checking the timesteps prior to crashing for high velocities or pressure jumps could help with this.

In my case with adjustable timeSteps I have a deltaT of 2e-7, but my mesh is about a factor of 1000 smaller.

Does your case crash at some point and do you have the error message at the end of the log at hand?

[mm66]

Originally Posted by TeresaT

You probably have checked the quality of your mesh and if it crashes when liquid flows to a particular area?

Checking the timesteps prior to crashing for high velocities or pressure jumps could help with this.

In my case with adjustable timeSteps I have a deltaT of 2e-7, but my mesh is about a factor of 1000 smaller.

Does your case crash at some point and do you have the error message at the end of the log at hand?

Thank for you prompt response...
I always check the quality of the mesh before running the simulations. Indeed, they have a very good quality but refining it has not resulted in convergence yet
I also ran the simulations with/without adjustableTimeStep. With it, the timestep goes down to 10^-100
I tried several meshes, every time refining it further, same error appears every time
Attached you can find the results I got from one of the simulations.
BTW, I ran these simulations on an HPC cluster. They only error that I got is:

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 1840 RUNNING AT n122.localdomain
=   EXIT CODE: 8                                                                                        =   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
   Intel(R) MPI Library troubleshooting guide:
      https://software.intel.com/node/561764
===================================================================================

Any ideas what might be causing this? Thank you very much.

Regards,
MJ

[mm66]

Originally Posted by TeresaT

You probably have checked the quality of your mesh and if it crashes when liquid flows to a particular area?

Checking the timesteps prior to crashing for high velocities or pressure jumps could help with this.

In my case with adjustable timeSteps I have a deltaT of 2e-7, but my mesh is about a factor of 1000 smaller.

Does your case crash at some point and do you have the error message at the end of the log at hand?

BTW, may I ask what you mean by "but my mesh is about a factor of 1000 smaller"?

[TeresaT]

Hi MJ,

I have a gridwidth of 2.5 µm, I think a grid study will show that I am still to big. The over all domain is about 1x0.2x0.175 mm during my tests. That's was what I meant with mine beeing smaller.

The error code of MPI is not very helpful, at least not for me - doesn't spark any ideas.
What about your log file? The legends of your plots aren't always visible you might want to change that for better understanding.

Have you checked if there are high velocity gradients somewhere? Something happens between the 1500th and the 2000th Iteration and looking there might spark some ideas.

Regards,
Teresa
Edit:
One more thing: What Boundary Conditions do you use?

[mm66]

Originally Posted by TeresaT

Hi MJ,

I have a gridwidth of 2.5 µm, I think a grid study will show that I am still to big. The over all domain is about 1x0.2x0.175 mm during my tests. That's was what I meant with mine beeing smaller.

The error code of MPI is not very helpful, at least not for me - doesn't spark any ideas.
What about your log file? The legends of your plots aren't always visible you might want to change that for better understanding.

Have you checked if there are high velocity gradients somewhere? Something happens between the 1500th and the 2000th Iteration and looking there might spark some ideas.

Regards,
Teresa
Edit:
One more thing: What Boundary Conditions do you use?

Hi Teresa,

Thanks for the clarification. My simulation domain is a quarter of a cylinder with 0.5 m height and 0.35 m in diameter.

The error code is not helpful at all. That is the only error that I have at the end of the log file...

I am terribly sorry for the legends. Attached please find the new plots with updated legends and titles.

I tried to access the files (from the server) to verify if there are high velocity gradients. But the server seems to be down... I will check this as soon as I can.

Regarding boundary conditions, there is no inlet or outlet to the domain, it is a confined solid material undergoing phase change, so I am using the following simple boundaries at the walls:

U: fixedValue uniform (0 0 0)
T: one side is fixedValue, the rest are zeroGradient
alpha.liquid: zeroGradient
alpha.solid: zeroGradient
p_rgh: fixedFluxPressure

I really appreciate your help.

Regards,
MJ

[TeresaT]

Hi MJ,
thanks for the new plots!

Maybe I only suggest the following because I don't understand you case completely but without an inlet/outlet or atmospheric boundary condition you might want to look at your pressure as well.

I hope looking at velocity/pressure and temperature fields will give some hints.

Good luck
Teresa

[TeresaT]

Originally Posted by wolfindark

I also try to simulate laser melting. but at the same time liquid metal evaporation. When I added liquid to gas evaporation, simulation diverges.

Today I tried to add condensation and encountered something strange. The simulation crashed at the start with an error:
kineticGasEvaporation does not exist....

Found out that these:
(solid to liquid) with positive model constant c
(liquid to solid) with negative model constant c
(liquid to gas) with positive model constant c
work fine and I see the corresponding effects But as soon as i try to add
(gas to liquid) with negative model constant c
kineticGasEvaporation becomes unknown.
However,:
(liquid to gas) with negative model constant c
works fine again. Now I just have to find out if the condensation really works.

No, condensation does not work yet. I see lot's of vapour with temperatures between 364.3 and 587.3 Kelvin while condensation should take place at 2743.

[flo(w)]

I also had some problems with a suddenly increasing Co number. I recognized that the problem comes from the interface between two phases. It seems as if this is a common problem of some solvers, that instabilities may arise at the interface between two phases. I could solve the problem by adding some porosity to the interface through something like this in the phaseProperties file:

interfacePorous
(
(solid and gas)
{
type VollerPrakash;
solidPhase alpha.solid;
Cu 1e7;
}
);