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

melting problem

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

Like Tree74Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 22, 2013, 15:12
Default History excurse
  #61
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
Hi

Well it is the enthalpy method. You use one eulerian grid for both, fluid and solid PCM. Oleinik [1] proposed a weak formulation of the Stefan problem in the 1960s. Based on this approach, Rose [2] and Solomon [3] developed a method based on enthalpy. This is the mother of all enthalpy methods developed since.
[1] O. A. Oleinik: A method of solution of the general Stefan problem, Soviet Mathematics, vol. 1, p. 1350–1354, 1960.
[2] M. E. Rose: A Method for Calculating Solutions of Parabolic Equations with a Free Boundary, Mathematics of Computation, vol. 14, p. 249–256, 1960.
[3] A. D. Solomon: Some Remarks on the Stefan Problem, Mathematics of Computation, vol. 20, p. 347–360, 1966.

Hope this is what you were searching for.

Regards

Fabian


fabian_roesler is offline   Reply With Quote

Old   August 28, 2013, 04:42
Default
  #62
Member
 
Thomas Vossel
Join Date: Aug 2013
Location: Germany
Posts: 45
Rep Power: 4
ThomasV is on a distinguished road
Thanks for your answer...

I guess I was a bit confused with the model and things like the buoyancy term and thus thought that there was some kind of multiphase-flow going on. As the model simplifies things into a one-phase flow there is nothing like that though...
ThomasV is offline   Reply With Quote

Old   September 12, 2013, 05:53
Default
  #63
New Member
 
Eric Pernot
Join Date: Apr 2013
Posts: 1
Rep Power: 0
EPernot is on a distinguished road
Hi everybody,

I would like to update the Rosler's solver to the version 2.2.0 of O.F.
Do you know what kind of modifications I have to made for that?

Thanking you in advance.

Eric
EPernot is offline   Reply With Quote

Old   September 19, 2013, 04:53
Default
  #64
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi,
I have an stupid question about post processing. I can plot the alpha (phase fraction) contour in paraview. But I'm wondering how to plot the graph (and export) for comparing as it is not possible in Paraview.
ahmmedshakil is offline   Reply With Quote

Old   September 19, 2013, 05:21
Default
  #65
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
Hello Mohammad,

well, first of all, there are no stupid questions :-).

A way I use is only to enable alpha to be shown > in Volume Fields. Then I do a "contour" like you with only one value, 0.5 for example. Split the view and click "Spreadsheet View" for the new one. For the spreadsheet view, only the eye in front of "contour" should be enabled. With the spreadsheet view still active, do "Save Data" and choose "csv" as file type, for example. You will obtain a column file with x,y,z,alpha (not necessarily in that order) that you can plot with gnuplot or other tools to compare the results.

If your cases do not vary in geometry, you could compare even a distorted contour plane in 3D using paraview. Only open both cases, create the contour for each case, colour them differently and let them be shown (active eye") in a single 3D view.

I hope, this does help you
Regards, Anja
apple-tree likes this.
AnjaMiehe is offline   Reply With Quote

Old   September 20, 2013, 04:50
Default
  #66
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Thanks AnjaMiehe.
ahmmedshakil is offline   Reply With Quote

Old   December 3, 2013, 11:23
Default Voller Source based enthalpy method
  #67
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi,
I'm trying to update the solid-liquid fraction as mentioned by Voller source based enthalpy method. First, for checking the case I have used just heat conduction. The way I'm trying to solve the energy equation is:

volScalarField gamma0 = gamma; //-storing the prev time value

do
{
gamma.storePrevIter();
volScalarField dFdT = 1/(Tl-Ts);
volScalarField Finv = gamma*(Tl-Ts)+Ts;
volScalarField Sp = -((rho*L)/runTime.deltaTValue())*dFdT;
volScalarField Sc = -Sp*Finv + ((rho*L)/runTime.deltaTValue())*(gamma0-gamma);

fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T) -fvm::laplacian(K,T) - fvm::Sp(Sp, T) = Sc

);
TEqn.solve();

gamma = gamma + dFdT*(T-Finv);
gamma = gamma = min(max(gamma,scalar(0)),scalar(1));
Cp = gamma*Cpl+(1-gamma)*Cps;
K = gamma*Kl +(1-gamma)*Ks;
//- residual calculation of g

}while (condition for convergence of g)

