CFD with OpenSource software

A course at Chalmers University of Technology

Taught by Håkan Nilsson

Project work:

An introduction to twoPhaseEulerFoam
with addition of an heat exchange model

Developed for OpenFOAM-2.3.x

Under Review by:

Author:
Surya Kiran Peravali
Alessandro Manni
Håkan Nilsson

Disclaimer: This is a student project work, done as part of a course where OpenFOAM and some
other OpenSource software are introduced to the students. Any reader should be aware that it
might not be free of errors. Still, it might be useful for someone who would like learn some details
similar to the ones presented in the report and in the accompanying files. The material has gone
through a review process. The role of the reviewer is to go through the tutorial and make sure that
it works, that it is possible to follow, and to some extent correct the writing. The reviewer has no
responsibility for the contents.

December 11, 2014

Introduction

twoPhaseEulerFoam is a ”solver for a system of 2 compressible fluid phases with one phase dispersed
(e.g. gas bubbles in a liquid) including heat-transfer” since the last version (the capability of manage
compressible fluids has been added in the 2.1.x version with compressibleTwoPhaseEulerFoam).
Assuming that different phases could be mathematically descripted as a continuum media(we are
excluding the lagrangian approach), we could use two different ways to derive the equation of the
multiphase flux: the averaging or the mixture approach. In the first one , known as Euler Euler,
the equations are considered and solved one by one for each phase. Conversely,in the mixture
one,well known as Volume Of Fluid (VOF,used in the family of solvers interFoam), is used a single
equation from which the results are taken from the costitutive equation of the considered phase.
To get an idea, both the approaches leads to a set of equations of the same kind(state equations,
interactions between phases,momentum mass and energy balances). In our case the Eu-Eu approach
is used,considering a solid phase with fluid like behaviour without particle tracking (typical of the
lagrangian framework) reducing the numerical effort required with a large number of particles. As a
consequence,will be necessary introduce a physical definition for the particle stress tensor to provide
a mathematical closure of the problem.

Figure 1: Example of Fluidized Bed

1
Chapter 1

Governing Equations

1.1 Mass Momentum and Energy Eqns

The Eu-Eu solvers, are based on the assumpion of treating both phases as interpenetrating continua.
In this case with two phases, the following assumption must be verified:

αs + αg = α1 + α2 = 1 (1.1)

Being α the phase volume fraction in the cell and the subscripts s and g related for example to a solid
and a gas phase.Hereafter we willl use the subscripts 1 and 2 to be in line with the notation used in
the last release of the code. The mass continuity equation for both the phases are(not considering
fvOptions capability added only recently in 23x that will be shown later):

∂(α1 ρ1 )
+ ∇ · (α1 ρ1 ~u1 ) = R12 = 0 (1.2)
∂t

∂(α2 ρ2 )
+ ∇ · (α2 ρ2 ~u2 ) = R21 = 0 (1.3)
∂t
The momentum equations for both the phases are:

∂(α1 ρ1 ~u1 )
+ ∇ · (α1 ρ1 ~u1 ~u1 ) = ∇ · S 1 + (α1 ρ1~g ) + I~1 (1.4)
∂t
∂(α2 ρ2 ~u2 )
+ ∇ · (α2 ρ2 ~u2 ~u2 ) = ∇ · S 2 + (α2 ρ2~g ) − I~2 (1.5)
∂t
where
2
∇ · S 1 = ∇ · (α1 τ 1 ) − α1 ∇p − ∇p1 with τ 1 = µ1 D1 + (λ1 − µ1 )(∇ · ~u1 )I (1.6)
3
is the solid stress tensor div in which ∇p1 is gradient of the granular pressure (see section 1.2.1) and
λ1 is the solid bulk viscosity and
1
D1 = [∇~u1 + (∇~u1 )T ] (1.7)
2
2
∇ · S 2 = ∇ · (α2 τ 2 ) − α2 ∇p with τ 2 = µ2 D2 − µ2 (∇ · ~u2 )I (1.8)
3
1.8 is the gas stress tensor div

I~1,2 = drag + vme + lif t + walllubri + turbulentdispersion (1.9)

is the sum of the possible momentum interactions.For example, the drag force, proportional to
the relative velocity of the phases and to a coefficient that could be evaluated with well known
correlations (Gidaspow,Ergun,Syamal and O’Brien, Wen-Yu, Hill and Koch etc...):

