[danylo]

Hello Foamers!

Hope everyone is doing well.

I am trying to model a phase change of water to ice, which includes an expansion of my system due to the density difference of the phases. The cooling comes from the wall on the right side and the idea is, that expansion is inversely proportional to ice present in the top cell of a column, so frozen regions dont expand any further and my system is bulge shaped in the end. The expansion is realized by applying a calculated Point-Displacement-value to the top boundary. I attached screenshots, where you can see the beginning of this happening.

Now my problem is, that once the first frozen region slows its expansion down, the pressure field shows a weird behaviour in the corner. The cell in the corner has a higher pressure than the system (1e5), while the one right beneath it has a lower one. This trend accelerates as the simulation goes on and leads to a crash. And the pressure field is not the only one showing weird behavior due to the mesh motion (e.g. ice fraction tends to go over 100 % in this case). I guess the reason for this is, that the CellMotion, which OpenFOAM calculates, tends to zero in the top right corner, but the cells beneath it still have a relevant CellMotion value due to the moving neighbour cells (Screenshot of this is attached), but I dont know how I could counteract this.

My question is, what are the internal mechanisms of the OpenFOAM mesh-motion, which could lead to this an what can I do to prevent it from happening? Any help would be appreciated!

Some extra Info on my case:
BCs:
U
Top: pressureInletOutletVelocity (0 0 0)
Rest: noSlip

p
Top: totalPressure 1e5
Rest: zeroGradient

motionSolver:
velocityComponentLaplacian

diffusivity:
quadratic inverse distance to bottom wall

Down below are screenshots of the CellMotion in y direction, Ice Fraction and temperature.

Thank you in advance!

[danylo]

Addition:

I am using a customized version of rhoPimpleDyMFoam on OpenFOAM 5.

[danylo]

Update on my question

I improved the code insofar as the equations work, when the Mesh Motion is distributed evenly across the Mesh. To come by this problem of having a CellMotion in the middle of the cell column and none at the top, I now try to work with directional diffusivity, where I disable the diffusivity in x and z direction. So the line in dynamicMeshDict looks like this:

diffusivity: directional (0 1 0);
(I attached the rest of dynamicMeshDict as a screenshot)

In the beginning this works out nicely and I have next to no CellMotion in the Frozen Regions. But after a while the the internal cells start to have some strange Motion, while the top layer has none at all. So it seems like there still is a diffusivity in the x direction even though it should not be there. Some cells even move with a negative velocity which should not be possible at all (Either the whole mesh expands or the whole mesh contracts) and there are opposing movements in the same column. This skews my mesh and leads to a solver crash once again

I attached screenshots of how that looks.

Does anybody have an idea what could cause this?