Phase transition phenomenon ==> with Only TEqn.H
 

Hello Foamers,


# Post 19 --> Mr. Fabian Roesler had constructed the solver for melting problem in OpenFOAM 2.1.1 based on Stefan Task .
I have modified this solver in simplest way using SIMPLE algorithm based on only Temperature (TEqn.H) as shown in the attachment. Here, I have neglected pressure (pEqn.H) and velocity (UEqn.H) terms to fit the conditions based on my solver as explained in next comment section.
In this case, by default 'alpha' is a function of temperature, 'lamda' and 'cp' is the function of 'alpha'

This solver explains the movement of interface from right to left for non-isothermal phase change where melting occurs. [Stefan Problem - moving boundary is met]

Phase transition phenomenon ==> with Only TEqn.H
 

Hello Foamers,

Following my above comment, I would like to explain my problem conditions to fit into the phase change solver. After solving the simple energy equation, I want to include the phase transition phenomenon - Boiling, Condensation and Convection [ATTACHMENT 1].

ATTACHMENT 2 explains the schematic 1D model with phenomenon. Here the main equation is only ENERGY equ. The problem is basically the coal pyrolysis. Initially, the coal is considered as wet and 2 phases of moisture and steam need to be implemented. After evaporation of moisture, once if wet coal is converted into dry ~ pyrolysis will take place, which is a different topic of discussion.

Position of interface (boiling plane) is explicitly defined based on mass and heat balance ~ as not defined in previous Fabian Roesler problem.
Heat balance of boiling and condensation expressed in unit of W/m2 ==> where it couldnot be added as source term in energy equation which carries the unit W/m3 ~ And I am confused here
(1) Boiling: When the heating wall (left side) reaches the temperature of 100deg, the moisture evaporates [the condition (alpha) w = 0] instantly. [ATTACHMENT 3]
(2) Condensation: The steam condenses (right side of chamber) until the temperature reaches 100deg on the right side. [ATTACHMENT 4]
(3) Convection: When T = 100deg, the steam flow after moisture evaporation. [ATTACHMENT 5]
How to fix such condition at interface quoted below ==>
REF - MANUSCRIPT which I am focussed with based on GOVERNING EQUATIONS and SOLUTIONS ~ https://sci-hub.tw/https://doi.org/10.1016/0016-2361(83)90225-9
Although the equations to solve looks simple, I couldn't able to figure out the exact path to find my solution (considering my queries mainly marked in "RED mark")
Kindly someone share their ideas please. Thank you !!!

Hello Foamers,
To make it simple, I will explain the main core of the problem about phase change.
alpha as a function of T in melting problem ~ as mentioned in above explanation
But in my case, alpha is not function of T, instead the position of the interface (boiling plane) is fixed with the condition as follows:
//************************************************** *****************//
Heat is transferred by wall conduction only. When the temperature reaches T=100deg (heating wall), the surface reaches above the saturation temperature, and so liquid (moisture) evaporates leading to the motion of vapor-liquid flat interface. If the vapor-liquid interface moves to the complete right, then the chamber is saturated with only evaporated vapor. Interface stays flat for 1D PROBLEM.
//************************************************** *****************//
Following the above condition, mass and heat balance are calculated as,
REF - Clear details here in this reference ~ https://sci-hub.tw/https://doi.org/10.1016/0016-2361(83)90225-9

Kindly someone help me by sharing ideas..

Thank you

Is there any CHT melting solver that you know?


Thanks!

Hello Everyone,

First of all, I thank everyone who has contributed to this wonderful thread.

I downloaded the solver posted by Typhian in post #263. I am able to compile it successfully in OF V6 and perform the test case too.

But, when I changed the mesh resolution to 132X100 in the blockMesh dictionary. The case gives an abnormal result.

After changing the entries int he blockMeshdict, I ran blockMesh and meltFoam.

I believe it has something to do with the cellToRegion file.

Someone please explain me the use of cellToRegion file. I mean what it does exactly.

Sorry if my question is too silly.

Thank You.