2
1.1. MASS MOMENTUM AND ENERGY EQNS CHAPTER 1. GOVERNING EQUATIONS

K12 (~u1 − ~u2 ) (1.10)

 ρ2 α2 α1 |~
u2 −~
u1 | −2,65
3

 4 CD d1 g α2 ≥ 0.8 Gidaspow
K12 =
150s (1−α2 )µ2 1.75ρ2 α1 |~
u2 −~
u1 |
+ α2 < 0.8 Ergun


α2 d21 d1

where: 
 24/Re(1 + 0.15Re0.687 ) per Re < 1000
CD =
0.44 per Re ≥ 1000


with:
ρ2 α2 |~u2 − ~u1 |d1
Re = (1.11)
µ2
The energy equation, on which we will ace in the last part of the report ,could be written in the
following two formulations depending on the user preferences:
∂[αi ρi (hi + ki )] ∂p
+ ∇ · [αi ρi (hi + ki )~ui ] = αi + ∇ · αi αef f ∇hi + Kht ∆T (1.12)
∂t ∂t
and, remembering that h=e+pv
∂[αi ρi (ei + ki )] ∂α1
+ ∇ · [αi ρi (ei + ki )~ui ] = −[ p + ∇ · α1 ~u1 p] + ∇ · αi αef f ∇hi + Kht ∆T (1.13)
∂t ∂t
where ki is the kinetic energy and αef f is the effective thermal diffusivity. In the following, we can
see the implementation of the equation in his fully conservative formulation (rearranged to be read
easier,here is shown the he1eqn):
fvScalarMatrix he1Eqn
(
fvm::ddt(alpha1, rho1, he1) + fvm::div(alphaRhoPhi1, he1)
+ fvc::ddt(alpha1, rho1, K1) + fvc::div(alphaRhoPhi1, K1)

// Compressibility correction
- fvm::Sp(fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1), he1)
- (fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1))*K1

==

+ (
he1.name() == thermo1.phasePropertyName("e")
? - fvc::ddt(alpha1)*p + fvc::div(alphaPhi1, p)
: + alpha1*dpdt
)

fvm::laplacian
(
fvc::interpolate(alpha1)
*fvc::interpolate(thermo1.alphaEff(phase1.turbulence().mut())),
+ he1
+ )

heatTransferCoeff*(thermo2.T() - thermo1.T())
+ heatTransferCoeff*he1/Cpv1
- fvm::Sp(heatTransferCoeff/Cpv1, he1)
);

3
1.2. KINETIC TH. AND FRICTIONAL MODELS CHAPTER 1. GOVERNING EQUATIONS

1.2 Kinetic Th. and Frictional Models

As said before, we need a mathematical model to close the global solid stress tensor. For this
purpose, different flow regimes should be taken into account:
• Kinetic Regime → particles present in a ”diluite” state
• Collisional Regime → particles present in a higher conc. with collisions

• Frictional Regime → particles present in a very high conc.with Coulombian friction

For the first twepo regimes , the Kinetic Theory of Granular Flows (by Gidaspow and Lun)gives a
closure introducing the concept of granular temperature, defined as the kinetic energy due to the free
motion of the particles. This granular temperature is then solved (with the aid of other submodels)
in a balance equation. For the frictional regime, as will be shown later, are used correlations taken
from the mechanics of soils.To better understand the theory, wi will pass through the definition of
the solid stress tensor , to come across the granular energy balance. You can find the implementation
of that theory in:

1.2.1 Solid stress Tensor-Kinetic Regime

Exists a lot of attempts to combine both the Kinetic/Collisional and Frictional regimes but is
common practice to use a value of switch α2∗ , in which we can assume the change of regime (from
here starts the problem of maximum packing αminF riction to αmax ):
p
 −pk1 I + τ 1 + p if α2 ≤ α2∗


S1 = (1.14)
f f ∗
−p1 I + τ 1 + p if α2 > α2