This the main code for energy equation, and implemented exactly the same case of Gallium error melting just for Test. I'm having the PROBLEM as,
PROBLEM
(1) Temperature is not updating with time advancing
(2) I have switch off (Sp and Sc )to check from where the problem is coming, I saw it's coming from Sp and Sc.
(3) I'm wondering whether the loop has any problem ???
I have attached the solver.....
PLEASE ADVICE ME TO COME OUT OF THE PROBLEM
Attached Files
File Type: gz enthalpyFoam.tar.gz (4.7 KB, 24 views)
ahmmedshakil is offline   Reply With Quote

Old   December 3, 2013, 11:49
Default
  #68
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
Hello Mohammad,
I don't have time right now to look into it, i will try it tomorrow.

However, there are some things that make me think about it already now:
- Could you provide the equations you based your code on here. I don't have the paper at hand an Voller wrote lots of them ;-)
- for dFdt I suppose a linear slope of the fraction of solid, right? But even thought, this is not 1/(Tl-Ts) everywhere, only within the mushy zone. Thus:
Code:
 dFdt=pos(Tl-T)*pos(T-Ts)/(Tl-Ts)
where pos(Tl-T) is 1 when Tl-T is positive, that is when T is smaller than Tl, and the same reasoning for the second pos(..) factor.
- The old time value is used gamma.oldTime(), no need to store it in an additional variable

So much for now,
Regards, Anja
AnjaMiehe is offline   Reply With Quote

Old   December 4, 2013, 01:18
Default
  #69
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi AnjaMiehe,
I have attached the equation that I am solving.
The paper that I am working on is: http://www.tandfonline.com/doi/abs/1...2#.Up63WkAW1hg
Attachment 27158

Quote:
Originally Posted by AnjaMiehe View Post
Hello Mohammad,
I don't have time right now to look into it, i will try it tomorrow.

However, there are some things that make me think about it already now:
- Could you provide the equations you based your code on here. I don't have the paper at hand an Voller wrote lots of them ;-)
- for dFdt I suppose a linear slope of the fraction of solid, right? But even thought, this is not 1/(Tl-Ts) everywhere, only within the mushy zone. Thus:
Code:
 dFdt=pos(Tl-T)*pos(T-Ts)/(Tl-Ts)
where pos(Tl-T) is 1 when Tl-T is positive, that is when T is smaller than Tl, and the same reasoning for the second pos(..) factor.
- The old time value is used gamma.oldTime(), no need to store it in an additional variable

So much for now,
Regards, Anja
Attached Images
File Type: jpg Voller.jpg (39.6 KB, 102 views)
ahmmedshakil is offline   Reply With Quote

Old   December 4, 2013, 07:51
Default
  #70
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
Hello Mohammad,

to begin with, I am not quite sure if you will like my solution. I have always understood OpenFOAM to code in using grad, div, laplacian ... and not trying to put the discrete equation together with a_p ... Using the Voller paper and your equations I did the following:
\rho\,c_p\,\frac{\partial T}{\partial t}= \nabla \left( K\nabla T\right) + S
with
S=-\delta\,H\frac{\partial g}{\partial t}\quad \text{ and } \quad \delta\,H=\rho\,\left(C_{pl}-C_{ps} \right) T + \rho\,Lf.
Knowing that
\frac{\partial g}{\partial t} = \frac{\partial g}{\partial T} \frac{\partial T}{\partial t} \quad \text{ and } \quad g=\frac{T-T_s}{T_l-T_s} \,\forall\, T_l > T > T_s \quad \text{ gives } \quad  \frac{\partial g}{\partial t} = \frac{1}{T_l-T_s} \frac{\partial T}{\partial t}
the source term S can be put together in the following way
S=-\frac{\rho}{T_l-T_s}\left[\left(C_{pl}-C_{ps}\right)T+Lf\right]\frac{\partial T}{\partial t} \,\forall\, T_l > T > T_s.
The equation to be coded is
\rho\,\left( C_p+ A\right) \frac{\partial T}{\partial t} =  \nabla \left( K\nabla T\right)
with
A= \begin{cases}  \frac{1}{T_l-T_s}\left[\left(C_{pl}-C_{ps}\right)T+Lf\right]  & \,\forall\, T_l>T>T_s \\      0 \quad & \text{else}      \end{cases}
The equation is coded in the following way
Code:
fvScalarMatrix TEqn
(
    (Cp+pos(Tl-T)*pos(T-Ts)*((Cpl-Cps)*T+Lf)/(Tl-Ts))*rho0*fvm::ddt(T)
    - fvm::laplacian(fvc::interpolate(K),T)
);
followed by the update of gamma
Code:
gamma = (T-Ts)/(Tl-Ts);
gamma = min(max(gamma,scalar(0)),scalar(1));
gamma.relax();.
Please find attached the code and a Stefan test case. For the test case, you need a working swak4Foam installation and gnuplot for the final comparison.

