Transient Phase Change Solver

[Roland]

Hi,
Since my last post the Crystal Growth model starts to work now. Here attached is a zipped animation (Crystal_Growth_Elmer.zip) showing the 1rst steps. The problem came from the fact that in the former version (explained in my former post), the solvers PS(solver 3: Transient Phase Change Solver) and CN (Solver 4: Compute Normal Solver) where duplicated twice in the case.sif file. After suppressing this duplication now it works with very small time steps (0.001s / with bigger time staps it does not converge) , but after a few hundreed time steps (about 0.3 seconds) the solver stops converging. Here attached is also the zipped project (CG_02_Elmer.zip).
Another thing which puzzles me : in this 1rst step model, there is a 0 heat flux on the horizontal free surface of the liquid (Body 1) and on the vertical right surface of the "solid crystal" (body 2). But when looking at the temperature Contours (here attached: temperature_Contours.png), we see that the contours are parallel to these boundary surfaces, which is not normal since, because of the 0 heat flux, these contours should be perpendicular to these boundary surfaces (!!).
It would be very kind if somebody could take a look at this and bring some ideas for fixing this problem.
Thanks in advance
Roland

[kevinarden]

I built a similar case with just the heat equation and the contours looked OK. I took your case and removed everything but the heat and the contours looked like yours. Then I realized in my case I had used steady state. So I too your modified case, heat only, and ran steady state, and then contours came out good. I then modified your case heat only back to transient and extended the time out. After about 0.5 seconds the contours started to look good. So it appears that early time there is very little heat flowing.

[kevinarden]

I then noticed that body 3 has no initial condition, so its starting temperature is 0 by default. Difficult to transfer heat into 0 in a short amount of time. I gave body 3 the same initial temperature as body 2, reset the time steps to your original values, 0.001 for 10 steps, and the contours immediately looked good on step 1. So I conclude the problem was body 3 is starting at 0 temperature.

[Roland]

Hello Kevin,
Many thanks for your reply! Yes you are completely right: I forgot to initiate the Body 3 air domain (my mistake!). Indeed this air domain was just forecast for introducing in a next step the magnetic field which must heat the liquid bath (Body 1). But for the moment, as the vertical right boundary is set by a Dirichlet BC at the hot temperature T_hot, this air domain Body 3 is not necessary and after suppressing it I came to the same conclusion as you, and now the temperature contours are Ok. I also disabled the Navier Stokes solver in the solid crystal domain Body 2, since it seems logical not to model this domain with Navier Stoke, as it is solid... Here attached is the corresponding zipped modified GUI project.
The solver works, but only with very small time steps (0.001s !) which takes a lot of time to reach 0.5s or 1s physical time, and strangely, after 0.5s the "solid/liquid" interface (boundary 4) does not seem to move...(for the moment the "v_cast" pulling velocity is set at zero, so that the interface should only move due to the Stefan phase change velocity. Here attached is a plot at 0.5s. As such a crystal growth process lasts physically several hours, I wonder how it could be possible to increase the time step in order to be able to model the process over for instance 500s or 1000s... Do you think that it is possible to change the solver settings so that the solver converges for instance with 1s (or 5s) time steps (instead of 0.001s)?
It would be very kind of you if you could take a look at this and give me your opinion about this strange behaviour. Say that I am not at all sure to set correctly the Transient Phase Change Solver and it would be welcome if you could explain how to manage this Phase Change Solver (I wonder also how it could be possible to set an imposed pulling velocity to control the triple point position / I also don't feel if the "v_cast" pulling velocity must be added to the "PhaseSurface" velocity at the interface or if it must rather be set as an imposed vertical upward velocity in the Navier Stokes BC at this interface boundary...).
Many thanks in advance for your help!
Roland

[kevinarden]

That is a lot to work through, but I will attempt.

Initial observation is that the the flow is turbulent. Maybe just barely, but once it is turbulent, Elmer is not very focused on turbulent flow. Turbulence is not just based on velocity, it is Reynolds's number, and you don't have to get very high for Elmer to start having convergence issues. There is a K-epsilon solver that can be added to allow for increased velocity. Small time steps help in turbulent flow, but heat flow takes a long time, competing objectives. If you can slow it down, increase the viscosity, the Reynolds's would be lower and convergence would improve.

I have a study on turbulent flow in Elmer here
https://github.com/mrkearden/turbulent

Some of the cases add the K-epsilon solver.

[kevinarden]

There is much interest and discussion on turbulent flow here;
viewtopic.php?p=23763&hilit=turbulent+flow#p23763

[Roland]

Hi Kevin,
Yes I agree with you concerning the turbulent flow. But in a first step I would like to have the model work and show a crystal growth (meaning have the solvers Compute_Normals for the liquid free surface and Phase_Change_Solver for the solid/liquid interface work) with a more or less low Reynolds number due to the fact that I already entered a bigger dynamic viscosity (0.7 Pa*s), say around 100 times bigger than the physical dynamic viscosity (we do that often in our Comsol models). But even with that "artificial bigger viscosity" (which can simulate turbulence on a simplier way than a real turbulent model like k-epsilon, ok?) the current model (attached to my former post) requires very small time steps and don't show a significative move of the interface, which I don't understand for the moment. Perhaps I made a mistake somewhere in the model....?
So, Kevin, it would be very nice if you could take a look at this model and tell me if there is something wrong somewhere.
Still many thanks for your help and we keep us informed.
Roland

[kevinarden]

This strategy ran out to 0.71 seconds fairly quickly. I did reduce the density and increased the viscosity. It also produces this warning
WARNING:: ListFind: Requested property: [Surface Smoothing Factor], not found
I also changed settings in stokes.
The mesh is starting to distort severely, and much farther the element formulations start to fail.

ycase.sif

(9.16 KiB) Downloaded 27 times

5thstep.png

(164.5 KiB) Not downloaded yet

last_step.png

(180.86 KiB) Not downloaded yet

[Roland]

Hi Kevin,
Thanks again for your efforts.
This thing seems not so easy... I think that it is better to model again the solid crystal with Navier Stokes and faking it with a very big viscosity(as I did it before). So I will enable again Navier Stokes in this body 2. In think that I will also add a small meniscus between the liquid and the solid, so that the triple point is free to go upward or downard depending of the Stefan velocity sign. What do you think of that? My hope is also to be able to increase the time step in order to run the solve over several hundreds of seconds. Have you some ideas about how doing that?
We keep us informed and again many thanks

[kevinarden]

I would concentrate on the stokes equation it is going to limit the time step size. I would turn every thing off and only model a flow using a desired assumed/anticipated velocity. Once the stokes is converging well in achievable time steps than you can add in the other features.