where p1 is the solid pressure (hydrostatic) and τsm is the deviatoric part of the stress tensor,that is
the viscous one.The apexes k, f are referred respectively to the kinetic and frictional regime. In the
kinetic regime, the kinetic energy of the medium flux, degrades in the energy associated with the
random fluctuation of the particles , and then dissipated into heat because of inelastic collisions. This
process is very similar to that of the turbulence in the monophase flows . The kinetic energy of the
fluctuations,is taken into account for this theory by a granular temperature , θ1 , which is obviously
different from the temperature of the solid phase , measuring the vibrational kinetic energy of that
phase.
In our case the terms of the viscous stress are based on a theory reported by BGM van Wachem
in his doctoral thesis developed by Lun in 1984 resulting from the kinetic theory of smooth spherical
particles , subject to inelastic collisions with no exchange of momentum ( due to inability to obtain
in this case a closed solution for granular temperature when the void fraction become preponder-
ant,about 1 ) .The resulting expressions for stresses are outlined below . The granular pressure is
given by:
pk1 = Ks1 α12 θ1 (1.15)
where
Ks1 = 2(1 + e1 )ρ1 g0 (1.16)
The granular stresses are given by:

τ k1 = 2µk1 D1 + λk1 tr(D1 )I (1.17)

where λk1 √
λk1 = Ks2 α1 θ1 (1.18)

4
1.2. KINETIC TH. AND FRICTIONAL MODELS CHAPTER 1. GOVERNING EQUATIONS

and where the constant Ks2 is given by:

4d1 ρ1 (1 + e1 )α1 g0 2
Ks2 = √ − Ks3 (1.19)
3 π 3

and the constant Ks3 is

√
4d1 ρ1 π 8α1 g0 (1 + e1 )
Ks3 = { [1 + 0.4(1 + e1 )(3e1 − 1)α1 g0 ] + √ } (1.20)
2 3(3 − e1 ) 5 π

in which the solid dynamic viscosity µk1 ,is given by:

√
µk1 = Ks3 α1 θ1 (1.21)

1.2.2 Solid stress Tensor-Frictional Regime α ' 0.6

The stress in the frictional regime , are usually described by adopting theories arising from soil
mechanics . The stresses , caused friction , are described by phenomenological models rather than
mechanistic models as happens in the case kinetic . the soil mechanics uses a set of relationships
between the components of the stress tensor and the strain tensor . this arbitrary function that
allows a certain compressibility of the solid phase , allows to define the pressure of the solid phase
in the frictional regime as:
pf1 = α1 p∗ (1.22)
where p∗ is calculated by a law of the type:

p∗ = A(α2∗ − α2 )n (1.23)

where A typically assumes values close to 1025 and n = 10. The deviatoric part of the stress tensor
is calculated using a formula proposed by Schaeffer (1987):

τ f1 = 2µf1 D1 (1.24)

where
p∗ sin φ
µf1 = √ (1.25)
2 I 2D
and I2D i the second invariant of the steess tensor (the subscript ”s” means 1,for the solid phase):
1
I2D = [(Ds11 − Ds22 )2 + (Ds22 − Ds33 )2 + (Ds33 − Ds11 )2 ] + Ds12
2 2
+ Ds23 2
+ Ds31 (1.26)
6
Typically the values of viscosity in the frictional regime assume very high values and techniques
agents on the frictional angle are sometimes used to avoid nonphysical behavior. In this brief
description it is limited to a discussion of the problem needed to understand the main idea. For a
complete discussion , please refer to specialized texts.

1.2.3 Granular energy equation

As seen before , the description of a flow of smooth spherical particles that collide inelastically , is
used for the derivation of the equations of the kinetic theory associated with them , producing a
stress tensor solid underlineS1 . The resulting constitutive equations contain the quantity θ1 called
granular temperature of the solid phase . The granular temperature is then proportional to the
energy of the continuous granular , granular where the energy is defined as the kinetic energy of the
specific component of the floating speed of the particles :
3 1
θ1 = C12 (1.27)
2 2

5
1.3. NUMERICAL TREATMENT CHAPTER 1. GOVERNING EQUATIONS

whereC1 is the floating component of the istantaneus velocity cm defined as:

~1
~c1 = ~u1 + C (1.28)

