CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Solid/liquid phase change

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   May 13, 2009, 09:35
Default Solid/liquid phase change
  #1
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 161
Rep Power: 8
fabian_roesler is on a distinguished road
Dear Foamers

For a project I am trying to get a solid/liquid phase change simulation running in OF 1.5-dev. A paraffin with non isothermal phase change is molten in a 2D rectangular test case scenario. I started from interFoam solver by changing especially the gamma equation to my needs:

Gamma equation:

dimensionedScalar pi = mathematicalConstant:i;
dimensionedScalar alpha = (Tl-Ts)/(scalar(2)*tan(scalar(0.49)*pi));
gamma = atan(scalar(1)/alpha*(T-(Tl+Ts)/scalar(2)))/pi+scalar(0.5);

Changig gamma from 0 (solid) to 1 (liquid) depending on the temperatures
Tl (upper phase change temperature) Ts (lower phase change temperature)

Momentum equation:

fvVectorMatrix UEqn

(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
);

UEqn.relax();

if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}

No changes at all, however no surface tension left

Energy equation:

volScalarField rhoCp("rhoCp", rho*cp);
surfaceScalarField rhoPhiCp("rhoPhiCp", rhoPhi*fvc::interpolate(cp));

fvScalarMatrix hEqn
(
fvm::ddt(rhoCp, T)
+ fvm::div(rhoPhiCp, T)
- fvm::laplacian(lambda, T)
);

solve
(
hEqn ==
- hs*fvc::ddt(rho, gamma) //Source term from Stefan condition
);

The source term derived from the Stefan condition takes care of the temperature gradient at the interface induced by the phase change.

Solver:

Like in interFoam, the solver works as follows:
  1. Solving gamma equation
  2. Solving coupled pressure and velocity equations
  3. NEW solving energy equation NEW
The solver is very unstable what, in my opinion is due to the couppling of gamma and energy. Please correct me when I am wrong. My questions are:

Is it really the coupling of gamma and energy, making the solver instable?
Is there a way to solve the coupled equations (relaxation)?

Thanks allot for your replies.

Best regards

Fabian
fabian_roesler is offline   Reply With Quote

Old   May 18, 2009, 09:51
Default Simulation stops
  #2
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 161
Rep Power: 8
fabian_roesler is on a distinguished road
Hallo Foamers

I managed to get my model running. I substituted the gamma equation into the energy conservation equation and it worked. No I encountered another problem. The solver runs fine until a certain time, where all iteration numbers equal zero and the results do not change anymore. The solver continues running until endTime but nothing changes anymore. Here are the two time steps from the log file where the error occurred in between:

PBiCG: Solving for T, Initial residual = 1.00244e-05, Final residual = 8.71675e-11, No Iterations 1
Liquid phase volume fraction = 0.00105084 Min(gamma) = 0.000769482 Max(gamma) = 0.998571
PBiCG: Solving for Ux, Initial residual = 0.000224831, Final residual = 1.0262e-16, No Iterations 1
PBiCG: Solving for Uy, Initial residual = 0.000162936, Final residual = 7.40847e-17, No Iterations 1
PCG: Solving for pd, Initial residual = 0.000982493, Final residual = 1.25233e-05, No Iterations 2
PCG: Solving for pd, Initial residual = 1.25111e-05, Final residual = 4.69217e-07, No Iterations 23
PCG: Solving for pd, Initial residual = 4.69216e-07, Final residual = 7.52164e-08, No Iterations 25
PCG: Solving for pd, Initial residual = 7.55978e-08, Final residual = 7.55978e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55978e-08, Final residual = 7.55978e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55978e-08, Final residual = 7.55978e-08, No Iterations 0
time step continuity errors : sum local = 2.05188e-15, global = 1.44131e-30, cumulative = -4.99715e-30
ExecutionTime = 287.72 s ClockTime = 288 s

Courant Number mean: 2.76148e-12 max: 7.1681e-10 velocity magnitude: 2.51386e-11
deltaT = 0.00990099
Time = 11.7822