Hello Pavithra,
I downloaded the solver and compiled in version 5.0 (post #263). U can find the same case in the attachment. I don't find any trouble while compiling and solving this case.
No such problems even if I changed the mesh size.
You have compiled in OF V6 - thats where the problem is, I hope so.
cellToRegion file is basically created by splitMeshRegions: This file defines all the patches of a particular region. cellToRegion( ) const --> Return region split.
cellToRegion file is like any other field file but not associated to a specific field. The BC type is either zeroGradient or calculated.

calculated (if defined) --> in case of coupling patches to other regions
zeroGradient (if defined) -->boundary patches (outside of domain)

Hope it helps.
Thank you ^^

Respected Sir,

Thank you so much for your kind reply.

I have attached the liquid fraction and temperature distribution at 2 seconds.

I find the liquid fraction and temperature distribution at 2 secs as abnormal. Kindly please give your comments on this.

So, in normal case, cellToRegion file does not have any use. The default case should produce the results even without the cellToRegion file. Sir, please correct me, if I am wrong.

Thank You.

In 1D case of mesh size (100 1 1) - the results looks reasonable.
However, for 2D case, the behavior is weird. Have you tried increasing the width size or decreasing the default time size less than 0.005 ?

Thank You Sir. I also believe that a smaller time step can reduce this error. I am running the case with a smaller time step.

But, Sir please help me understand that how increasing the width can resolve this issue.

Thank You

Hello Pavithra,
- Yes, no need of cellToRegion here.
I thought it could be corner effects - butt thats not the case here actually.
I am not proficient in multiphase problems. My field of work is different. However, I am sharing ideas to the best of my knowledge.
In the attachment given below, both images are captured at 170s. When the temperature range is between 300-310K, the contour plot looks apparent. But when the range is rescaled btw 300-340K, there found strange behavior - bottom left corner temperature rises. Time reduction doesn't helps. I am not sure where could be the problem.

~Have you gone through this post completely from the beginning. I am sure it will give you the spark.
~Have you tried the extended version coding by Fabian Roesler [post #81]
~Posts #130 #151 #154 draw contour results - might be helpful for you to validate
~Fabian Roesler - Ph.D thesis on post #180
~Ole Richter posts from #201-205 - about incompressible two phase mixture might help too..
Keep updating your progress. You will find a way one day!!

Thank you

Sir,



Thank you so much for your detailed reply. I will go through all the posts carefully and will update my progress here.



Thank you so much Sir.

Hello Everyone,

I am trying to use the solver by Fabian posted in this thread. Could someone please clarify me the following points?


1) What is Tdim in constant/transportProperties?


2) From his paper, I can see that Tdim = (T_liquidus) - (T_solidus). If this is true, for isothermal melting, Tdim = 0. If I set Tdim = 0, the solver crashes due to division by 0.


3) Moreover the value of Tdim affects the results, drastically. Smaller values of Tdim leads to accelerated melting and viceversa.


Please help me in choosing a right value of Tdim for isothermal melting. Fabian's solver and link to his paper are attached.

https://link.springer.com/article/10...231-011-0866-9

Thank You.

Hello Everyone,

Finally, I have managed to modify buoyantBoussinesqPimpleFoam to simulate melting problems. I have attached the solver and Gallium melting testcase. The solver is tested in OF v6.

Any feedback or suggestions are welcome.

Also, please give me a direction or any guidance to parallelize this solver.

The solver works well when running on a single core. But, blows up while running in parallel.

Thank You.

- Pavithra.

Hi Everyone,

I am attaching the results of some validation cases, herewith.

I have attached the solver, cases and results. Any suggestions for improvements are welcome.

Thank You.

- Pavithra

Kummi,


For a finite volume solver shouldn't the energy balance always be written in W/m^3? In this case there could be some confusion due to the length attributed to dx/dt (the velocity of the interface, correct me if I am mistaken). I think you want to divide dx/dt by an appropriate length scale for your vapor-liquid interface transport phenomena, then the source term can be added to the TEqn with correct units.


It seems simple but in my experience careful dimensional analysis can solve/avoid many issues when writing equations for a finite volume solver.


Cheers

Hello Luca,
Thank you for your comments. Energy balance will be written in terms of W/m3 in finite volume approach basically. My problem is a kind of 1D pyrolysis (Coal) modelling, which is focussed mainly on heat loss - thus only Energy equation is included in my work.