then to calculate the value of such energy granular will be necessary, as said before , define a transport
equation of the granular temperature for the solid phase 1 :
3 ∂α1 θ1
ρ1 [ + ∇ · (α1 θ1 ~u1 )] = S 1 : ∇~u1 + ∇(kθ1 ∇θ1 ) − γθ1 + E12 (1.29)
2 ∂t
where the first term on the second member is the production of fluctuation due to the deviatoric
stress , the second a diffusive term of fluctuating energy , the third is due to dissipation of inelastic
collisions while the last term represents an exchange between the phases. With regard to the first
term is necessary to define the constitutive equations to close the governing equations :
2
S m = (−p1 (k, f ) + ηλ1 ∇~u1 ) + µ1 [(∇~u1 + (∇~u1 )T ) − ∇~u1 ] (1.30)
3
for the correlations of the various terms used in the above equation, such as the bulk viscosity λ1 ,
that of the solid phase µ1 , conductivity kθ1 correlations to the terms of exchange and collisional ,
refer to the correlations reported in the subfolders viscosityM odel and conductivityM odel of the
folder kineticT heoryM odels .To finish this brief introduction to the physics of the problem ,is shown
the algebraic solution of the granular temperature derived from the equality of terms of production
and energy dissipation granular same :
p
 −Ks1 α1 tr(D1 ) 2 tr(D )2 α2 + 4K α [K tr(D )2 + 2K (D D T )]
K1m 1 1 s4 1 s2 1 s3 1 1
θ1 = + (1.31)
2α1 Ks4 2α1 Ks4
where the terms Ks1 , Ks2 , Ks3 assumes the same meaning seen before and the term Ks4 is egual to:

12(1 − e21 )ρ1 g0

Ks4 = √ (1.32)
d1 π

1.3 Numerical treatment

An issue of the resolution of the mass and momentum equations, is that generally the α is present
everywhere, having a null equation when we are in a monophase zone(only fluid phase). Henry
Weller, in the first releases of the solver, solved that problem with a non conservative form of the
momentum equation (Oliveira and Issa, 2003) expanding the derivate where α1 is present,to come out
the mass conservation equation (=0!) and dividing all for α1 . In this way , the solution could be not
so accurate with steep gradients of the solid phase.In fact the non-conservative form of the momentum
equation presents difficulties when strong property transfer (Park et al., 2009) or strong particle
pressure gradients are present. In the last release (23x 7 months ago) has been implemented the fully
conservative formulation of the energy and momentum equation solving the problem and increasing
the complessive robustness(Passalacqua, Fox 2011). This development required the additions to
fvmDdt/fvcDdt to support time derivatives of triples, as you can see in the U Eqns.H. Should be
noticed that the convergence of the energy equations is more difficult when using the conservative
form especially with cases involving MRF. Usually however, convergence is achieved iterating over
the system of equations using a sufficiently small time-step. As suggested in the latest tutorial in
2.3.1 just released ,is recommended to set the number of outer correctors up to 3 (also 20), and
use a conservative value of the maximum Courant number (0.2 - 0.4 surely ¡0.5). It is also a good
practice to setting the residualAlpha to 1.0e-4 everywhere and checking it by add residualControl
to fvSolution/PIMPLE.
We can now start to look at the code,in particular the file twoPhaseUlerFoam.C that is the main:

cd $FOAM_SOLVERS/multiphase/twoPhaseEulerFoam/
vi twoPhaseEulerFoam.C

6
1.3. NUMERICAL TREATMENT CHAPTER 1. GOVERNING EQUATIONS

// --- Pressure-velocity PIMPLE corrector loop

while (pimple.loop())
{
fluid.solve(); //starts solving the continuity equation
fluid.correct();
volScalarField contErr1
(
fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
- (fvOptions(alpha1, rho1)&rho1)
);
volScalarField contErr2
(
fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
- (fvOptions(alpha2, rho2)&rho2)
);
#include "UEqns.H" //include the momentum equation,not solved
#include "EEqns.H" //solves the energy equation
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H" //construct a combined continuity eq.
and correct the velocities
}
#include "DDtU.H"
if (pimple.turbCorr())
{
fluid.correctTurbulence(); //calculate the turbulent viscosity
}
}

we can go now through the code, to understand all the passages included in the just seen piece
of code

1.3.1 Volume fraction Solution

the core of the solution is the MULES ”Multidimensional Universal Limiter with Explicit Solution”
that is an iterative implementation of the Flux Corrected Transport technique (FCT) , used to
guarantee boundedness in the solution of hyperbolic problems. MULES compute a corrected flux
between an high and a low order scheme solution with a weighting factor λ
FC = FL + λA where A is the anti − dif f usive f lux = (FH − FL ) (1.33)