PBiCG: Solving for T, Initial residual = 9.99292e-06, Final residual = 9.99292e-06, No Iterations 0
Liquid phase volume fraction = 0.00105084 Min(gamma) = 0.000769482 Max(gamma) = 0.998571
PBiCG: Solving for Ux, Initial residual = 1.10625e-12, Final residual = 1.10625e-12, No Iterations 0
PBiCG: Solving for Uy, Initial residual = 6.0714e-13, Final residual = 6.0714e-13, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55979e-08, Final residual = 7.55979e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55979e-08, Final residual = 7.55979e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55979e-08, Final residual = 7.55979e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55979e-08, Final residual = 7.55979e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55979e-08, Final residual = 7.55979e-08, No Iterations 0
PCG: Solving for pd, Initial residual = 7.55979e-08, Final residual = 7.55979e-08, No Iterations 0
time step continuity errors : sum local = 2.05188e-15, global = -4.11243e-31, cumulative = -5.40839e-30
ExecutionTime = 287.81 s ClockTime = 288 s

Courant Number mean: 2.76148e-12 max: 7.16811e-10 velocity magnitude: 2.51386e-11
deltaT = 0.00990099
Time = 11.7921

I think all residuals are OK and the Courant number is way to small to have an effect. Any ideas what could cause this problem? Thanks in advance for your replies.

Regards
fabian_roesler is offline   Reply With Quote

Old   June 2, 2009, 12:00
Default
  #3
Member
 
Rachel Vogl
Join Date: Jun 2009
Posts: 48
Rep Power: 8
Rachel is on a distinguished road
Hi Fabian,

My guess is that the reducing the number of pd iterations might solve the issue. I guess its the non-orthogonal pressure corrections.

In case this does not work, consider reducing the residual to say e-10 or relTol parameters.

Regards
Rachel
Rachel is offline   Reply With Quote

Old   June 3, 2009, 03:37
Smile PISO in multiphase phase change problems
  #4
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 161
Rep Power: 8
fabian_roesler is on a distinguished road
Hi Rachel

Thank you for your reply. You are right: It's the residual causing this problem. The initial residual

PBiCG: Solving for T, Initial residual = 9.99292e-06, Final residual = 9.99292e-06, No Iterations 0

is allready below the one beeing requested by fvSolution. I set the residual to 1e-10 and the simulation runs without problems.

However I have still some problems with continuity. Like in the cavitation multiphase solvers, the PISO pressure equation has to be changed in a certain way. In interPhaseChange Foam for example, pEqn is changed from

fvScalarMatrix pdEqn
(
fvc::div(phi) - fvm::laplacian(rUAf, pd)
);

to

fvScalarMatrix pdEqn
(
fvc::div(phi) - fvm::laplacian(rUAf, pd)
+ (vDotvP - vDotcP)*(rho*gh - pSat) + fvm::Sp(vDotvP - vDotcP, pd)
);

I have to have a closer look on the PISO scheme and the pressure equation to learn how to derive the equation suited to my model.
Does anybody have some tips on how to use PISO in phase change multiphase problems? Thanks allot in advance.

Regards

Fabian
fabian_roesler is offline   Reply With Quote

Old   June 4, 2009, 10:53
Default substitution of equation into the energy conservation equation
  #5
Member
 
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 8
haghajani is on a distinguished road
Dear Fabian,

Would you please let me know, how do you substituted the gamma equation into the energy conservation equation. Please drop some lines, how this could solve the problem arose?

Thanks in advance,

Hamed
haghajani is offline   Reply With Quote

Old   June 5, 2009, 06:18
Default
  #6
New Member
 
Join Date: May 2009
Posts: 7
Rep Power: 8
dakos is on a distinguished road
Hi Fabian,

I'm also interested in solidification simulation, but i have some problems in defining gamma as a function of temperature.
Can you tell me how you did?

Dario
dakos is offline   Reply With Quote

Old   June 5, 2009, 06:45
Default Continuous gamma function
  #7
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 161
Rep Power: 8
fabian_roesler is on a distinguished road
Hi Hamed, hi Dario

My gamma equation is in contrast to the literature a continuous function. The literature defines gamma as 0 in the solid phase, 1 in the liquid phase and linear in between. I am using a arcus tangent function.

dimensionedScalar alpha = (Tl-Ts)/(scalar(2)*tan(scalar(0.49)*pi));
gamma = atan(scalar(1)/alpha*(T-(Tl+Ts)/scalar(2)))/pi+scalar(0.5);

The source term in the energy equation

- hs*fvc::ddt(rho, gamma)

can be written as

- hs*
(
rho*fvc::ddt(gamma)
+ gamma*fvc::ddt(rho)
);

the time derivative of gamma can be directly included here so the two coupled equations form one equation. The gamma field can be calculated after solving the energy equation, however this is not necessary for the solver any more.
Anybody can tell me more about the PISO scheme in multiphase solvers?