ENERGY EQU: rho*Cp*(dT/dt) = del/delx[K*dT/dx ] + rho*dQ/dT + r*cp*dT/dx (W/m3)
Unsteady conduction = Diffusion term + SOURCE TERM 1 + SOURCE TERM 2
SOURCE TERM 1 = heat loss due to pyrolysis
SOURCE TERM 2 = phase change of moisture (moisture embedded in the wet coal)

My query is all about the SOURCE TERM 2. Concerning it, I posted my query here - because this post discusses phase change problem - thought of gaining some assistance here for my work.
Pyrolysis - Arrhenius-like degradation chemistry
Drying - surface modelling at boiling plane (drying when moisture boils at 100deg)
Since, drying is based on surface modelling, the mass and heat balances unit is calculated at the surface. The calculated mass balance (r) is multiplied with cp and dT/dx as SOURCE TERM 2 gives the unit of W/m3.
REF - Clear details here in this MANUSCRIPT ~ https://sci-hub.tw/https://doi.org/1...361(83)90225-9 (ATKINSON, B., & MERRICK, D. (1983). Mathematical models of the thermal decomposition of coal4. Heat transfer and temperature profiles in a coke-oven charge. Fuel, 62(5), 553–561)
How the source terms can be included in such problems? I have contacted certain people, few replied as they are not familier with pyrolysis, others asked me to look into interFoam solvers. But the approach of interFoam seems to different with separate transport equation for volume fraction, which is not the same in my case. I'm trying to figure out and learn as how the problem can be attacked based on algorithm and importantly how such problems can be approached step by step in OpenFOAM//

Please share your ideas, it will be highly helpful.
Thank you.

I should preface by saying I am relatively new to OpenFOAM and am currently working on a cht problem with solidification and melting, hence my monitoring of this thread. I do have more extensive experience with CFD modeling in general.



Based on your last post I don't see a problem. Perhaps I am misunderstanding, but it looks as if you have a valid energy equation in W/m^3. The 2d surface model is an additional function used to define your source term, the units of this don't necessarily matter so long as they are correct within the function definition. Do you know dx/dt based on results from the previous timestep? The way I see it, you should be able to calculate (r) with a function object, but I am not sure if you could then use that term to your TEqn.

Hello Luca,

Sorry for my late response.
Not sure, I'm looking forward to find a way for it. Yae, dx/dt are based on previous time steps. Function objects usually calculates post processing output results right (also using swak4Foam sometimes). Correct me if Im wrong
Thank you

Moisture evaporation - if-else loop
 

Hello Luca,
In previous message, I have given a condition for moisture evaporation (phase change) quoted in a box.I tried to evaluate the above condition separately. I mean, without pyrolysis, only phase change, in order to understand that how such condition can be implemented in the OpenFOAM.
I'm hereby attaching the solver (OF211). In this solver, under TEqn.H, I solved a if-else loop based on the above condition.However, the alpha (moisture content) is not solved properly with appropriate outcome. In my case, alpha has no expression (field operation), alpha is only solved based on the loop. So, do we need an expression to solve alpha basically? any ideas ?

If you have any ideas about it, kindly do share.

It will be helpful. Thank you

Solidification of pure metal
 

I was trying to make solidification model of pure metal.
can I use erfConvectiveMeltingPimpleFoam solver for solidification just by changing initial and boundary conditions or do I need to make any changes to the solver part?

making appropriate required physical changes should do

Sir, can I know what are all the changes to be done?

first check your solver for any test case.. if its working fine then all you need to do is supply physical BC & IC along with thermophyscial parameter..

Sir, i have checked the solver for melting problem and it's working fine. But, when it come to Solidification problem, i'm unable to capture the natural convection effects. so, could you please help me with the initial and boundary conditions.

post your results and error ...

Hi everyone,
I need to simulate the solidification of steel.
At the moment I'm using the solver solidificationFoam, but there are some missing things:
-Communication between solid and liquid domain, ingot and ingot mold
-Thermal Shrinkage
-Viscosity variation as a function of temperature

What solver do you recommend?

Thanks in advance!

meltingFoam works fine in OpenFOAM (v1912)
 

Thanks Pavithra for sharing the solver meltingFoam. It works fine in the present version of OpenFOAM (v1912).

Regards