To determine the weighting factor, Zalesak (1979) underlined that the value of the net flux in
a cell , must be neither greater than a local maximum nor lesser than a local minimum, correcting
the flux in order to not create new maximum or minimum. A possible question for a careful reader
could be the necessity of it’s widely usage with multiphase flow. Looking at the eqn 1.1, we find the
solution. In multiphse flows,we need to set global extrema for the problem In our case α1 + α2 = 1
and other limits imposed by the users.
Parameters are the variable to solve, the normal convective flux and the actual explicit flux of
the variable which is also used to return limited flux used in the bounded-solution. In the following
you can see pieces of the implementation MULES , starting from the templated class in which are
implemented the methods for the explicit solution (first piece of the code) and the limiting capability.

src\finiteVolume\fvMatrices\solvers\MULES MULESTemplates.C

EXPLICIT SOLUTION
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,

7
1.3. NUMERICAL TREATMENT CHAPTER 1. GOVERNING EQUATIONS

const SpType& Sp,

const SuType& Su,
const scalar psiMax,
const scalar psiMin
)
{
const fvMesh& mesh = psi.mesh();
const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
psi.correctBoundaryConditions();
limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
LIMIT
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::limit
.....
surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;
.....
limiter
(
allLambda,
rDeltaT,
rho,
psi,
.
.
.
);
if (returnCorr)
{
phiCorr *= lambda;
}
else
{
phiPsi = phiBD + lambda*phiCorr;
}
}

in the same class is then reported how the explicit solution is really evaluated:
void Foam::MULES::explicitSolve
( const RdeltaTType& rDeltaT,const RhoType& rho,.......)
{ Info<< "MULES: Solving for " << psi.name() << endl;
const fvMesh& mesh = psi.mesh();
scalarField& psiIf = psi;
const scalarField& psi0 = psi.oldTime();
psiIf = 0.0;
fvc::surfaceIntegrate(psiIf, phiPsi);
if (mesh.moving())
{ }
else
{
psiIf =
(
rho.oldTime().field()*psi0*rDeltaT
+ Su.field()
- psiIf
)/(rho.field()*rDeltaT - Sp.field());
}
psi.correctBoundaryConditions();
}

as you can see, in the else brackets is reported the formulation of the explicit solution:

∂(α1 ρ1 ) αn − α10
+ ∇ · (α1 ρ1 ~u1 ) = 1 + psiIf =
∂t δt

8
1.3. NUMERICAL TREATMENT CHAPTER 1. GOVERNING EQUATIONS

α01
+ Su + psiIf
= Sp α1n + Su → α1n = δt
1 (1.34)
δt − Sp

where psiIf is the net flux,while Sp α1n + Su is the source term (if you want i can add here an
explanation of the terms Sp and Su and their relationship with the dgdt that you can see in the
code (for example in the pEqn.H)).Must be noticed that in the MULES algorithm the density is
considered constant.
Before executing MULES, for the mass conservation,a trick is used for the solid fraction solu-
tion(Rusche PHD thesis). The velocities aren’t considered as ~u1 and ~u2 but are defined the mixture
and relative velocities ( with the respective fluxes)so that in this way, a laplacian term appear (div
grad in 1.38) This approach is very useful in the solution of the pressure equation.

u = α1 u1 + α2 u2 ur = u1 − u2 → u1 = u + ur (1 − α1 ) (1.35)
∇p1 ∇p1 ∇p1 α1 ∂p1
u = (u + )− = u∗ − → u∗ = u + ∇α1 (1.36)
λs λs λs λs ∂α1
α1 ∂p1 1
u∗ = u + ∇α1 with λs = (1.37)
λs ∂α1 As + Kρ12
s

where As condensates the diagonal coefficients of the matrix of the momentum equation (Passalac-
qua,Fox 2011).
Then the solid phase continuity equation is finally solved:
∂α1 α1 ∂p1
+ ∇ · (α1 ~u∗1 ) − ∇ · [ ∇α1 ] = Sp α1 + Su (1.38)
∂t λs ∂α1

1.3.2 Momentum Equation