Sorry, that I cannot assist you with the exact implementation of the Voller paper. If you want to move on to more complex correlations of gamma to temperature, you could try to read in an interpolation table. See this link for more information: Read temperature dependent thermophysical properties from a file - boundaries false

Regards, Anja
Attached Files
File Type: gz enthalpy2Foam.tar.gz (2.9 KB, 45 views)
File Type: gz Stefan_Test.tar.gz (46.1 KB, 41 views)
mm.abdollahzadeh likes this.
AnjaMiehe is offline   Reply With Quote

Old   December 4, 2013, 08:32
Default
  #71
Member
 
Thomas Vossel
Join Date: Aug 2013
Location: Germany
Posts: 45
Rep Power: 4
ThomasV is on a distinguished road
Well I'm currently writing my masters thesis about this topic and I got something running that works (although has to be revamped at certain places to include additional features). My problem is about the flow and solidification of an alloy though but as I'm right now changing the temperature calculations in order to correctly take care of the eutectic solidification I can post how I'm doing things right now. My solver's design also is based on Voller's enthalpy-porosity-model but I'm not really copying a certain paper 1:1...

But oh well here is the code of my TEqn.H file which calculates the new temperature field:

Code:
// initialize and solve energy equation
    fvScalarMatrix TEqn
    (
        cpEff * fvm::ddt(Temp)
      + fvm::div(phib * fvc::interpolate(cpMean), Temp, "div(phib*cpEff,Temp)")
 
      - fvm::laplacian(lambda / rhob, Temp, "laplacian(lambda|rhob,Temp)")
    );

    TEqn.relax();
    TEqn.solve();
I guess the variables should be pretty much self explanatory. In addition to that I'll give you my code for calculating the cp values and how my data based fraction solid calculation method look like (which currently uses a fake function which is going to be replaced with actually measured data)...