please
 

Can i have the paper？ Thanks！！！

request
 

Hello Som.
I want to cool a heat sink by butting it on the wax, just like your problem.
May I have your code?
would you please send it to my email?
aerospace1365@gmail.com

request
 

Dear Sir
I have a case that should be solved by your solver. I have a heat sink (fins) located in wax. the wax will melt and cool the fins.
do you think this solver can solve this problem?
my Gmail is aerospace1365@gmail.com.
I am looking forward to hearing from you.
regards

OF8 solidificationMeltingSource with erf
 

Hello all,

I want to port the solidification solver to OF8 to utilize the tabulated thermophysicalProperties in my bachelor's thesis. "solidificaitionMeltingSource" looked like a good starting point to me, so I changed the "update" and "apply" functions according to the formulation in Fabian Röslers Solver and made some minor changes like deleting alpha1e and relax factor. The code is attached in the tar file.

---
Question 1: Is there missing a negative sign in the paper in equation 10? In the uploaded solver post 19 there is an additional negative sign in the upper numerator of TEqn, as noted in the equation snippet below. Implementation according to the paper throws a floating point exception error.


---

I ran some cases but the results are terrible as can be seen in the attached picture (as in the paper the timesteps 2, 6, 10 and 17 min are plotted). The melt front is way too fast and doesn't show curving like the results from the paper. Out of curiosity I tested the standard implemented linear fvOption with some cranked up numerics and ran the same settings again for the erf fvOption:

• linear fvOption differs from the results in the paper in melt front position and shape increasingly with time. Adjusting the model parameters helps a bit in later time steps with the shape but the initial deviation stays the same.
• erf fvOptoin with cranked up numerics moves the melt front a bit in the right direction, but is far from accurate.

---
Question 2: Does anybody have an idea where this behavior comes from? I even programmed it a second time from the scratch to avoid typos, but the results were the same. And this doesn't explain the behavior of the linear model as well.
---

For post processing of the melt front run "pvpython writeMeltFront.py" in the case directory and point plotMeltFront.py to meltFrontPosition directory in corresponding case directory. I left the results in the myCases cases for an easy first impression.

Thanks in advance and have a nice day!

Hello Kai,


Regarding your first question: the derivative of is , so I think yes there should be a minus in equation 10.

And regarding your second question. The erf-Solver is rather mediocre. Even if everything is implemented correctly you will need a very small time step, since the very strong non-linearity of the equation is not handled properly. You can take a look at: https://www.mdpi.com/1996-1073/12/5/868 for a better implementation.


Best,
Moritz

Hello again forum!

Moritz posted an interesting approach, I will look further into it. But in the mean time, does anyone of you have an Idea where especially the deviation of the linear model comes from?

Thanks in advance!

Hi,
This is my first post since years. The erf function was the first more or less running version we developed at the thermodynamics chair in Bayreuth. The more sophisticated version is in my PhD thesis. The latest version seems to be the one in Moritz's article. I had no time to reed through the article yet and I'm out of this topic since 7 years :).
What deviation do you mean?
Cheers
Fabian

Hi Fabian,

thanks for your reply, I received your PhD thesis yesterday and it looks very helpful from what I've read until now! Especially the parts about sinking solid and literature review.

The deviation I'm talking about are the more narrow expansion in x of the linear model compared to yours and Brent and Voller's simulations and what is the most significant difference in my opinion is the shape of the phase front. The more "wave-like" character up until 360 s and the pronounced bulge from then on.

But I missed a paper Moritz sent me in an E-Mail which explains these results, shame on me :eek: Hannoun et al. published "Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side" and state, that Brent and Voller's results are not grid-converged. So the positions of vortexes are wrong and in the early stages of the simulation, when vortexes are small, they are not resolved at all.
Further the experiment of Gau and Viskanta drained the fluid portion to measure the phase position, so we have different bc's in the experiment and in the simulation.

I will look into your thesis for the theory behind phase change but most likely adopt Moritz's approach because it is the most advanced and computational cost may be reduced.

kOmegaSSTIDDES in OpenFOAM-v2012(pisoFoam)
 

Hi:

Has anyone encountered the following problems？
I guess the boundary conditions are set incorrectly? But I have tried many, still have this problem.