In this file, the momentum equations are constructed but still not solved.It’s important to notice
that are not included the terms related to the pressure and the solid pressure .
UEqns.H
U1Eqn =
(
fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
- fvm::Sp(contErr1, U1)
+ mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
+ phase1.turbulence().divDevRhoReff(U1)
==
- liftForce
- wallLubricationForce
- turbulentDispersionForce
- virtualMassCoeff
*(
fvm::ddt(U1)
+ fvm::div(phi1, U1)
- fvm::Sp(fvc::div(phi1), U1)
- DDtU2
)
+ fvOptions(alpha1, rho1, U1)
);
U2Eqn =
..........................

1.3.3 Pressure equation

Summing the continuity eqns 1.39,1.40, we obtain the combined continuity equation 21 :
Di ∂
( where Dt = ∂t + ui · ∇)
∂α1 α1 D1 ρ1
+ ∇ · (α1 ~u1 ) = − (1.39)
∂t ρ1 Dt

9
1.3. NUMERICAL TREATMENT CHAPTER 1. GOVERNING EQUATIONS

∂α2 α2 D2 ρ2
+ ∇ · (α2 ~u2 ) = − (1.40)
∂t ρ2 Dt
α1 D1 ρ1 α2 D2 ρ2
∇ · ~u = − − (1.41)
ρ1 Dt ρ2 Dt
where the sum of the terms ∂α ∂α2
∂t and ∂t is ugual to zero for obvious reasons. As well known, the
1

momentum equation could be written in short in the semi discrete form:

As Usp = Hs (~u1 ) − α1 ∇p (1.42)

Ag Ugp = Hg (~u2 ) − α2 ∇p (1.43)

Combining togheter the equations 1.41,42 and 43 we obtain the pressure equation :
Hs Hg α1 α2 α1 D1 ρ1 α2 D2 ρ2
∇ · ~u = ∇ · (α1 + α2 ) − ∇ · [(α1 + α2 )∇p] = − − (1.44)
As Ag As Ag ρ1 Dt ρ2 Dt

that is solved in pEqn.H . Is useful to remark that Hi /Ai is the velocity predictor.

1.3.4 Turbulence modeling

The issue of compressibility has been managed for many years using two distinct turbulence mod-
elling frameworks depending on the constant/variable density flows. In the case of multiphase
turbulence the phase-fraction must be included into all transport and source terms of the turbulence
property equations making necessary the duplication of the code. Therefore, only for twoPhaseEuler-
Foam in this version, a new single, general templated turbulence modelling framework was created
that covers all four of the combinations(compr/uncompr/monoph/multiph). The framework can
be found in F OAM SRC/T urbulenceM odel and includes the incompressible, compressible, phase-
incompressible and phase-compressible options including support for laminar, RAS and LES simu-
lations with wall-functions and other specialised boundary conditions. This new framework is very
powerful and supports all of the turbulence modelling requirements needed so far.

1.3.5 Solution procedure: the PIMPLE algorithm

It’s a combination of PISO and SIMPLE with the power of quick solution for unsteady flows coming
from the PISO , and the possibility to use under− relaxation (look in the equation files):
• Start the mass conservation loop

– Solve the eqn 1.38 (without the particle pressure contribution) using MULES
∂p1
– Calculate ∂α s
, correct the continuity equation for a new α1 .For the p1 it’s used the
formulation related to value of α1 (k or f regime as seen in the sections 1.2.1 and 1.2.2)
– Calculates the α2
– Iterations

• Calculate and update the coefficients needed in the term I~ of the momentum equation ( drag
vme lift ecc....)
• Solve a predicted velocity field with the last calculated pressure
• Evaluates the energy equation

• Pressure equation loop

• Iterations until is reached the convergency criteria

10
Chapter 2

How to add an heat exchange

model

In the energy equations, we have seen the term Kht ∆T . Add an heat exchange model means to
add another way to calculate Kht depending on the conditions in the reactor(for instance the value
of the Re). In this last version,is present only the Ranz Marshall correlation. In this quick guide
we will add togheter the Gunn the Tomiyama and LiMason(2000)models that you can find in the
attached files. The just mentioned coefficient is
k2 N u1
Kht = (2.1)
d1
where k2 is the thermal fluid condictivity, the Nusselt number the ratio of convective/conductive
heat transfer,and d1 is the solid diameter. The Gunn correlation is more suitable for granular flows
and very high Re105 . The LiMason includes the Ranz Marshall for Re < 200 and presents two other
correlations in the range 200/1500 and over 1500,while the Tomiyama is similar to the RM.