Code:
forAll(fracSol, celli)
    {

    // update cpMean, cpEff and lambda
        if (Temp[celli] < TbSol.value() )
        {
            cpMean[celli] = cpbSol.value();
            cpEff[celli] = cpMean[celli];

        }

        else if (Temp[celli] > TbLiq.value() )
        {
            cpMean[celli] = cpbLiq.value();
            cpEff[celli] = cpMean[celli];
        }

        else 
        {
            cpMean[celli] = fracSol[celli] * cpbSol.value() + (1.0 - fracSol[celli]) * cpbLiq.value();
            cpEff[celli] = cpMean[celli] + Lb.value() * ( 1.0 / TbDelta.value() ); // currently Lb gets multiplied with constant slope of fake linear function for fracSol 
                                                                                   // -> replace with real slope term later on (plus in setInitialFracSol.H (now probably not used anymore) )
        }

        lambda[celli] = fracSol[celli] * lambdabSol.value() + (1.0 - fracSol[celli]) * lambdabLiq.value();


    // calculate fracSol depending on chosen method
        if (fracSolMethod == "dataBased")
        {
//            Info << "fraction solid calculated via dataBased method" << endl;

            // only do the calculations when temperature is between TLiq and TSol
            // otherwise the values 0 and 1 can be set directly
            if (Temp[celli] < TbSol.value() )
            {
                fracSol[celli] = 1.0;
            }

            else if (Temp[celli] > TbLiq.value() )
            {
                fracSol[celli] = 0;
            }

            else
            {
                // fracSol right now uses a "fake" linear function meant to be replaced later on by a method which grabs fracSol(T) values from a table
                fracSol[celli] = (-1.0 / TbDelta.value() ) * (Temp[celli] - TbSol.value() ) + 1; 
            }
        }
Oh and if you get confused with a "b" in one of the variables' names - that just refers to the phase b in my code...
ThomasV is offline   Reply With Quote

Old   December 7, 2013, 05:28
Default
  #72
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi,
I have implemented the method that I mentioned. It's working now and giving the correct value. Now I have some things to dig out:
(1)Anyone please suggest me, how can I find out the interface velocity ? (I want to calculate the interface velocity in the solver )
(2) Have anyone worked with the non equillibrium solidification considering the under cooling ?

cheers
ahmmedshakil is offline   Reply With Quote

Old   December 9, 2013, 04:34
Default
  #73
Member
 
Thomas Vossel
Join Date: Aug 2013
Location: Germany
Posts: 45
Rep Power: 4
ThomasV is on a distinguished road
Quote:
Originally Posted by ahmmedshakil View Post
(1)Anyone please suggest me, how can I find out the interface velocity ? (I want to calculate the interface velocity in the solver )
Have a look at the fvc::interpolate(...) and fvc::reconstruct(...) commands in order to get surface values and reconstruct them into cell values. You can see them in use in for example the twoPhaseEulerFoam solver...
ThomasV is offline   Reply With Quote

Old   December 11, 2013, 07:10
Default
  #74
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi,
I'm little bit of puzzle of a term. Please help me to correct the definition:
In the energy equation:
d/dt(rho*Cp T) + div(rho U Cp T) = laplacian(K T) + ST, I am a bit of concern about the source term, whether it be...
(1) ST = rho*L*d/dt(gamma), ONLY THE LATENT HEAT PRODUCT OR
(2) ST = dH*d/dt(gamma), where dH = rho*(Cl-Cs)*(T-Tref)+rho* L. Where, Cl = liquid specific heat and Cs is the solid specific heat.
ahmmedshakil is offline   Reply With Quote

Old   December 11, 2013, 09:18
Default
  #75
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
Hello Mohammad,

you probably find it confusing but both is right in its own way. If you start form the energy conservation equation in form of enthalpy, the question is, how you define the relation between temperature and energy.
  1. For you first choice, the definition is dh = cp\, dT which implies that the derivatives in the energy equation are only applied on T, not on cp. The specific heat should therefore be outside the derivatives like cp\,\frac{d (\rho\, T)}{dt} + cp \, div(\rho\, U\, T) = laplacian(K\, T) + ST (where you put \rho is a question of conservative or non-conservative definition of the equation -- both is right)
  2. For the second choice you give, you define the enthalpy by the integral h = h_0 + \int \limits_{T_{ref}}^T cp\, dT. Then, you obtain the additional term with cps, cpl after the integration but only if you have two values, one constant one per phase.
Bird, Steward, Lightfoot has a very nice long table on that, with all possibilities and they claim the first version to be right when using the energy conservation equation in terms of temperature. By the way, there is an additional part of the source term coming from the convection term. For that, take a look at Fabian Röslers code and paper.

#####

For you question how to get the interface velocity: take a look at the definition of the VoF function. There it is \frac{d \gamma}{dt} + u_{int}\,\cdot\,\nabla \gamma = 0, where 0 is only true for a system without phase change. If the phases are liquid and solid of one material, for example, the right hand side of the equation yields the rate of phase change. u_int is than:
\vec{n}\cdot\,u_{int} = \frac{1}{|\nabla S|}\,\frac{\partial\,S}{\partial\,t}
which is the normal component of the interface. I guess, you aim onto the cooling rate, if so, you only need the normal component. Take care when coding, avoid to divide by zero outside the mushy zone where f_l, f_s (both can be inserted for S) do not change and the gradient in the denominator is zero.

By the way, you get this equation style using [ math ] ... [ / math ] without the spacings, and you simply write your equation in Latex style. Simply click on one of the equations to see that.

Have a nice day,
Anja
AnjaMiehe is offline   Reply With Quote

Old   December 17, 2013, 05:44
Default Question about Fabian Test Case (Gallium)
  #76
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi,
I have a query about the Fabian Gallium test case. I have seen that, in the fvSolution file, the solver has been chosen for p_rgh and T , but not for U. Can any one please explain the reason behind it ?
in fvSolution file:
Code:
solvers
{
    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-8;
        relTol          0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    T
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }

    TFinal
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-10;
        relTol           0;
    }
}
ahmmedshakil is offline   Reply With Quote

Old   December 17, 2013, 09:55
Default
  #77
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
Hello Mohammad,

do please take a look at posts No #31,#32 of this thread. Fabian is giving an idea about it there.

Regards, Anja
AnjaMiehe is offline   Reply With Quote