--> FOAM FATAL ERROR: (openfoam-2012)
Attempt to cast type calculated to type nutWallFunction at index 2

From To& Foam::refCast(From&, Foam::label) [with To = const Foam::nutWallFunctionFvPatchScalarField; From = const Foam::fvPatchField<double>; Foam::label = int]
in file /home/user2/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/lnInclude/typeInfo.H at line 162.

FOAM aborting

#0 Foam::error::printStack(Foam::Ostream&) at ??:?
#1 Foam::error::exitOrAbort(int, bool) at ??:?
#2 Foam::nutWallFunctionFvPatchScalarField::nutw(Foam ::turbulenceModel const&, int) at ??:?
#3 Foam::omegaWallFunctionFvPatchScalarField::calcula te(Foam::turbulenceModel const&, Foam::List<double> const&, Foam::fvPatch const&, Foam::Field<double>&, Foam::Field<double>&) at ??:?
#4 Foam::omegaWallFunctionFvPatchScalarField::calcula teTurbulenceFields(Foam::turbulenceModel const&, Foam::Field<double>&, Foam::Field<double>&) at ??:?
#5 Foam::omegaWallFunctionFvPatchScalarField::updateC oeffs() at ??:?
#6 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::Boundary::updateCoeffs() at ??:?
#7 Foam::kOmegaSSTBase<Foam::LESModels::DESModel<Foam ::IncompressibleTurbulenceModel<Foam::transportMod el> > >::correct() at ??:?
#8 ? at ??:?
#9 __libc_start_main in /lib64/libc.so.6
#10 ? at ??:?
Aborted (core dumped)



Thanks in advance

Thesis
 

Hello Fabian, is it okay to have a copy of your thesis. Would it be better if it is written in English but if not, then I will try to find a way to translate it. Thanks

I have solved this problem, my boundary conditions are set incorrectly, mainly because I have to be careful about the setting of inlet, outlet and wall.

Question about the parallel calculation
 

Hi everyone
I want to cpuple interfoam and the solidification model to simulate the droplet freezing. Thus I put the T equation in the interFoam solver and compiled. The solver works fine in serial mode. But when I run in parallel. There is error:
//////////////////////////////////////////////////////////////////////////////
Fatal error in MPI_Recv: Message truncated, error stack:
MPI_Recv(224).....................: MPI_Recv(buf=0x1da60c0, count=6400, MPI_BYTE, src=8, tag=1, comm=0x84000002, status=0x7ffec83b2e40) failed
MPIDI_CH3U_Receive_data_found(131): Message from rank 8 and tag 1 truncated; 19200 bytes received but buffer size is 6400
Fatal error in MPI_Recv: Message truncated, error stack:
MPI_Recv(224).....................: MPI_Recv(buf=0x18748a0, count=6400, MPI_BYTE, src=10, tag=1, comm=0x84000002, status=0x7ffe00c99a10) failed
MPIDI_CH3U_Receive_data_found(131): Message from rank 10 and tag 1 truncated; 19200 bytes received but buffer size is 6400
Fatal error in MPI_Recv: Message truncated, error stack:
MPI_Recv(224).....................: MPI_Recv(buf=0x236adb0, count=6400, MPI_BYTE, src=11, tag=1, comm=0x84000002, status=0x7ffd32387010) failed
MPIDI_CH3U_Receive_data_found(131): Message from rank 11 and tag 1 truncated; 19200 bytes received but buffer size is 6400
Fatal error in MPI_Recv: Message truncated, error stack:
MPI_Recv(224).....................: MPI_Recv(buf=0x2a4f080, count=6400, MPI_BYTE, src=9, tag=1, comm=0x84000002, status=0x7ffd606c04e0) failed
MPIDI_CH3U_Receive_data_found(131): Message from rank 9 and tag 1 truncated; 19200 bytes received but buffer size is 6400
//////////////////////////////////////////////////////////////////////////////////
the attachmen is the T Equation file. Can someone know the reason why this error happens

Error in parallel runnnig
 

Hello, everyone, have you succeeded in running cases in parallel when using the fabian_roesler code. I have followed fabian_roesler's code and it is successful in serial running but failed in paralleling running. I do not know the reason. Have someone met and fixed this problem?