Regards Fabian
fabian_roesler is offline   Reply With Quote

Old   May 27, 2010, 10:11
Default
  #8
Member
 
Jitao Liu
Join Date: Mar 2009
Location: Jinan , China
Posts: 64
Rep Power: 8
awacs is on a distinguished road
Quote:
Originally Posted by fabian_roesler View Post
Hi Hamed, hi Dario

My gamma equation is in contrast to the literature a continuous function. The literature defines gamma as 0 in the solid phase, 1 in the liquid phase and linear in between. I am using a arcus tangent function.

dimensionedScalar alpha = (Tl-Ts)/(scalar(2)*tan(scalar(0.49)*pi));
gamma = atan(scalar(1)/alpha*(T-(Tl+Ts)/scalar(2)))/pi+scalar(0.5);

The source term in the energy equation

- hs*fvc::ddt(rho, gamma)

can be written as

- hs*
(
rho*fvc::ddt(gamma)
+ gamma*fvc::ddt(rho)
);

the time derivative of gamma can be directly included here so the two coupled equations form one equation. The gamma field can be calculated after solving the energy equation, however this is not necessary for the solver any more.
Anybody can tell me more about the PISO scheme in multiphase solvers?

Regards Fabian

Hi Fabian,

I am trying to simulate melt filling process in injcetion molding. The polymer melt solidifies and stops flowing at cavity wall due to the cooling effect of cold cavity. As a result, the melt/air and melt/solid interfaces need to be fixed in simulation. I have added the energy equation into interFoam. The non-newtonian viscosity is take into account, too. The problem is that how can I determine the varied melt/solid interface.

Please give me some tips. Thanks in advance.

Best regards,
Jitao
awacs is offline   Reply With Quote

Old   May 28, 2010, 08:45
Smile Gamma function for injection molding
  #9
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 161
Rep Power: 8
fabian_roesler is on a distinguished road
Hi Jitao

the phase change is taken into account by adding a source term to the energy conservation equation. Internal energy is no longer

cp*T

but

cp*T + gamma*hs

where gamma is the phase fraction and hs is the latent heat of fusion. The method is called the enthalpy method. It can be used when phase change occurs over a temperature range and without crystallisation effect.

Regards Fabian
fabian_roesler is offline   Reply With Quote

Old   September 1, 2010, 18:01
Default
  #10
Member
 
Join Date: Nov 2009
Posts: 48
Rep Power: 7
farhagim is on a distinguished road
Hi Fabian,

First of all Thank you for sharing your information. Second, I have some question regarding the solidification (Ice accretion ) with interFoam.
1- does your mesh update in each iteration?
2- have you managed to get good result?
3- can you explain it to how does it works, I mean I did not see any condition in your solid -liquid phase for temperature. I guess you involve the condition in your gamma, am I right??
Finally would you please give me some hint or any progress in your implementation?


Thanks

Mehran

Quote:
Originally Posted by fabian_roesler View Post
Hi Jitao

the phase change is taken into account by adding a source term to the energy conservation equation. Internal energy is no longer

cp*T

but

cp*T + gamma*hs

where gamma is the phase fraction and hs is the latent heat of fusion. The method is called the enthalpy method. It can be used when phase change occurs over a temperature range and without crystallisation effect.

Regards Fabian
farhagim is offline   Reply With Quote

Old   December 24, 2012, 07:37
Default
  #11
New Member
 
mamadreza
Join Date: Mar 2011
Posts: 22
Rep Power: 6
mamyjooooon is on a distinguished road
hi
Have you succeeded in solidification in OpenFoam?
I am facing a similar problem and I want to model the aircraft in-fight ice accretion using openfoam or fluent
please help me about this problem
ThankS
mamyjooooon is offline   Reply With Quote

Reply

Tags
coupling, interfoam, phase change, piso, solid/liquid

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Thermal phase change model Piti CFX 0 February 9, 2009 10:27
Two phase flow with phase change Ahmad Al-Zoubi CFX 1 November 26, 2008 04:59
How to model the solid-fluid phase change in CFX ohrmond CFX 2 May 26, 2006 06:27
thermal phase change question CFDflying CFX 1 February 18, 2004 05:10
compressible two phase flow in CFX4.4 youngan CFX 0 July 1, 2003 23:32


All times are GMT -4. The time now is 04:41.