[Kummi]

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]

[Kummi]

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.

Mass balance (r) = rho * alpha * (dx/dt)
Heat balance (-K dT/dx) = rho * latent heat
There is only TEqu.H is my case, then how to define (dx/dt) velocity. If I consider phi as velocity, then phi is surfaceScalarField interpolated from U.
So, how to correct the mass balance based on velocity (dx/dt)?

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 ==>

@T=100deg (moiture content_w=0)

>> Tl (liquid) = Tv (vapor) [Looks like single point phase change scheme over the specified temperature by Voller and Cross]
>> H2O(l) = H2O (v)

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 !!!

[Kummi]

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

alpha = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5);

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:

Initial moisture content alpha = 0.1 (10%)
@BOILING TEMP, T=100deg
(moiture content) alpha=0 (near the heated wall)

//************************************************** *****************//
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,

Mass balance (r) = rho * alpha * (dx/dt)
Heat balance (-K dT/dx) = (r) * latent heat
Heat balance expressed as W/m2 --> looks like surface tracking problem based on STEFAN 1D phase change CONDITION

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

[PositronCascade]

Is there any CHT melting solver that you know?


Thanks!

[Pavithra]

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.

[Kummi]

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.

But, when I changed the mesh resolution to 132X100 in the blockMesh dictionary.

No such problems even if I changed the mesh size.
You have compiled in OF V6 - thats where the problem is, I hope so.

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

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 ^^

[Pavithra]

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.

[Kummi]

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 ?

[Pavithra]

Originally Posted by Kummi

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

[Kummi]

Hello Pavithra,

So, in normal case, cellToRegion file does not have any use.

- Yes, no need of cellToRegion here.

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

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

[Pavithra]

Originally Posted by Kummi

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.

[Pavithra]

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.

[Pavithra]

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.

[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

[luks1910]

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

[Kummi]

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.

CONDITION: T= 100 ==> moisture content =0
Based on it, mass and heat balances are calculated at the boiling plane (surface) as,
Mass balance (r) = rho * alpha * (dx/dt) [kg/m2.s]
Heat balance (-K dT/dx) = (r) * latent heat [W/m2]

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.

[luks1910]

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.

[Kummi]

Hello Luca,

Sorry for my late response.

The 2d surface model is an additional function used to define your source term

Not sure, I'm looking forward to find a way for it.

Do you know dx/dt based on results from the previous timestep?

Yae, dx/dt are based on previous time steps.

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.

Function objects usually calculates post processing output results right (also using swak4Foam sometimes). Correct me if Im wrong
Thank you

[Kummi]

Hello Luca,
In previous message, I have given a condition for moisture evaporation (phase change) quoted in a box.

CONDITION: T= 100 ==> moisture content =0
Based on it, mass and heat balances are calculated at the boiling plane (surface) as,
Mass balance (r) = rho * alpha * (dx/dt) [kg/m2.s]
Heat balance (-K dT/dx) = (r) * latent heat [W/m2]

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.

forAll(mesh.cells(), celli)
{
if (T[celli] >= Tsat.value())
{
alpha[celli] = scalar(0.0);
}
else
{
alpha[celli] = alphaTemp[celli];
// alpha[celli] = alpha[celli];
// alpha[celli] = alpha[nearestCellIndex];
}
}

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

[dp32]

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?