• We can start start copying the following lines on the terminal

OF23x
cd $WM_PROJECT_DIR
cp -r --parents applications/solvers/multiphase/twoPhaseEulerFoam \
$WM_PROJECT_USER_DIR
cd $WM_PROJECT_USER_DIR/applications/solvers/multiphase
mv twoPhaseEulerFoam myTwoPhaseEulerFoam
cd myTwoPhaseEulerFoam
./Allwclean

• Now you should have a clean src code in your project user dir

• Copy the AMpresentation.tar.gz into the /heatTransferModels folder , then

cd interfacialModels/heatTransferModels/
tar xzf AMpresentation.tar.gz
ls

• Now we have to simply modify the make/files

cd ..
cd Make
vi files

• And insert after the line 31 the following lines

heatTransferModels/Gunn/Gunn.C
heatTransferModels/LiMason/LiMason.C
heatTransferModels/Tomiyama/Tomiyama.C

11
 CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

• before to compile the code we have to modify other few things:

cd ..
cd ..
mv twoPhaseEulerFoam.C myTwoPhaseEulerFoam.C
cd Make
sed -i s/twoPhaseEulerFoam/myTwoPhaseEulerFoam/g files
sed -i s/FOAM_APPBIN/FOAM_USER_APPBIN/g files
cd ..
cd interfacialModels/Make
sed -i s/FOAM_LIBBIN/FOAM_USER_LIBBIN/g files
cd ..
cd phaseCompressibleTurbulenceModels/Make
sed -i s/FOAM_LIBBIN/FOAM_USER_LIBBIN/g files
cd ..
cd twoPhaseSystem/Make
sed -i s/FOAM_LIBBIN/FOAM_USER_LIBBIN/g files

• Now we can compile!

cd ..
./Allwmake

• now we need a fresh tutorial to test

run
cp -r $FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed .
cd fluidisedBed/constant
sed -i s/RanzMarshall/ASD/g phaseProperties
cd ..
blockMesh
myTwoPhaseEulerFoam

• now check the error!

--> FOAM FATAL ERROR:
Unknown heatTransferModelType type ASD
Valid heatTransferModel types are :
4
(
Gunn
LiMason
RanzMarshall
Tomiyama
)

Another way to add submodels like these, is to compile only the libraries(wmake libso ,to build
a ”dynamic library”), without compiling the whole code, and then adding the following lines in the
control dict, related to the just compiled libraries:
libs ("myGunn.so")
libs ("myLiMason.so")
libs ("myTomiyama.so")

having obviously changed the name of the files before compile.

12
2.1. MODELS ADDED CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

2.1 Models added