Old   December 18, 2013, 06:58
Default melting/solidification with temperature dependent thermo-properties
  #78
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi Anja Miehe,
Sorry it's me again. In my solver, I want to account for temperature dependent thermo-physical properties during melting/solidification. The governing equation that I'm solving is similar to Fabian and as follows,
cp\,\frac{d (\rho\, T)}{dt} + cp \, div(\rho\, U\, T) = laplacian(K\, T) + ST

and the source term is updated as iterative procedure as stated by Voller, like Sp and Sc. And at each iteration I'm updating the liquid fraction. Until now it works fine, and I have also validated the case with Gallium melting. The result is good. But the problem I'm facing now is to work with the temperature dependent thermo-physical properties. The code diverges after 2/3 time steps. I'm sure, the way I am working with the temperature dependent properties are not correct.
Please give me some advice to deal with the temperature dependent thermo-physical properties in the case of melting/solidification problem.
**My solver and case file is attached.

Thanks in advance
Attached Files
File Type: gz meltEnthalpyFoam.tar.gz (7.3 KB, 46 views)
File Type: gz gallium.tar.gz (2.8 KB, 31 views)
fabian_roesler likes this.
ahmmedshakil is offline   Reply With Quote

Old   December 18, 2013, 12:15
Default
  #79
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
Hello Mohammad,

no need to be sorry, I am glad when I can help. Good to see that you found a way to code the Voller scheme like you intended. Well done.

I did not try to compile your code as I am already in my Christmas Holidays and have only my little home laptop with me. Nonetheless, there are some things striking me.
  1. I do not really agree to where you placed your "include TEqn.H". When you take a look at bouyantBoussineqPimpleFoam, the TEqn is included directly after the UEqn. The advantage is that when you do several pimple iterations, the TEqn is updated as well. Would you solve without phase change and only vary cp and lambda, everything is still fine as the UEqn is not affected by the results for T. In case of phase change and when varying the viscosity as well, both equations are highly coupled and should both be within the Pimple Iteration.
  2. Line 18 in liquidFraction.H seems strange to me. Why are you adding 1 or 0 which are the only values pos() will give? This is not multiplied by anything.
  3. You should update cpS, cpL, lambdaS and lambdaL also within the while slope. They greatly affect the properties.
  4. You use alpha.storePrevIter but you do not relax the fraction of liquid. This is necessary the more complicated the chosen model for the liquid fraction and the properties are, as they are calculated using the fraction as well. In liquidFraction.H add alpha.relax() at the end of the file.
  5. If your properties do vary highly non-linear, try "Gauss harmonic uncorrected" instead of linear in the laplacienSchemes entry in fvSchemes. The second word accounts for the interpolation scheme used, and in the case described, linear would not be sensible.
  6. For the relaxation factors add the following code at the end of fvSolution.
Code:
relaxationFactors
{
  fields
  {
    "(alpha).*"     0.1;
    "(p_rgh).*"     0.3;
  }
  equations
  {
    default            0.7;
    "(T).*"            0.7;
    "(U).*"            0.7;
  }
}
Use an entry in "equations" whenever you do XEqn.relax() {TEqn.relax(), UEqn.relax()} and an entry in fields when it says X.relax() {p_rgh.relax(), alpha.relax() } in the code. Of course, you should choose the factors to your needs.

I hope, this will do it already.
Regards and happy Xmas,
Anja
apple-tree likes this.
AnjaMiehe is offline   Reply With Quote

Old   December 20, 2013, 09:44
Post solid/liquid phase change solver convMeltFoam
  #80
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
Hi all

First of all, sorry for backing out from the very nice discussion in this thread in the last few months. I was packing together my thesis. But now it is officially submitted and two months of vacation are over. I appreciate that the solid/liquid phase change topic is still strongly represented in this forum and I am happy that Anja and Mohammad shared their solvers with us.

I went through all of your posts and there is the topics of the iterative correction of the liquid fraction that should be discussed. The method used in Anjas Solver (Post 70) is iterating the energy conservation equation but is not correcting the liquid fraction.
I recommend this early Voller paper to understand the necessity of the iterative update:

Quote:
V. R. Voller: Implicit Finite-difference Solutions of the Enthalpy Formulation of Stefan Problems, IMA Journal of Numerical Analysis, Bd. 5, Nr. 2, S. 201–214, 1985.
The energy conservation equation

\rho c\frac{\partial T}{\partial t}=\nabla \cdot \left( \lambda \nabla T \right)-\rho L\frac{\partial {{\gamma }_{l}}}{\partial t} (1)

is discretized as follows:

\left[ {{a}_{p}}+{{\left( \rho cV \right)}_{p}} \right]T_{p}^{{}}={{\left( \rho cV \right)}_{p}}T_{p}^{old}+\sum\limits_{nb}{{{a}_{nb}}T_{nb}^{{}}}+{{\left( \rho cV \right)}_{p}}\left[ \gamma _{l,p}^{old}-\gamma _{l,p}^{{}} \right] (2)

As gamma is calculated from temperature, this algebraic equation is non linear. Voller linearized it using an iterative method:

\left[ {{a}_{p}}+{{\left( \rho cV \right)}_{p}} \right]T_{p}^{m+1}={{\left( \rho cV \right)}_{p}}T_{p}^{old}+\sum\limits_{nb}{{{a}_{nb}}T_{nb}^{m+1}}+{{\left( \rho cV \right)}_{p}}\left[ \gamma _{l,p}^{old}-\gamma _{l,p}^{m} \right] (3)

This equation can now be solved with any code (like OpenFOAM). However, as the source term was linearized, temperature,liquid fraction and enthalpy don not necessarily have to correlate. For this purpose, within the iterative procedure, a correction of the liquid fraction has to be performed. Have a look at the attached picture.
When enthalpy is conserved, the following is correct:

\int\limits_{T_{p}^{m}}^{T_{p}^{m+1}}{\frac{dh}{dT}}dT=\int\limits_{T_{p}^{m}}^{T_{p}^{corr}}{\frac{dh}{dT}}dT or h_{p}^{pre}=h_{p}^{m+1} (4)

Substituting the enthalpy temperature correlation

h\left( T \right)=cT+{{\gamma }_{l}}L (5)

into (4) leads to the equation to calculate the corrected liquid fraction:

cT_{p}^{m+1}+\gamma _{l,p}^{m}L=cT_{p}^{corr}+\gamma _{l,p}^{m+1}L or \gamma _{l,p}^{m+1}=\gamma _{l,p}^{m}+\frac{c}{L}\left( T_{p}^{m+1}-T_{p}^{corr} \right) (6)

The steeper the slope of the enthalpy-temperature during phase change, the more unstable gets this method. Thus, equation (6) has to be under-relaxed leading to

\gamma _{l,p}^{m+1}=\gamma _{l,p}^{m}+\omega \frac{c}{L}\left( T_{p}^{m+1}-\left. \gamma _{l}^{-1} \right|_{p}^{m} \right) (7)

with omega in the range of 0.6 to 0.8. However, one can consider the whole term in front of the brackets as under-relaxation factor which is much smaller. During iteration, the bracket will turn to zero as

T_{p}^{m}=T_{p}^{corr}=\left. \gamma _{l}^{-1} \right|_{p}^{m} (8)

So this is the iterative correction of the liquid fraction. However, this method will lead to nonsense when the energy conservation equation is not solved at least once per iteration step. In standard OpenFOAM, one can not give a minimum number of iterations in the fvSolution dictionary. So there are two possibilities to still achieve at least one iteration:

  1. Set the final residual to 1e-16. This value will never be reached and thus, OpenFOAM will always solve the conservation equation. However, this will lead to unnecessary iterations of the energy conservation equation in the first iterations of the correction method.
  2. The matrix solvers can be changed to take a minimum number of iterations into account. This is the method I used and it is already implemented in OpenFOAM-extend.

Never the less I would like to share the basic solver on which the solver from my thesis is founded as well as the necessary changes in the ldu matrix solvers to use the code in an efficient way.

Extract the data inside the convMeltFoam.tar.gz and run the install script. This will provide the changes inside the ldu matrix solvers, compile OpenFOAM src to make the changes and install the solver. However, be careful as this script will only run with OF 2.2.x.

Cheers
Attached Images
File Type: png Correction.png (16.9 KB, 96 views)
fabian_roesler is offline   Reply With Quote

Reply

Tags
melting openfoam

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
conduction problem venkataramana OpenFOAM 3 December 1, 2013 08:30
natural convection problem for a CHT problem Se-Hee CFX 2 June 10, 2007 06:29
Adiabatic and Rotating wall (Convection problem) ParodDav CFX 5 April 29, 2007 19:13
Melting Problem M FLUENT 0 April 29, 2007 16:07
Is this problem well posed? Thomas P. Abraham Main CFD Forum 5 September 8, 1999 14:52


All times are GMT -4. The time now is 15:57.