Here are simply listed the .C and the .H of the just added models.
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Gunn.H"
#include "phasePair.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace heatTransferModels
{
defineTypeNameAndDebug(Gunn, 0);
addToRunTimeSelectionTable(heatTransferModel, Gunn, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::heatTransferModels::Gunn::Gunn
(
const dictionary& dict,
const phasePair& pair
)
:
heatTransferModel(dict, pair)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::heatTransferModels::Gunn::~Gunn()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::heatTransferModels::Gunn::K() const
{
volScalarField Nu((scalar(7.0) -10.0*pair_.phase1() +
5.0*sqr(pair_.phase1()))*(scalar(1.0)
+ 0.7*pow(pair_.Re(), 0.2)*pow(pair_.Pr(), 1.0/3.0))
+ (scalar(1.33) - 2.4*pair_.phase1()
+ 1.2*sqr(pair_.phase1()))*pow(pair_.Re(), 0.7)*pow(pair_.Pr(), 1.0/3.0));
return
6.0
*max(pair_.dispersed(), residualAlpha_)
*pair_.continuous().kappa()
*Nu
/sqr(pair_.dispersed().d());
}
// ************************************************************************* //

13
2.1. MODELS ADDED CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::heatTransferModels::Gunn
Description
SourceFiles
Gunn.C
\*---------------------------------------------------------------------------*/
#ifndef Gunn_H
#define Gunn_H
#include "heatTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class phasePair;
namespace heatTransferModels
{
/*---------------------------------------------------------------------------*\
Class Gunn Declaration
\*---------------------------------------------------------------------------*/
class
: Gunn
public heatTransferModel
{
public:
//- Runtime type information
TypeName("Gunn");
// Constructors
//- Construct from components
Gunn
(
const dictionary& dict,
const phasePair& pair
);
//- Destructor
virtual ~Gunn();
// Member Functions
//- The heat transfer function K used in the enthalpy equation
tmp<volScalarField> K() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace heatTransferModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

14
2.1. MODELS ADDED CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "LiMason.H"
#include "phasePair.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace heatTransferModels
{
defineTypeNameAndDebug(LiMason, 0);
addToRunTimeSelectionTable(heatTransferModel, LiMason, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::heatTransferModels::LiMason::LiMason
(
const dictionary& dict,
const phasePair& pair
)
:
heatTransferModel(dict, pair)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::heatTransferModels::LiMason::~LiMason()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::heatTransferModels::LiMason::K() const
{
volScalarField Nu = neg(pair_.Re() - 200.0)*(scalar(2.0) +
0.6*pow(pair_.phase1(), 3.5)*sqrt(pair_.Re())*pow(pair_.Pr(), 1.0/3.0))
+ pos(pair_.Re() - 200.0)*neg(pair_.Re() - 1500.0)*(scalar(2.0)
+ pow(pair_.phase1(), 3.5)*pow(pair_.Pr(), 1.0/3.0)*(0.5*sqrt(pair_.Re())
+ 0.02*pow(pair_.Re(), 0.8))) + pos(pair_.Re() - 1500.0)*(scalar(2.0)
+ 0.000045*pow(pair_.phase1(), 3.5)*pow(pair_.Re(), 1.8));
return
6.0
*max(pair_.dispersed(), residualAlpha_)
*pair_.continuous().kappa()
*Nu
/sqr(pair_.dispersed().d());
}
// ************************************************************************* //

15
2.1. MODELS ADDED CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::heatTransferModels::LiMason
Description
SourceFiles
LiMason.C
\*---------------------------------------------------------------------------*/
#ifndef LiMason_H
#define LiMason_H
#include "heatTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class phasePair;
namespace heatTransferModels
{
/*---------------------------------------------------------------------------*\
Class LiMason Declaration
\*---------------------------------------------------------------------------*/
class
: LiMason
public heatTransferModel
{
public:
//- Runtime type information
TypeName("LiMason");
// Constructors
//- Construct from components
LiMason
(
const dictionary& dict,
const phasePair& pair
);
//- Destructor
virtual ~LiMason();
// Member Functions
//- The heat transfer function K used in the enthalpy equation
tmp<volScalarField> K() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace heatTransferModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

16
2.1. MODELS ADDED CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Tomiyama.H"
#include "phasePair.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace heatTransferModels
{
defineTypeNameAndDebug(Tomiyama, 0);
addToRunTimeSelectionTable(heatTransferModel, Tomiyama, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::heatTransferModels::Tomiyama::Tomiyama
(
const dictionary& dict,
const phasePair& pair
)
:
heatTransferModel(dict, pair)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::heatTransferModels::Tomiyama::~Tomiyama()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
Foam::heatTransferModels::Tomiyama::K() const
{
volScalarField Nu(scalar(2) + 0.15*pow(pair_.Re(),0.8)*pow(pair_.Pr(),0.5));
return
6.0
*max(pair_.dispersed(), residualAlpha_)
*pair_.continuous().kappa()
*Nu
/sqr(pair_.dispersed().d());
}
// ************************************************************************* //

17
2.1. MODELS ADDED CHAPTER 2. HOW TO ADD AN HEAT EXCHANGE MODEL

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::heatTransferModels::Tomiyama
Description
SourceFiles
Tomiyama.C
\*---------------------------------------------------------------------------*/
#ifndef Tomiyama_H
#define Tomiyama_H
#include "heatTransferModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class phasePair;
namespace heatTransferModels
{
/*---------------------------------------------------------------------------*\
Class Tomiyama Declaration
\*---------------------------------------------------------------------------*/
class Tomiyama
:
public heatTransferModel
{
public:
//- Runtime type information
TypeName("Tomiyama");
// Constructors
//- Construct from components
Tomiyama
(
const dictionary& dict,
const phasePair& pair
);
//- Destructor
virtual ~Tomiyama();
// Member Functions
//- The heat transfer function K used in the enthalpy equation
tmp<volScalarField> K() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace heatTransferModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

18