CFD Simulation of Multicomponent Gas Flow Through Porous Media

CFD Simulation of Multicomponent gas flow through porous media

Master Thesis by Chethan Mohan Kumar

at Lehrstuhl fu omungsmechanik ¨ r Str¨ Bergische Universit¨ at Wuppertal

Advisors: Prof. Dr.-Ing. habil. Uwe Janoske Dipl.-Math. Markus Bu ¨ rger

11.03.2013

Declaration Hereby I declare that I have worked out the theme CFD simulation of Multicomponent gas flow through porous media entirely on my own without using foreign help, that I have only used the indicated sources and utilities and that all citations, having been taken literally or by content, are designated.

(Date)

(Signature)

Abstract The objective of this thesis is to develop a generic CFD solver to simulate multicomponent gas transport involving multiscales. A comparative study between the mixture and Eulerian approach for multicomponent flows is done. Eulerian approach where every component has its own characteristic velocity is used in this work with reasons stated. Inter-component momentum exchange term has been modelled using Maxwell-Stefan relations. A semi-heuristic drag term for modelling porous drag is used. Temperature transport for the mixture is developed by considering an ensemble averaging method of all components for the mixture. Volume averaged form for all the equations is applied by considering fluid in presence of porous media as a pseudo-homogeneous medium. All the equations and terms mentioned are implemented in the solver. An open source CFD software OpenFOAM has been used to develop the solver. Capability of the solver to simulate diffusion dominated mass transfer has been established by using a Loschmidt tube which involves diffusion of a ternary mixture. Accuracy of 0.5% is observed for the case in comparison with analytical solution for the problem in one dimension. Validation of porous drag and energy transport has been done by using a fully developed laminar flow over two parallel plates and comparing the Nusselt numbers. The Nusselt numbers of 8.85 and 9.13 for porous and non-porous zones were observed compared to the literature values of 7.54 and 9.8. The reason for the deviations are stated. Keywords: multicomponent, multiphase, porous media, Maxwell-Stefan, volume averaging, multiphaseEulerFoam, Eulerian, MULES, mixture, Nusselt, ternary, diffusion, mass transfer, multicomponentPorousFoam, Darcy, Knudsen, incompressible, Loschmidt tube

TO MY PARENTS

Acknowledgements I would like to express my deepest gratitude to Professor Uwe Janoske for his constant encouragement and guidance throughout this work. His concern and enthusiasm towards the progress of this thesis has been of great help. I take this opportunity to thank Markus B¨ urger for his suggestions and guidance without which this work would not have been the same. I must mention that his patience over my supposedly two minute discussions which often extended to prolonged conversations has been remarkable. His knowledge and approach to problem solving has influenced me to a great extent. I extend my gratitude to Shilpa for her constructive criticism over the semantics and writing style of this thesis and my family and friends for their encouragement and support. Finally, I wish to thank my parents for their unprecedented support throughout my study.

CONTENTS

Contents 1 Motivation

1

2 Multi-component gas flow

3

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

2.2

Approaches to Multicomponent modelling . . . . . . . . . . . . . . .

5

2.2.1

Volume Of Fluid approach . . . . . . . . . . . . . . . . . . . .

5

2.2.2

Mixture approach . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2.3

Eulerian approach

2.2.4

Model comparisons . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3

. . . . . . . . . . . . . . . . . . . . . . . . 10

Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Mathematical Modelling 3.1

14

Drag modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.1

Effects of Diffusion . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1.2

Maxwell-Stefan relations for modelling drag . . . . . . . . . . 16

3.2

Porosity Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.3

Energy transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.4

Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Numerical Simulation

28

4.1

Open Source CFD: OpenFOAM . . . . . . . . . . . . . . . . . . . . . 29

4.2

Solver: multiphaseEulerFoam

4.3

. . . . . . . . . . . . . . . . . . . . 29

4.2.1

Solution Method: Phase fraction . . . . . . . . . . . . . . . . 31

4.2.2

Solution Method: Pressure-Velocity coupling . . . . . . . . . . 36

4.2.3

Semi-implicit treatment of Drag term . . . . . . . . . . . . . . 37

Modified solver: multicomponentPorousFoam . . . . . . . . . . . 38

i

CONTENTS

4.4

4.3.1

Source term: Maxwell-Stefan drag . . . . . . . . . . . . . . . . 39

4.3.2

Source term: Porous drag . . . . . . . . . . . . . . . . . . . . 40

4.3.3

Energy transport . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.3.4

Solution algorithm . . . . . . . . . . . . . . . . . . . . . . . . 42

Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5 Validation 5.1

5.2

5.3

45

Case1: Loschmidt tube . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.1.1

Case Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.1.2

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

Case2: Flow through a channel . . . . . . . . . . . . . . . . . . . . . 49 5.2.1

Case Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5.2.2

Validation of temperature transport . . . . . . . . . . . . . . . 51

5.2.3

Validation of porosity effects . . . . . . . . . . . . . . . . . . . 51

Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

6 Conclusions

55

A Loschmidt tube

57

A.1 Scilab Program to calculate 1-D analytical solution . . . . . . . . . . 57 A.2 Animation of initial diffusion of gases . . . . . . . . . . . . . . . . . . 60 B Scilab program to calculate Nusselt number

63

Bibliography

64

ii

LIST OF FIGURES

List of Figures 2.1

Examples of multiphase flows . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Multicomponent flow in porous media . . . . . . . . . . . . . . . . . .

4

2.3

Volume of Fluid approach . . . . . . . . . . . . . . . . . . . . . . . .

5

2.4

Solution procedure of mixture modelling approach . . . . . . . . . . . 10

3.1

Fundamental types of pore diffusion . . . . . . . . . . . . . . . . . . . 15

3.2

Momentum interaction in multicomponent systems . . . . . . . . . . 16

3.3

Collision Scenario of two molecules . . . . . . . . . . . . . . . . . . . 17

3.4

Interaction of molecules in a control volume . . . . . . . . . . . . . . 18

3.5

Movement of molecules of a single phase in porous region

4.1

Flow chart of multiphaseEulerFoam solver . . . . . . . . . . . . . 30

4.2

Schematic to demonstrate MULES convective-only transport solution

4.3

Procedure to solve phase volume fraction

4.4

Solution procedure of phase continuity equation . . . . . . . . . . . . 35

5.1

Experimetal set up of loschmidt tube . . . . . . . . . . . . . . . . . . 45

5.2

Simulation model of loschmidt tube

5.3

Phase fraction ofArgon in left tube: Analytical vs. Calculated . . . . 47

5.4

Phase fraction of CH4 in left tube : Analytical vs. Calculated . . . . 48

5.5

Phase fraction of H2 in left tube: Analytical vs. Calculated

5.6

Calculation of Nusselt number . . . . . . . . . . . . . . . . . . . . . . 49

5.7

Case set up for a channel flow . . . . . . . . . . . . . . . . . . . . . . 50

5.8

Velocity profile of a non porous zone 2 . . . . . . . . . . . . . . . . . 51

5.9

Temperature profile of a non porous zone 2 . . . . . . . . . . . . . . . 52

. . . . . . 21

33

. . . . . . . . . . . . . . . 34

. . . . . . . . . . . . . . . . . . 46

. . . . . 48

5.10 Temperature and velocity variation (non porous zone 2) at x= 0.1m

52

5.11 Velocity profile for a non porous zone 2 . . . . . . . . . . . . . . . . . 53

iii

LIST OF FIGURES 5.12 Temperature profile for a porous zone 2 . . . . . . . . . . . . . . . . . 53 5.13 Temperature and velocity variation (porous zone 2) at x= 0.1m . . . 53 A.1 Initial diffusion of Argon . . . . . . . . . . . . . . . . . . . . . . . . . 60 A.2 Initial diffusion of Methane . . . . . . . . . . . . . . . . . . . . . . . . 61 A.3 Initial diffusion of Hydrogen . . . . . . . . . . . . . . . . . . . . . . . 62

iv

LIST OF TABLES

List of Tables 2.1

Comparison of mixture model and Eulerian model . . . . . . . . . . . 12

4.1

Phase properties and functions to access them . . . . . . . . . . . . . 31

4.2

Member function in multiphaseSystem . . . . . . . . . . . . . . . . . 31

4.3

Computing inter-component drag force by maxwellStefan() function . 39

4.4

Computing porous drag multiplier by porousSource.H . . . . . . . . . 41

4.5

Variables and functions required for energy transport . . . . . . . . . 42

4.6

Solution procedure for temperature of the mixture . . . . . . . . . . . 42

4.7

multicomponentPorousFoam algorithm . . . . . . . . . . . . . . . . . 43

5.1

Initial composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

v

Nomenclature

Nomenclature Roman Symbols m ˙

Mass transfer

kg

~υqx,por Velocity of phase q in porous region

ms−1

~υM q

Drift velocity of a phase from reference mass of mixture

ms−1

~qm

Effective conduction heat flux

cp

Specific isobaric heat capacity(fluid)

JKg−1 K−1

cs

Specific heat capacity(solid)

JKg−1 K−1

dp

Pore diameter

Jm−2 s−1

m

Fdrag Drag Force

N

Fq,por Porous Drag force of a phase q

N

Flift,q Lift force

N

Fvm,q Virtual mass force

N

T

Viscous stress tensor

Nm−2

A

Surface area

c

Mass fraction

D

Diffusion coefficient

d

Particle diameter

E

Energy

J

H

Height

m

h

Specific enthalpy

K

Permeability

k

Thermal conductivity

kn

Knudsen number

m2 1 m2 s−1 m

JKg−1 m−2 Wm−1 K−1 1

vi

Nomenclature L

Length

M

Molar mass

m

Mass

n

Amount of substance of the gas

Nu

Nusselt number

p

Pressure

R

Ideal gas constant

T

Temperature

V

Volume

m Kgmol−1 Kg mol 1 Pa Jmol−1 K−1 K m3

Greek Symbols α

Volume fraction

1

µ

Dynamic viscosity

ǫ

Porosity

1

ρ

Density

kg m−3

σpq

Molecular shock diameter

τ

Tortuosity

υ

Velocity

Pa s

m 1 ms−1

Subscripts p, q

Phases

1,2

Molecules

f

Fluid

m

Mixture

por

Porous region

Superscripts dar

Darcy

i

initial

k

Knudsen

vii

1. Motivation

1. Motivation Impact of soot is quite adverse on the environment. There is a need for the reduction of emission which contributes to soot . A recent study by the American Geophysical Union in their article “Bounding the role of black Carbon in the climate system: A scientific assessment” [4] gives an objective point of view and sheds light on the importance of soot reduction. It indicates that the effect of soot on the change in climate is twice as much as the previous estimates. It further states that fully mitigating soot can save 1 − 2 million lives by avoiding diseases caused from soot pollution. Complete statistics and details about the hazardous effects of soot can be viewed in the cited article. The use of Diesel Particulate Filters (DPF) to reduce emission is widespread in automotive and related industries. A first step towards soot reduction is to thorougly understand various phenomena occuring in a DPF. Diesel filter involves a multicomponent flow consisting of Particulate matter, Carbon-Di-Oxide, Oxygen, Nitrogen etc. The mixture encounters the filter which is a porous region. The presence of porous region further influences the motion of the components in the filter. Hence it is important to model the presence of the porous media with the mixture and capture the interaction effects among individual components and the components with the porous matrix. Approaches used to calculate gas transport is generally based upon parameterized and simplified models like the dusty gas model, mean transport model or binary friction model which predominantly depend on emperical parameters to model the effects of porosity. Another approach is to resolve the porous media and to consider the porosity effects. Though this is quite accurate in predicting actual flow phenomena, it is not generic in nature and is a trade off between having a generic model to represent porous media and computation time. Hence a need for an independant model is imminent which does not base itself upon emperical parameters to a large extent and is generic without compromising the accuracy of the solution. Simulating the behaviour of a multicomponent gas mixture in porous domains which is an abstract form of flow phenomena occuring in DPF, using Computational Fluid Dynamics (CFD) is undertaken in this thesis. The momentum exchange between individual components is modelled. Momentum source due to the presence of porous

1

1. Motivation structure in the domain is also accounted for. A semi-heuristic term for porous drag force which considers both Darcy effect and Knudsen diffusion in micro porous domains is used. Volume-averaging method of Slattery is used to model porous media. The article by Piesche and Goll [27] gives a direction to the modelling of the required terms. In this work, the development of the terms are explained and modelled. The modelled terms are then implemented into the open source field manipulation software OpenFOAM and a solver multicomponentPorousFoam is developed. Finally, the implementation of modelled terms and capability of the solver to simulate the flow scenario is substantiated.

2

2. Multi-component gas flow

2. Multi-component gas flow 2.1

Introduction

Transport of more than one fluid is imperative in fluid mechanics given the burgeoning applications involving many fluids interacting in a domain. Multiphase flows do not restrict the presence of different phases only in terms of solids, liquids and gases but is percieved in a broader sense. Eventually all the different phases are expressed in numbers in terms of density, viscosity, molecular weight, molecular diameter and other properties of individual phases. Hence, it is quite obvious that a group of particles having the same properties is represented as one phase rather than abiding by the conventional definition. This means that this group has similar dynamic responses. Some examples of multiphase flows are smoke from a chimney dispersing into the atmosphere, bubble flow due to aeration and blood flow. Further examples are shown in the schematic 2.1. Figure 2.1 gives an inital perspective of multiphase flows. In bubble flow, the presence of bubbles induce drag into the flow domain. In a slurry, the behavior of the dispersed phase is influenced by the continuous medium. Further, in this kind of flow, if the dispersed particles can be divided into two groups which largely differ in their respective diameters, then the number of phases must be treated as three instead of two since the change in diameter means a change in dynamic response by a group of identical particles. The spray dryer involves a fluid stream comprising of a solute and a solvent which is sent into a chamber containing hot vapour. This results in the vaporization of the fluid resulting in the separation of solute particles which is collected at the bottom of the chamber. These examples give an insight and emphasizes the importance of understanding multiphase flow phenomenon. Multiphase flows and multicomponent flows are two terminologies that varies in a trifle manner but demands to be distinguished. Multiphase flows involve fluids mixed above molecular level whereas multicomponent flows comprises fluids mixed at the molecular level. Multiphase flows require interface capturing. Multicomponent flows do not require interface capturing since various components are mixed at the molecular level. Exhaust gas treatment which involves the presence of diesel filter defined as a porous matrix defines the scope of this work.

3

2. Multi-component gas flow Bubble column

Slurry

Spray dryer Slurry

Hot vapour Air

Dried powder

Water

Figure 2.1: Examples of multiphase flows Porous media

Particulate Matter + CO + O2 + N2

Atmosphere

Figure 2.2: Multicomponent flow in porous media The figure 2.2 abstracts the whole process that occur in the treatment of exhaust gases. A multitude of gases and particulate matter mixed at the molecular level enters a channel. This multicomponent mixture encounters the porous media. The characterisitc property of a porous media is a randomised distribution of pores along space and varying pore diameters. At this juncture, it is important to note all the components flow through as one mixture and there is only a small flux due to molecular diffusion that drifts from the main mixture. “This is because of Knudsen numbers(kn) being greater than 1 × 10−3 because of elevated temperatures, low operating pressures and pore sizes on the microscale”[27]. The methods used to model multiphase or multicomponent flow is addressed in the following section.

4


2.2

Approaches to Multicomponent modelling

The two distinct treatment of multiphase flows are Lagrangian and Euler-Euler. The former is a particle tracking technique used when there is a dispersed phase interacting with a continuous media. It can largely be remarked that this method works well when there is a one way coupling between the interacting phases. Consider n particles with negligible mass being carried on by a fluid stream. This results in a one way coupling since the flow of fluid influences the particles but not the vice versa. We can go ahead and state that there is a momentum lost from the fluid through drag to the particles but not the other way around. Thus the whole process is controlled by the dominant continuous phase. But in flow domains where no phase is dominant, a need arises to formulate a method where all the phases involved are equally weighed. The Euler-Euler method models multiple phases with this approach. As formerly mentioned, the interest of this work deals with multicomponent modelling, we can directly decide to model it with the Euler-Euler method since there is no phase which is dominant when the components are mixed at the molecular level. The coupling between the phases is achieved by volume fraction. The Volume fraction of phases varies across space and over time. Simultaneous occupation of a phase by two fluids is not possible. Hence closure is obtained by a condition that the sum of volume fractions of all the phases to be unity.

2.2.1

Volume Of Fluid approach

αAir=1

αwater=1

Figure 2.3: Volume of Fluid approach Volume of Fluid (VOF) model is a surface-tracking technique applied to a fixed Eulerian mesh. “It is designed for two or more immiscible fluids where the position of the interface between the fluids is of interest”[1]. Surface tracking is done by evaluating the values at individual contol volumes. Consider two phases interacting as shown in the figure 2.3. A cell is considered to be completely filled by air when the phase fraction is one and completely empty when the phase fraction is zero.

5

2. Multi-component gas flow When the volume fraction is between unity and zero in a particular cell, then it is marked with the presence of interface. This demands an equation in order to solve the phase fraction of individual phases which is given by:

1 ρq

"

n

X ∂ (αq ρq ) + ∇ · (αq ρq υ~q ) = Sαq + (m˙pq − m˙qp ) ∂t p=1

#

.

(2.1)

The above equation 2.1 contains the substantive derivative of mass of individual phases. m˙pq is the mass transfer from phase p to phase q and m˙qp is the mass transfer from phase q to phase p. Sαq is source terms for mass of individual phases. Only volume averaged field variables are solved. Hence the individual phases share this common field variables proportionate to their volume fraction in all the cells. The volume averaged momentum equation is:

∂ (ρ~υ ) + ∇ · (ρ~υ~υ ) = −∇p + ∇ · [µ(∇~υ + ∇~υ T )] + ρ~g + F~ . ∂t

(2.2)

The energy equation shared by the phases is:

∂ (ρE) + ∇ · (~υ (ρE + p)) = ∇ · (kef f ∇T ) + Sh . ∂t

(2.3)

Energy and temperature are treated as mass averaged variables in the above equation 2.3. Sh is the heat source term. Energy shared by individual phases is proportionate to its mass fraction in the whole mixture. Surface tension can also be modelled by adding a source term to the momentum equation. Thus the VOF approach models the required effects in case of immiscible fluids where the need of interface tracking exists.

2.2.2

Mixture approach

The mixture modelling is one of the methods to model multiphase systems. In this approach, a single continuity equation and momentum equation is solved for the mixture. The mixture is assumed to be the averaged field properties of the field variables. Though the averaging is quite similar to the VOF approach, it has to be noticed that there was no equation solved for the conservation of mass in VOF approach since the interest was mainly in the interface. Contrastingly, in the mixture approach, the continuity equation for the mixture is solved. Hence it can model phases moving with different velocities. It is assumed that the individual phases with different velocities diffuse away from the whole mixture with a drift flux. This drift is represented as the relative velocity between the mixture and the individual phases. The generic equations for the mixture approach ([1],[18]) is explained subsequently . The continuity and momentum equations for a single phase q can be written as:

6


∂ (αq ρq ) + ∇ · (αq ρq υ~q ) = Γq , ∂t

(2.4)

∂ (αq ρq υ~q ) + ∇ · (αq ρq υ~q υ~q ) = −αq ∇pq + ∇ · [αq (τq + τT q )] + αq ρq g + Mq , (2.5) ∂t where Γq represents the rate of mass generation of phase q at the interface. Momentum exchange between the interface is represented by the source term Mq . τq is the average viscous stress tensor and τT q is the turbulent stress tensor. To obtain the equations for the mixture as a whole, all the equations of individual phases are added. The continuity equation for the mixture is: n

n

n

X X ∂ X (αq ρq ) + ∇ · (αq ρq ~υq ) = Γq . ∂t q=1 q=1 q=1

(2.6)

The mass in the whole system is conserved. Hence, the sum of Γq over individual phases must be zero; n X

Γq = 0 .

(2.7)

q=1

The continuity equation thus becomes: ∂ (αm ρm ) + ∇ · (αm ρm~υm ) = 0 , ∂t

(2.8)

where the mixture density is:

ρm =

n X

α q ρq ,

(2.9)

q=1

and mixture velocity is:

~υm =

n 1 X αq ρq ~υq . ρm q=1

(2.10)

The drift flux needs to have a reference velocity which is represented by the mixture velocity at the mass centre. The momentum equation for the mixture is obtained by adding the momentum equation of all phases:

7


n

n

n

n

X X X ∂ X (αq ρq ~υq ) + ∇ · (αq ρq ~υq ~υq ) = −αq ∇pq + ∇ · [αq (τq + τT q )] ∂t q=1 q=1 q=1 q=1 +

n X

α q ρq g +

q=1

n X

Mq . (2.11)

q=1

The second term on the left hand side equation 2.11 can be rewritten in terms of mixture density and mixture velocity:

∇·

n X

αq ρq ~υq ~υq = ∇ · (ρm~υm~υm ) + ∇ ·

n X

αq ρq ~υM q ~υM q .

(2.12)

q=1

k=1

Here ~υM q is the drift velocity of a phase from the reference mass of the mixture;

~υM q = ~υq − ~υm .

(2.13)

The momentum equation in terms of mixture variables is then:

∂ (ρm~υm ) + ∇ · (ρm~υm~υm ) = −∇pm + ∇ · ((τm + τT m )] − ∂t n X ∇· αq ρq ~υM q ~υM q + ρm g + Mm ,

(2.14)

q=1

where τm , τT m represent the average viscous and turbulent stress because of the slip of a phase. The pressure of the mixture is:

∇pm =

n X

αq ∇pq .

(2.15)

q=1

But more often the partial pressure and the mixture pressure is taken as the same barring a few cases. The term Mm is the contribution of the sum of surface tension from individual phases. The formulation for relative(slip) velocity is needed since the continuity and momentum equations are expressed in terms of it following which the continuity equation for phase fraction can be solved. In this approach, an algebraic formulation for slip velocity is constructed. Density differences of individual phases has a major influence on particles with different physical properties resulting in a varied dynamic response. The relative velocity between two phases q and j is further coupled to the drift velocity by:

8


n X

~υM q = ~υqj −

(ci~υji ) ,

(2.16)

i=1

where ci is the mass fraction of a phase;

ci =

α i ρi . ρm

(2.17)

The algebraic formulation for the relative velocity(slip) between two phases generally considers that the mass of a particle in the mixture is very small and that the particle reaches a terminal velocity over a short distance. In short, a local equilibrium is assumed. The additional force due to the relative velocity of the particle with respect to a fluid is represented by the drag force on the particle by the fluid. Thus, the relative velocities of a particular phase is obtained which can be used in the phase continuity and momentum equations subsequently. The expression of Manninen et al.[18] is generally used:

~υqj =

τq Fdrag

ρq − ρm α ~ ρq

,

(2.18)

where α is the secondary phase particle acceleration:

α ~ = ~g − (~υm · ∇)~υm −

∂~υm ∂t

(2.19)

and τq is the particle relaxation time ρq d q 2 τq = , 18µq

(2.20)

where dq is the particle diameter of phase q. The continuity equation for a phase which is incompressible and for a process where no phase changes occur is: ∂ (αq ) + ∇ · (αq ~υm ) = −∇ · (αq ~υM q ) . ∂t

(2.21)

An overview of the complete solution procedure is shown in figure 2.4. Thus, the mixture approach sums up the continuity and momentum equations of individual phases and circumvents solving for the velocity of phases by algebraic formulation for the relative velocity between phases and coupling drift velocity to relative velocity. It assumes local equilibrium over short spatial scales.

9

2. Multi-component gas flow timen

find slip velocities from (2.18)

find the phase fraction of all phases

Update mixture density (2.9)

find the drift velocity from drift velocity-slip velocity coupling (2.16)

solve for velocity-pressure using momentum(2.14) and continuity(2.8) equations iteratively

End time

Figure 2.4: Solution procedure of mixture modelling approach

2.2.3

Eulerian approach

Eulerian approach relies on solving the momentum and continuity equations for all the phases. “In the Euler-Euler approach, the different phases are treated mathematically as inter-penetrating continua”[1]. The coupling between the phases is established by an interaction term in the momentum equation of each phases which contains the relative velocity of a phase with respect to all the other phases. Unlike the former approaches, a number of phases can be modelled which can move with a particular velocity. The concept of volume fraction is considered and a particular space cannot be occupied by two phases. The volume occupied by a phase is proportionate to its phase fraction similar to the mixture theory approach. Though continuity and momentum equations are solved for individual phases, the effect of all these phases has to be considered in the multiphase framework. Hence, characteristic properties of viscosity and density of a particular phase have to be modified to fit the multiphase framework. It is straightforward to predict that the effect of physical property from a phase in the multiphase domain is proportional to its phase fraction. Hence the effective density used in the equations are:

10


ρˆq = αq ρq

(2.22)

µˆq = αq µq .

(2.23)

and viscosity

If all terms related to mass transfer across interfaces are set to zero, the continuity equation for a phase is

∂ (αq ρq ) + ∇ · (αq ρq υ~q ) = 0 . ∂t

(2.24)

The momentum balance equation is:

n

X ∂ ~ pq ) (αq ρq ~υq ) + ∇ · (αq ρq ~υq ~υq ) = −αq ∇p + ∇ · T + αq ρq~g + (R ∂t p=1 +(F~q + F~lif t,q + F~vm,q ) .

(2.25)

~ pq is modelled as an interphase drag term which is a The interaction term for R function of relative velocity of the phases:

n X p=1

R~pq =

n X

Kpq (υ~p − υ~q ) .

(2.26)

p=1

Kpq here is the interphase momentum exchange coefficient. Lift and virtual mass forces are not in the scope of this work and will not be discussed exhaustively. Various models of drag, heat transfer or molecular diffusion between the phases can ~ pq . The phase be modelled using the generic interphase momentum transfer term R fraction is solved from the continuity equation and the effective densities of phases are updated. Eulerian modelling is used in complicated flow phenomenon where the formerly stated approaches do not work. The complexity involved is the courtesy of allowing the phases to have its individual velocities and maintaining the boundedness of phase fractions. Additionally, it is computationally expensive since n continuity and momentum equations are solved in every iteration and demands larger memory requirements. The equations in this section are referred from [27] and [1].

11

2. Multi-component gas flow Feature Velocity

Mixture model mixture velocity

Pressure

mixture pressure

Coupling between phases

relative velocity between phases, drift velocity, mixture velocity

Wall slip modelling

mixture slip at wall

Computation reasonably same cost as single phase model

Eulerian Model individual phase velocities shared by all the phases

interphase momentum exchange term which accounts for the slips between phases phase slip at wall

increases with the number of phases

Remark

mixture pressure equal to partial pressure in most cases solving for phase velocity circumvented and never solved in mixture approach velocity of phase known in Eulerian model Eulerian model demands large memory requirements

Table 2.1: Comparison of mixture model and Eulerian model

2.2.4

Model comparisons

A suitable multiphase modelling approach has to be selected which suffices and represents all the processes that occur in case of flow through the porous media. It can directly be analysed that VOF approach is out of the purview of our work since interface tracking in not an important requirement. Now we are left with two approaches which needs diligent comparison. Table 2.1 shows the capabilities of the two modelling approaches. Though the mixture approach seems an attractive option for our flow scenario, it assumes local equilibrium between the phases and the relaxation time of a particle must be lesser than the timescale of the problem. This means that a strong coupling between the phases must be present. A sudden acceleration of a few phases may occur as the flow encounters porous media which might weaken the coupling. Another drawback of using mixture approach for our flow scenario is the requirement of the slip boundary condition for velocity at the wall. “The usual boundary condition, that the velocity of a gas bounded by a wall should be equal to that of the wall at every point of its surface, is not exactly fulfilled when there exists a velocity gradient perpendicular to the wall, a temperature gradient parallel to the wall or-in case of a mixture- a concentration gradient parallel to the wall ”[17]. Both temperature and concentration gradients exist parallel to the wall as a result of porous media in the flow domain. The porous media has its characteristic thermal conductivity which is different than the multiphase mixture. Hence there is a temperature gradient over the porous domain parallel to the wall. Existence of concentration gradient parallel to the wall can be attributed to the permeability of the porous media. Because of the permeability and a knudsen number greater than

12

2. Multi-component gas flow 0.01 in the porous domain amounts to the presence of varied mass across its sides. There is relatively more mass towards the inlet than the outlet. The aforementioned reasons establishes the fact that the mixture theory though with its lower computation time fails to fulfill the present requirements. Though with its considerably higher computation time, we proceed with the Eulerian multiphase model as it is capable of modelling the desired effects.

2.3

Closure

The three main approaches of multiphase modelling; VOF, mixture theory and Eulerian were introduced in this chapter. A broader perspective of considering multiphase media and the slight distinction between multiphase and multicomponent flows were discussed. Though this work concerns multicomponent flows, the term multiphase has been used to explain the approaches in order to retain the semantics of prevalent literatures on the same. Further, the characteristics of all the approaches were evaluated and compared. The Eulerian multiphase modelling was selected for modelling our flow problem with reasons stated.

13

3. Mathematical Modelling

3. Mathematical Modelling Mathematical modelling is imperative in order to express a physical system in a materialized form. A mathematical representation of systems or processes generates a possibility to study and analyse the influence of different components in the system. It starts with transforming the various processes and interactions into governing equations. The governing equations have to be modelled in such a way that, it must posses the capability to express all the important phenomena that may occur. The equations must also convene the application of boundary conditions to realistically simulate the physics of the system. In fluid systems, the continuity or the mass conservation, momentum conservation and the energy conservation are the three fundamental forms of equations that require mathematical modelling. Differential equations are used quite often in order to express such processes where there is a continuous variation of a quantity over space and time. A solution of the differential equation in terms of time and domain thus results in a possibility of monitoring the behaviour of a quantity over the whole domain and its variation with time. Most of the mathematical systems results in non-linear equations. The non-linearities pose some difficulties in finding an analytical solution to the modelling systems. The Navier-Stokes equation models the mass and momentum conservation of a fluid system and will be used in this work. The momentum interaction and porous drag terms are developed and are considered as source terms in the momentum equation. Energy equation for the pseudo homogeneous medium is modelled where the mixture and the porous matrix is assumed to have the same temperature and a local thermodynamic equilibrium exists between the individual phases. The subsequent sections deals with the development of mathematical models to express the flow phenomena of our system.

3.1 3.1.1

Drag modelling Effects of Diffusion

Diffusion is a prime mechanism responsible for a variety of processes involving mass transfer. This broad spectrum involves bulk gas or liquid diffusion, knudsen diffusion

14


Bulk diffusion

Knudsen diffusion

Micropore diffusion

Figure 3.1: Fundamental types of pore diffusion in pores and molecular diffusion inside micro pores. The driving force for diffusion is ascribed to the concentration gradients that exist in space. The tendency is always to flow from a region of higher concentration to a region of lower concentration. The subject of interest here is porous diffusion. To distinguish between different kind of porous diffusion mechanisms, we need to first distinguish between the different kinds of pores which are macro pores(size > 50nm), micro pores(size < 2nm) and meso pores(size between 2nm-50nm). The three types of diffusion mechanisms are shown in the figure 3.1 referred from [24]. In high pressure systems the number of molecules are so high that collision occurs more frequently amongst the molecules than with the wall. Large pore sizes also accommodate for such free molecular collisions. This kind of diffusion is termed as Bulk diffusion. Knudsen diffusion arises when the mean free path of a molecule is lesser than the charateristic length of the domain. Hence, the molecule travels lesser than the diameter in case of a circular tube before collision which results in a high probability of molecule-wall collisions. Surface diffusion occurs primarily in micro pores when the molecular species are adsorbed on the surface of porous media. “Bulk and Knudsen diffusion occur together and it is prudent to take both mechanisms into account rather than assume that one or other mechanism is ‘controlling’ ”[24]. The physical meaning of interaction between molecules of different components in a multicomponent system can be related to drag effects. The interpretation is that, molecules travelling with different velocities in a system exchange momentum among them due to relative velocity. It is visualised in the figure 3.2. The inter component

15

3. Mathematical Modelling U1 1

friction 1-2

U2

friction 12

friction 2-

Figure 3.2: Momentum interaction in multicomponent systems momentum exchange is directly proportional to the relative velocity. A suitable model which considers the stated effects of diffusion and replicates the appropriate interaction effects between components must be selected and integrated to the multicomponent framework. Diffusion of ideal gas mixtures is considered in this work.

3.1.2

Maxwell-Stefan relations for modelling drag

Consider a collision scenario of two molecules with mass m1 and m2 travelling with velocities u1 and u2 respectively. The assumption of collision is completely elastic. There is no loss of energy during the process. The other important assumption is that the molecules are spherical in shape. An example of such a process is collision between two billiard balls. The two molecules carry a momentum of m1 υ1 and m2 υ2 . If υ1 ′ and υ2 ′ are velocities after the collision and the path traced by the particles before and after collision is assumed to be the same, then according to the law of conservation of momentum:

m1~υ1 + m2~υ2 = m1~υ1′ + m2~υ2 ′ ,

(3.1)

m1 (~υ1 − ~υ1′ ) = −m2 (~υ2 − ~υ2′ ) .

(3.2)

16

3. Mathematical Modelling Plane of collision Path of collision

V1

V1'

1

1

Before collision

V2'

V2

2

After collision

After collision

2

Before collision

Figure 3.3: Collision Scenario of two molecules m1 (~υ1 −~υ1′ ) is the momentum transferred from molecule 1 to 2. The average velocity after collision υ1′ in such conditions is the velocity of the center of mass of the pair of molecules[25].

uc =

m1~υ1 + m2~υ2 = ~υ1′ . m1 + m2

(3.3)

The momentum transferred from molecule 1 to molecule 2 is:

m1 (~υ1 −

~υ1′ )

= m1~υ1 −

m1~υ1 − m2~υ2 m1 + m2

=

m1 m2 (~υ1 − ~υ2 ) . m1 + m2

(3.4)

The inference from equation 3.4 is that momentum is transferred only between molecules of different components. The amount of momentum transfer is proportional to the relative velocity between components. In short no intra-component but only inter-component momentum transfer takes place. With the mechanics of collision understood, we move on to the Maxwell-Stefan relations. In the control volume shown in the schematic 3.4, movement of molecules of type 1 is obstructed by molecules of type 2 which can be interpreted as the drag effect. Momentum transfer between the phases take place in such interactions. The application of Newton’s second law implies that Sum of forces ∝ Rate of change of momentum in the system The momentum exchange between molecules depends on the frequency of their collisions. It is obvious that the probability of molecules of two components increases

17


dz

Figure 3.4: Interaction of molecules in a control volume with the number of molecules present. For example, if 10 moles of Oxygen, 20 moles of Nitrogen and 4 moles of Argon is present in a contol volume, then the frequency of collision between Oxygen and Nitrogen is more extensive than other pairs of collision. The collision frequency thus depends on the phase fraction of the two components. Collision frequency between 1 − 2 ∝ α1 α2 Assuming that the control volume moves with the flow at the same velocity, the net flux becomes zero. Then in a system of unit volume in unit time, for molecules 1 and 2 the rate of change of momentum of molecule 1 is proportional to the product of average momentum exchanged in one collision and the collision frequency. Applying force balance on molecules of type 1 in the control volume and neglecting gravity forces and forces due to the gradient of velocity and considering only force due to pressure on the system:

− ∇p1 ∝ α1 α2 (~υ1 − ~υ2 ) ,

(3.5)

where p is the absolute pressure. The proportionality constant which is the drag coefficient f12 is introduced to the equations:

− ∇p1 = f12 α1 α2 (~υ1 − ~υ2 ) . f12 is defined as an inverse drag coefficient and D12 = p/f12

18

(3.6)


− ∇p1 = p

α1 α2 (~υ2 − ~υ1 ) α1 α2 (~υ1 − ~υ2 ) = −p . D12 D12

(3.7)

System pressure is often taken as constant. Hence,

− ∇p1 = −p∇α1

(3.8)

p1 = α1 p

(3.9)

since

The final relation becomes:

− p∇α1 = −p

α1 α2 (~υ2 − ~υ1 ) . D12

(3.10)

α2 α1 (~υ1 − ~υ2 ) . D21

(3.11)

The force balance for component 2 is,

− p∇α2 = −p

Extending the relation to a multi component system,

− p∇αq = −p

n X αq αp (~υp − ~υq )

Dqp

p=1

p∇αq = p

n X αq αp (~υp − ~υq )

Dqp

p=1

.

.

(3.12)

(3.13)

Thus to implement Maxwell-Stefan drag into the multicomponent model, equation 3.13 has to be added to the momentum equation. The pressure terms cancel with each other but it is convenient to have it in this form to synchronize with dimensional units in the momentum equation. The procedure of drag modelling in the section is referred from [25] and further details can be found in the same.

19


3.2

Porosity Modelling

The presence of porous media in the domain offers resistance to the flow because of drag effects. Porous structures are characterized by voids or pores through which the fluid flows, labyrinth or tortousity factor which defines interconnectivity of the porous structure and porosity which is the ratio of the amount of fluid present in the porous structure to the total volume including fluid and solid volumes. It can be directly deduced that the drag force increases with increase in the surface area of the porous matrix. The surface area of the matrix can further be divided into two parts. Area normal to the fluid flow and internal surface area of the matrix structure which induces viscous drag on the fluid flowing through it. In the experiment conducted by Darcy, it was realized that the bulk shear stresses were dominant as the length of the porous media and hence the internal surface used was many times greater than its normal surface area. Darcy then found the velocity by measuring flux of an incompressible of the fluid in the porous matrix ~υqdar x,por and dividing it by fluid flowing through a cylindrical tube at the outlet ρA~υqdar x,por ρA. Pressure difference was the driving force which was directly proportional to . He realized that the the flux and hence the velocity in the porous region ~υqdar x,por proportionality constant could be expressed in terms of permeability K of the porous media and viscosity µ of the fluid. The darcy relation is [20]:

∂p µ , = − ~υqdar ∂x K x,por

(3.14)

K ∂p . µ ∂x

(3.15)

which is rewritten as:

=− ~υqdar x,por

“Permeability accounts for the interstitial surface area, the fluid particle path as it flows through the matrix and other related hydrodynamic characteristics of the matrix ”[20]. In micro porous matrices, the flux at the outlet is generally more than those measured by Darcy’s experiment since an additional Knudsen diffusion flux appears because of slip flow of the molecules with the porous structure walls. Hence the additional knudsen flux has to be modelled to appropriately represent the flow condition. Figure 3.5 shows different forms of transport of molecules in a porous matrix. As seen, momentum is lost by the molecules due to their interaction with the surface of the porous matrix. This is termed as the Darcy drag. A few molecules which move in micro pores experience knudsen slip as mentioned earlier contributing to the additional Knudsen flux. The remaining molecules retrace their path because of the inter connected pores in the porous matrix. This is accounted for by the tortuosity term. Consider a fluid q flowing through a porous matrix. From the kinetic theory of gases, the mean molecular velocity for the component can be written as [20]:

20


Px

Px+dx

Molecules that undergo Darcy drag from internal surface area of porous matrix Molecules that are transported due to Knudsen slip Molecules that are retraced due to labyrinth/tortuosity factor of the porous matrix

Figure 3.5: Movement of molecules of a single phase in porous region

~υqm,por =

s

8RT . πMq

(3.16)

R is the universal gas constant, T is the temperature in Kelvin and Mq is the molecular weight in Kg mol−1 . Considering a flow through a tube which is the idealised form of a pore with diameter dp and applying momentum conservation, Shear stress at the radius of the pore is[20]:

~υqk 3 dp dp τ(dp /2) = nq kB T x,por = − . 4 ~υqm,por 4 dx

(3.17)

From the above equation, an expression for ~υqkx,por is obtained,

~υqkx,por = −

4 dp dp ~υq . 4 dx 3nq kB T m,por

(3.18)

Ideal gas behavior is assumed:

pq = nq kB T .

(3.19)

Substituting ~υqm from equation 3.16:

~υqkx, por

dp =− 3pq

s 21

8RT dp , πMq dx

(3.20)

3. Mathematical Modelling which can be rewritten as,

~υqkx,por

Dq k dp , =− pq dx

(3.21)

s

(3.22)

where

Dq k

dp = 3

8RT πMq

is the Knudsen diffusion coefficient of a gas component. Thus for flows in 1-dimension and from 3.21 , 3.15 we have:

~υqx,por = ~υqdar + ~υqkx,por , x,por

~υqx,por = −

Dq k dp K ∂p − , pq dx µ ∂x

(3.23)

(3.24)

and extending it to 3-dimensions, ~υq,por

[Dq k ] [K] =− ∇p . ∇p − pq µ

(3.25)

In this work, porosity is modelled for non deformable and isotropic media. The permeability expression for permeability K is modelled according to the velocity profile of a Hagen-Poiseuille flow[27]. It is dependant on the diameter of the pore and is:

K=

dp 2 . 8

(3.26)

The dimensionless tortuosity or the labyrinth factor τ is introduced into the equations to compensate the extended path a fluid travels when it encounters closed paths due to interconnectivity of the porous media. Thus the expression for pressure gradient which is the driving force for the flows in porous media is:

− ∇p = τ

2

Dq k K + pq µ 22

−1

~υq,por .

(3.27)

3. Mathematical Modelling In order to model this momentum source due to porous structure in the fluid equations, local volume averaging is used. It is a technique introduced by Whitaker[28] and Slattery[13]. The approach is to define a quantity φ which is in a fluid volume Vf by averaging over a volume V. The condition for averaging is that the volume over which φ is averaged must neither be large compared to the flow domain nor very small to represent the geometry of the porous structure. “The local representative elementary volume is chosen such that the smallest differential volume that results in statistically meaningful local average properties (such as local porosity) and when this volume is appropriately chosen, adding extra pores(extra volume) around it will not result in changes to the values of these local properties”[20]. A penalty approach is used in order to define the void distribution fuction or the solid distribution function. Let us name this distribution function as gp defined as,

gp (x) =

1 if x is in void region,Vf 0 if x is in solid region,Vs

(3.28)

Local porosity is defined as:

1 ǫ(x) = V

Z

Vf V

(3.29)

φdV = ǫφ .

(3.30)

gp (x)dV = V

where V = Vf + Vs . Volume average of a quantity φ is, 1 < φ >= V

Z

Vf

A comprehensive explanation of volume averaging including the general form of Reynolds transport theorem can be found in [20] and will not be explicitly mentioned here. Thus by applying the stated averaging technique to the modelled source porosity terms in equation 3.30 the final term that is added to the fluid momentum equation is:

Fq,por = −ǫ∇p − gp τ

2

Dq k K + pq µ

−1

~υq,por .

(3.31)

The first term on the right hand side can replace −∇p in the Eulerian momentum equation since ǫ is 1 in the absence of porous media and accounts for both porous regions and non porous regions.

23


3.3

Energy transport

Energy transport is modelled by the mixture approach. In case of non-isothermal flow conditions, the change of temperature over space and time has to be considered. Local thermodynamic equilibrium is assumed so that a single energy equation is sufficient for the gas mixture as a whole[22]. Consider the general Reynold’s transport equation for a variable φ for a specific mass: ∂(ρφ) + ∇ · (ρ~υ φ) = ∇ · (γ∇φ) + Sφ . ∂t

(3.32)

It comprises the unsteady term, convective term, diffusive and source terms. The variable in case of energy transport is the enthalpy. Since we are solving the energy equation for the mixture, the equation becomes[27]: ∂p ∂(ρm hm ) + ∇ · (ρm hm~υm ) = −∇ · (~qm ) + . ∂t ∂t

(3.33)

Heat conduction is written normally in terms of mass averaged velocity of the mixture[29]: n 1 X ~υm = ρq ~υq . ρm q=1

(3.34)

The effective heat conduction of the mixture is:

~qm = −km ∇T +

n X

ρq hq (~υq − ~υm ) .

(3.35)

q=1

km is the molecular thermal conductivity of the gas mixture[22] and is obtained by:

km =

n X

αq kq .

(3.36)

q=1

αq and kq are the phase fraction and thermal conductivity of a component respectively. Now substituting 3.35 in 3.33: n

X ∂(ρm hm ) + ∇ · (ρm hm~υm ) = ∇ · (km ∇T ) − ∇ · (ρq hq ~υq ) ∂t q=1 +∇ · (ρm hm~υm ) +

24

∂p ∂t

(3.37)


n n X X (ρq hq ) = (αq ρq i )(αq hq i ) = ρm hm q=1

(3.38)

q=1

The above equation then reduces to the following form: n

X ∂(ρm hm ) ∂p +∇· . (ρq hq ~υq ) = ∇ · (km ∇T ) + ∂t ∂t q=1

(3.39)

Further the enthalpy can be related to temperature by:

dhq = Cp,q dT ,

(3.40)

where Cp,q is the specific heat of a component. Then equation 3.39 transforms into a scalar transport equation by using the relation from equation 3.40: n

X ∂p ∂(ρm Cp,m T ) +∇· . (ρq Cp,q ~υq T ) = ∇ · (km ∇T ) + ∂t ∂t q=1

(3.41)

“When porous media is considered, in the absence of local heat sources or high frequency temperature changes, local thermodynamic equilibrium between gas and solid can be assumed ” [20]. This in the concept of volume averaging implies that in the averaged volume, fluid and solid both have the same temperature. To describe the transport of energy in such cases, a single equation is setup which considers the pseudo homogeneous medium. The temporal changes of local enthalpy and heat conduction have to be expressed in terms of pseudo-homogeneous gas/solid medium [27]: 1 < ρCp T >= Vf

Z

1 (ρm Cp,m T )dV + Vs Vf

Z

(ρs Cs T )dV Vs

= [ǫ(ρm Cp,m ) + (1 − ǫ)(ρs Cs )]T .

(3.42)

By the application of volume averaging of Slattery’s theorem[20] of volume averaging: km − ks < ∇ · (k∇T ) >= ∇· < k∇T > + V

Z

~nT ′ dA .

(3.43)

Af,s

Af,s is the fluid-solid interface area. The second part on the right hand side of the above equation demands information on geometry of porous structure and is hence neglected following [15, 19] and [2]. The first term of equation 3.43 is then similarly transformed by following equation 3.42:

25


< k∇T >=< k > ∇T = [ǫkm + (1 − ǫ)ks ]∇T

(3.44)

The final equation for the pseudo-homogeneous medium is: ∂ [ǫ(ρm Cp,m ) + (1 − ǫ)(ρs Cs )]T + ∇ · ∂t

"

n X

(ρq Cp,q T ~υq )

q=1

= ∇ · [ǫkm + (1 − ǫ)ks ]∇T + ǫ

#

∂p . ∂t

(3.45)

Equation 3.45 is used to model energy transport in the present work. Following the solution of the equation, temperature dependant properties of viscosity µ and Diffusion coefficient Dpq in the system are updated using the Chapman-Enskog relation. 5 1 µq = 16 σq2 ΩV

r

kB Mq T . πNA

(3.46)

νq is the viscosity of a component, σq is the molecular shock diameter, NA is the Avagadro number and kB is the Boltzmann constant. ΩV is the collision integral for viscosity based on Lennard-Jones-Potenzial[3]. Dpq is the diffusion coefficient between two components: 3

Dpq

RT 2 3 = 2 p 8σpq ΩD

s

kB . πMpq NA

(3.47)

In equation 3.47, ΩD is the collision integral for diffusion. Mpq and σpq are the effective molar mass and molecular shock diameter of two components respectively. They are given by: Mpq = 2(Mp −1 + Mq −1 )

σpq =

σp + σq . 2

−1

,

(3.48)

(3.49)

The density of the component can also be updated by using the equation of state for ideal gases:

ρq =

pMq . RT

(3.50)

Thus the energy transport model takes into account all the properties that change with temperature. The implementation of the equation into the CFD solver will be discussed in the following chapter. It can be summarized that though the momentum and continuity equations are solved using the Euler-Euler approach for reasons mentioned earlier, energy equation is solved by the mixture approach in this work.

26


3.4

Closure

In this chapter, the effects of diffusion were stated and the drag term was modelled using Maxwell-Stefan relations. A flow scenario of molecules in a control volume was used to model the additional driving force due to concentration gradients of a phase and the drag term. Both terms were modelled and the final term Fdrag from equation 3.13 has to be added to the momentum equation of the Eulerian model. The effect of porous media in the flow domain was treated by considering both the viscous drag or Darcy drag effects and the addition Knudsen flux due to slip of the molecules with walls of micro porous domains. The modelled porosity term can represent macro properties of micro porous domain. Further, volume averaging was used in order to represent averaged properties of the pseudo-homogeneous media of solid and gas. Permeability and Knudsen diffusion coefficient are considered for isotropic porous media though the derived momentum source terms due to porous effects Fq,por equation 3.30 is inversed and can be easily extended to other types of porous structures. The Knudsen diffusivity Dq k depends on the ideal pore diameter dp (3.24). Permeability K which also depends on the pore diameter and tortuosity is modelled semi-heuristically. Finally, the energy transport for the whole domain is modelled by volume-averaging (3.45). The mixture approach is used for solving the equation by assuming local thermodynamic equilibrium between the phases. The solid and liquid phase present in the averaged volume that represent porous media share the same temperature. Temperature dependance of viscosity (3.48), diffusion coefficient (3.49) and density (3.50) are also considered in the modelling framework.

27

4. Numerical Simulation

4. Numerical Simulation Computational methods for Fluid Dynamics (CFD) offers many possibilities and flexibility to envision the influence of various parameters on the flow. Experimental methods though sturdy, do not offer this and takes a longer time and higher costs. “While the first computers were built in 1950s performed only a few operations per second, machines are being designed to produce 1012 floating point operations per second ”[14]. With computing power increasing, it is not a surprise that CFD has found its popularity for flow simulations. The previous chapter on mathematical modelling contained many differential equations which are mainly used to describe transport phenomena. Analytical solution to these differential equations is restricted to simple flows. Hence numerical methods must be used to discritize and solve complicated differential equations. The three main approaches for discritization are Finite Element Method(FEM), Finite Difference Method(FDM) and Finite Volume Method(FVM). All three methods follow the saying “How do you eat an elephant - one bite at a time”, by dividing the whole domain into smaller parts and solving the discritized form of the differential equation for every part and in turn finding the solution for the whole domain. FVM is the most commonly used discritization technique in CFD. The small “bites” are defined by a control volume(C.V) and two C.Vs share one common face among them. The calculation domain is divided into a number of non overlapping control volumes so that there is one control volume surrounding each grid point[21]. This works well for flow simulations as volume based properties like pressure, velocity can be distinguished with face centered properties such as flux. The discritized equation constructed in each C.V is of the form: A~x = ~b

(4.1)

where A is a coefficient matrix, ~x is the unknown vector and ~b is the source vector. The value of unknown variable is known at all locations of the domain by interpolation. The interpolation schemes can be linear or higher order schemes based on requirements.

28


4.1

Open Source CFD: OpenFOAM

OpenFOAM is an open source software used for Field Operation And Manipulation. It has a good potential to solve CFD problems as it houses many algorithms and interpolation schemes to solve discritized equations and interpolation schemes. Standard Dirichlet and Neumann along with more commonly used boundary conditions are built in the software and can be directly used. Coupling betwen pressure and velocity is handled by algorithms such as SIMPLE, PIMPLE and PISO which use the iterative approach for solution. Object Oriented Programming is the language used for programming and hence gives us enormous possibilities to extend and customize or build even a solver to suit our requirements. “Majority of fluid dynamics can be described using the tensor calculus upto rank 2, i.e. scalars, vectors and second-rank tensors” [9]. Tensors are stored as lists such as Field, Field, Field which are termed as scalarField, vectorField and tensorField for easy identification. Volume centred and surface centred properties can be further distinguished by adding a prefix ‘vol’ or ‘surface’ such as volScalarField in case of pressure and surfaceScalarField in case of flux. Storing a variable as a list offers an easy way to operate or manipulate variables with each other. Focusing on developing a solver in OpenFOAM for our flow scenario, the solver multiphaseEulerFoam in OpenFOAM 2.1.1 which solves for systems with many incompressible phases is identified as a good starting point.

4.2

Solver: multiphaseEulerFoam

multiphaseEulerFoam in OpenFOAM 2.1.1 is based on the Eulerian method for solving multiphase flows. It treats individual phases as incompressible. The default mechanism used for the transport of individual phases is the single phase transport model based on the viscosity model. Unsteady terms are included in the equations and hence transient flow cases can be simulated. The presence of individual phases and the averaged mixture in the system demands a separate consideration during programming to deal with the properties of phases and mixture. multiphaseEulerFoam handles this requirement by two distinct libraries, phaseModel and multiphaseSystem. Many models for drag, heat transfer, kinetic theory are already present and can be used by invoking them from the case directory according to the requirements. Additionally, it handles every phase/component in an idealized way by defining it as spherical particles and the diameter model is used for this purpose. Turbulence for the mixture can be solved by using LES method. MULES algortihm solves for phase fraction. Pressure-velocity coupling is dealt by PIMPLE algortihm using a small variation to projection methods. “Methods which first construct velocity field that does not satisfy the continuity equation and then correct it by subtracting something (usually a pressure gradient) are known as projection methods”[14]. The variation and the method will be discussed in detail subsequently. Flow chart in 4.1 shows the main parts and working of multiphaseEulerFoam solver. The required attributes and methods are intialized with phaseModel, multiphaseSystem and createFields. Following initialization, the pressure-velocity coupling is

29

4. Numerical Simulation Interfacial Models:

phaseModel:

Library for drag and heat transfer models

Start

Library to define phase properties and functions to operate on individual phases

multiphaseEulerFoam multiphaseSystem: KineticTheoryModels: Library with functions to operate on fluid properties or calculations involving more than

Library for kinetic theory models

one phase

createFields: Define mixture properties

Pimple loop for Pressure Velocity correction: 1. Update mixture turbulence 2. Solve phase fraction using phase continuity equation:fluid.solve() 3. Update mixture density 4. Calculate phase volumes of zones when cell zones are present in the domain 5. Set up the velocity equation for individual phases 6. Solve for pressure using pressure corrector pimple loop 7. Reconstruct phase velocities

Update substantive derivative for individual phases

End

Figure 4.1: Flow chart of multiphaseEulerFoam solver handled by a PIMPLE loop. Turbulence for the mixture is first corrected and updated. Then, the phase fraction for individual phases is solved and the phase fraction for the current time step is updated. The phase continuity equation is solved using the old values of velocity. Density of the mixture which depends on phase fraction is then updated. Next step in the loop involves calculating phase volumes of a cell zone when required. The velocity for individual phases is now set up which do not satisfy the continuity equation. Then pressure is solved by a PIMPLE loop and velocity is reconstructed. The reconstructed velocities now satisfy the continuity equation. After executing the PIMPLE loop, the substantive derivative of individual phases is updated. Table 4.1 and 4.2 lists the attributes and methods available in the solver to operate on component and mixture based parameters. The list defines the operating space and potential of the solver and gives a direction to extend the variables and functions

30

4. Numerical Simulation required for our flow case. The solver is already validated for a few cases which can be found in the tutorials of OpenFOAM 2.1.1 and will not be discussed here. A detailed veiw on the solution method of the solver will be discussed next. Variable Name nu kappa Cp rho U phiAlpha alpha

Description Kinematic Viscosity Thermal Conductivity Heat capacity Density Velocity Flux Volume fraction

Access function nu()

Field Type

Units

dimensionedScalar

m2 s−1

kappa()

dimensionedScalar

kgms−3 K −1

Cp() rho() U() phiAlpha() alpha()

dimensionedScalar dimensionedScalar volVectorField surfaceScalarField volScalarField

m2 s−2 K −1 kgm−3 ms− −1 m3 s−1 1

D ∂ Substantive Derivative: Dt = ∂t + ~υ · ∇ All variables in the table are phase specific and always appears with the phase subscript. for example, nuq ,ρq . Viscosity of a phase q can be accessed by q.nu()

Table 4.1: Phase properties and functions to access them Name phases() rho()

Return type Pointer Dictionary volScalarField

nu()

volScalarField

Cvm(phase) volScalarField

Svm(phase) volVectorField

Function All the phases can be iterated Returns the mixture density Returns the mixture laminar viscosity Returns virtual mass coefficient for a given phase Returns the virtual mass- source for the given phase

Actual term

ρm = νm =

Pn

q=1

Pn

q=1

α q ρq αq nuq

Cvm = Pn q Cvm(q, i) ∗ α i i=1 Svm Pn q

= Cvm(q, i)α α i q i=1

Cvm(q,i) is user defined value obtained from the case directory. To access the properties, an instance of the class multiphaseSystem is created (mulitphaseSystem fluid) and then accessed (fluid.rho()). Surface tension and drag model properties are not mentioned as it is out of the frame of this study.

Table 4.2: Member function in multiphaseSystem

4.2.1

Solution Method: Phase fraction

Phase fractions are solved from the phase continuity equation.The solution procedure bases itself on the two fluid model developed by Weller [7]. Assuming there is no mass transfer across phases, the equation is of the form[27]:

31


∂ (αq ρq ) + ∇ · αq ρq ~υq = 0 . ∂t

(4.2)

Since each phase is incompressible, the operating equation can be rewritten as:

∂αq + ∇ · αq ~υq = 0 . ∂t

(4.3)

Continuity is enforced by the phase fraction of each components and demands a strong boundedness property for phase fractions. In his work, Weller[7] introduced the conservative form of the equation for a system of two fluids a and b:

∂αa + ∇ · αa~υ + ∇ · (~υr αa (1 − αa )) = 0 , ∂t

(4.4)

where ~υ = αa~υa + αb ~υb and ~υr = ~υa − ~υb . This approach also couples the two phases more implicitly through the presence of relative velocity ~υr in the third term[26]. The discritised form of the equation is:

∂[αa ] + ∇ · αa f (φ,S) φ + ∇ · αa f (φra ,S) φra = 0 , ∂t

where φra = αbf(−φ

r ,S)

(4.5)

φr and φr = φa − φb .

The general approach to solve the phase fractions of a multiphase system in multiphaseEulerFoam is similar but there are some differences in handling the equations, or it can be said that it is “evolving”. MULES algorithm handles the boundedness property by first limiting the flux transport and then solves for phase fraction. The limiting of flux transport is important since large transport of fluxes from a cell may drive the volume fraction in a particular cell below zero during one time step. Hence the flux is at first limited and then is passed on to the MULES solver to solve for phase fraction. The constructor of MULES function is given below: MULES::explicitSolve (const RhoType& rho, volScalarField& psi, const surfaceScalarField& phiPsi, const SpType& Sp, const SuType& Su) It is seen that density has to be passed as an argument to the constructor. In case of incompressible phases, it is passed as geometricOneField() which is a unit value field. The second argument psi is the variable to be solved which is phase fraction in our case. The limited normal convective flux is the next argument which is earlier

32

4. Numerical Simulation solved explicitly by MULES:limiter. The next two terms are the explicit and implicit source terms in the continuity equation which arise when mass transfer across the phases or reaction source terms exist. In our case we pass both the arguments as zeroFields() since there is no mass transfer. MULES algorithm solves for phase fraction with explicit consideration of convective flux of phase fraction. It must also be mentioned that, the transport mechanism considered is convective-only. Consider equation 4.3. The flux of phase fraction is written as:

φq = αq ~υq .

(4.6)

Equation 4.3 in numerical form is written as:

αq new − αq old + ∇ · φq = 0 , δt

(4.7)

αq new = αq old − (∇ · φq )δt .

(4.8)

But in equation 4.8, the right hand side has two different types of fields and the surfaceScalarField of flux has to be converted to an equivalent volScalarField in order to perform algebraic operations.

(αqUq)in

(αq)old E

(αqUq)out

Figure 4.2: Schematic to demonstrate MULES convective-only transport solution Figure 4.2 shows a cell volume in 2 dimensions. Let us consider the flow dominant only in one direction. At the new time step, we need to know the effective amount of αq , phase fraction added or removed from the cell. If the mechanism of transport over the cell faces is only convective, then the quantity of αq at each of the cell faces going out or coming in from the neighbouring cell is αq ~υq . The normal convective fluxes are already known from the MULES:limiter function. Hence the divergence of flux is known, which means that the effective amount of αq added or removed is known. Applying the Gauss divergence theorem and from [11], the divergence operator can be transformed into a sum of all faces by:

33


Z

∇ · αq dV = Vp

X Z f

(∇ · αq )Vp =

X

dS · αq f

,

(4.9)

S · α qf .

(4.10)

f

The theorem can be used to convert a surfaceScalarField back to a volScalarField. This means that by knowing the convective fluxes, precisely limited convective fluxes, we can find αq added for the new time. Thus phase fraction is updated for the current time using equation 4.8 solveAlphas()

Calculate the flux of the mixture : Φc

for all pha

q=0

Initialize limited bounded flux phiAlpha=phiAlphaCorr

for all pha Interpola Initialize phiAlpaq=0

Initialize phiAlphaCorr

Solve for MULES:explicitSolve() phiAlpaq)

phiAlphaq += phiAlpha

+=

for all pha

Calculate relative flux be phir=phiq-phip

q

End

Update phiAlphaCorr phiAlphaCorr= phirapaq

Limit the flux tran bounde Initialize

e MULES:limit() =0

Figure 4.3: Procedure to solve phase volume fraction The solution procedure is shown in figure 4.3 In the first step, relative fluxes of a phase with respect to all other phases are found out. Then fluxes for individual

34

4. Numerical Simulation phases are limited to ensure boundedness. Next step involves passing the arguments of limited fluxes(MULES:limit) formerly found out and phase to solve the phase fraction for the current time step using MULES:explicitSolve(). Th interpolation of fluxes from centre to face is done by Upwind Differencing(UD). Using UD might be quite diffusive, however using a higher order scheme instead reduces the numerical diffusion , but might compromise the boundedness of the solution[26]. Finally, all the phase fractions are solved for the current time step and sum of all the phases is calculated. )

alpha correction nalpha

for all pha

Set old alpha value αqold

q=0

Initialize

while na

Solve αq and

q

for all pha

=

q

end

Figure 4.4: Solution procedure of phase continuity equation Solution for phase fraction is invoked by fluidSolve() function. The schematic of the function operation is shown in figure 4.4. Correctors for phase fraction is read from the pimple directory in the case file. If the number of correctors are greater than one, then phase fraction at old time is stored and fluxes for all phases are set to zero. Non-orthogonality for phase fractions is solved by the function solveAlphas()4.3. Then the total contribution of fluxes from individual phases is calculated. If the

35

4. Numerical Simulation number of non orthogonal correctors is zero, then solve for phase fraction and total flux by calling solveAlphas().

4.2.2

Solution Method: Pressure-Velocity coupling

Projection methods are used to handle the pressure velocity coupling. The method involves setting up an equation for velocity by projecting out pressure related terms. At this stage, the velocities do not satisfy the continuity equation. Equation for pressure is then set up using the continuity equation and then solved. Since pressure is now known, the contribution of pressure term is removed from velocity to find the correct velocities which satisfy the continuity equation. This method works for an arrangement where pressure is stored at cell centres and velocity at cell faces; in case of staggered grids. But since OpenFOAM uses collocated grid arrangement where both pressure and velocities are stored at cell centres, appropriate changes have to be made for the procedure to work. In his work, Weller[8] reformulated the solution procedure for collocated grid arrangement where velocities are reconstructed from flux fields. The approach of momentum prediction and pressure correction transforms to flux prediction and pressure correction. It is because flux, which is stored at the cell face now becomes the primary variable in the solution procedure which represents velocity. “The cell centred velocity is merely regarded as a secondary variable which is used in the construction of momentum equation” [26]. The solution procedure of Weller[8] is as follows: Consider the semi-discritized momentum equation:

AU U = {b} − ∇p ,

(4.11)

where b is a vector of all source terms except pressure gradient. AU is the coefficient matrix of veclocity vector. If it is written as sum of diagonal and off diagonal elements as,

AU = AD + AN .

(4.12)

Considering the off diagonal elements explicitly, the discritized equation is of the form,

AD U = AH − ∇p ,

(4.13)

where AH = {b} − AN U . The velocity equation is:

U=

AH ∇p − . AD AD

The velocity at the face is then:

36

(4.14)


Uf =

AH AD

− f

∇p AD

.

(4.15)

f

Flux is obtained by multiplying velocity with the area of cell face. Hence velocity at the face is obtained by:

φ = |S|Uf .

(4.16)

Equation for flux is then set up in agreement with equation 4.15 as:

∗

φf = φf −

∇p AD

.

(4.17)

f

Applying the continuity equation on equation 4.17:

∇·

"

AH AD

− f

∇p AD

#

=0.

(4.18)

f

The final poisson equation for pressure is:

∇2 p AD

=∇· f

AH AD

.

(4.19)

f

Since the pressure is solved at this stage, flux is corrected using equation 4.17 and velocities are reconstructed using equations 4.16 and 4.15

4.2.3

Semi-implicit treatment of Drag term

Drag force by and large is a function of relative velocity. If ~υp and ~υq are velocities of phases p and q, drag force is written as:

FDrag = C(~υp − ~υq ) .

(4.20)

The above expression can be treated in three ways; explicit, semi-implicit and fully implicit. In multiphaseEulerFoam ,it is treated as semi- implicit source term. If n phases are interacting, the drag term for phase q is:

FDrag,q =

n X

Cpq (~υp − ~υq ) .

p=1

37

(4.21)

4. Numerical Simulation Cpq is the known coefficient of the drag term between two interacting phases. The above equation is rewritten as:

FDrag,q =

n X

Cpq (~υp ) −

p=1

n X

Cpq (~υq ) ,

(4.22)

p=1

FDrag,q = (FDrag,q )explicit + (FDrag,q )implicit .

(4.23)

Hence while setting up the velocity equation 4.15, the terms are split into two as shown above into fully explicit and fully implicit. The semi-implicit treatment of the whole drag term, does not affect the diagonal dominance of the coefficient matrix but in fact enhances it.

4.3

Modified solver: multicomponentPorousFoam

Source terms required to model multicomponent flows through porous media were discussed and derived in chapter 3. Maxwell-Stefan relation is used for modelling momentum transfer between different components. An additional source term to model porous drag has been added. Volume-averaged energy equation needs to be implemented in the solution procedure. Modifications and extensions to the existing solver needs to be done. The following are the equations of continuity and momentum in the existing multiphaseEulerFoam :

∂αq + ∇ · (αq ~υq ) = 0 , ∂t

(4.24)

∂αq ~υq + ∇ · (αq ~υq ~υq ) = −αq ∇p + ∇ · T + αq ρq~g . ∂t The energy transport equation was developed in section 3.3 and is rewritten here:

" n # X ∂ [ǫ(ρm Cp,m ) + (1 − ǫ)(ρs Cs )] T + ∇ · (ρq Cp,q T ~υq ) ∂t q=1 = ∇ · [ǫkm + (1 − ǫ)ks ]∇T + ǫ

∂p . ∂t

(4.25)

Modifications and extension to the momentum equation along with the implementation of the complete energy equation is done in the new solver multicomponentPorousFoam .

38


4.3.1

Source term: Maxwell-Stefan drag

The terms added to the momentum equation for modelling drag is derived in section 3.1.2 and is rewritten here:

Fdrag = −p∇αq + p

n X αq αp (υp − υq ) p=1

Dqp

= Sp1 + Sp2 .

(4.26)

Both the source terms Sp1 and Sp2 are not present in the existing equations 4.25 of the solver. Hence they are added as source terms to the momentum equations. Equation 4.26 is added to the momentum equation as:

− p ∗ fvc :: grad(αq ) + fluid.maxwellStefan(phase,p,Tfluid) .

(4.27)

Both terms are considered explicitly in this work though the second term can be easily extended to a semi-implicit treatment as shown in section 4.2.3. As seen in 4.27, the calculation of drag force is done by the function maxwellStefan.() which is a method of the multiphaseSystem and can be accessed by instances of the multiphaseSystem class (fluid.multiphaseSystem). The constructor takes three arguments which is the phase name, pressure and temperature. The execution of the function is enumerated below for phase q.

1. Read phase name, pressure and temperature. 2. Set FDrag,q = 0 3. Loop over all other phases p = 1 : n excluding q • Calculate temperature and pressure dependant diffusion coefficient for two interacting phases Dpq from equation 3.47. • Calculate Relative velocity between p and q • Calculate drag for two interacting phases F(Drag,pq) • Calculate total drag for phase q , F(Drag,q) = F(Drag,q) + F(Drag,pq) , Sp2 in equation 4.26 4. Return total drag F(Drag,q) Table 4.3: Computing inter-component drag force by maxwellStefan() function The source terms added creates an explicit dependance of velocity on total pressure of the system. It is not a cause for concern since the drag term in our case is linearly related to relative velocity unlike the usual non-linear dependance. The solution divergence is controlled using lesser value for time steps to avoid instabilities.

39

4. Numerical Simulation Lessening the time steps may act as an under-relaxation factor for the iterative method. The presence of phase fraction ∇αq poses another restriction. Phaseintensive formulation of momentum equations suggested by Weller[7] is not used and hence cases involving sudden reductions of phase fractions to zero must be avoided which may drive the momentum equation to singularity or cause the presence of exceptionally large numbers due to high gradients of phase fractions resulting in floating point exceptions.

4.3.2

Source term: Porous drag

Momentum source term due to porous media derived in section 3.2 is rewritten here,

Fq,por = −ǫ∇p − gp τ

2

Dq k K + pq µ

−1

~υq,por = Sp3 + Sp4 .

(4.28)

Before proceeding further, three fields are introduced in multicomponentPorousFoam ; porosity(ǫ), tortuosity(τ ) and porousPlug(gp ). The three fields facilitate switching on the required terms when porous media is used in the solver and turning off porous drag effects when an open solid is used. All three fields are dimensionless and must be present in the starting conditions folder in the operating case file. The fields mentioned are set by using setFields application in the flow domain similar to the setting up of phase fractions and must adhere to the following conditions,

porosity =

tortuosity =

(

Vf (Vf +Vs )

1

in porous region, in non-porous region,

tortuosity in porous region, 1 in non-porous region,

porousPlug =

1 in porous region, 0 in non-porous region,

(4.29)

(4.30)

(4.31)

In order to introduce porosity effects in the equations, the momentum equation 4.25 is multiplied by porosity (ǫ) along with source terms Sp1 and Sp2 added due to drag. To accommodate path elongations in porous media, the source term Sp2 is multiplied by tortuosity(τ ). The porous drag term Sp4 is added to the momentum equation. The generic momentum equation that accomodates both porous and non-porous effects is then:

ǫ(

∂αq ρq ~υq + ∇ · (αq ρq ~υq ~υq ) = −αq ∇p + ∇ · T + αq ρq~g + ∂t Sp4 Sp1 + τ Sp2 + gp ) ǫ 40

(4.32)

4. Numerical Simulation All terms are implemented into equation 4.33 and can represent the desired effects. It can be verified that, in case of a non-porous region, according to 4.29 , 4.31 and 4.31, the equation represents a flow in non-porous region and the multiplicity of ǫ is also marginalized as it becomes 1. Sp4 is modelled implicitly in the solver as:

−

g p

ǫ

∗ fvm :: Sp(pow(tortuosity, 2) ∗ phase.porousDragMultiplier(), Uq ) . (4.33)

porousDragMultiplier() is a function that returns the coefficient of porous drag force of the phase. Coefficient or multiplier is the drag force term without the velocity of phase since the porous drag force term is considered implicitly, velocity is not known and has to be calculated for the current time. It is programmed as a phase property since it does not depend on many fluid properties which was in case of maxwellStefan() function. It is invoked as shown in 4.33 when setting up the momentum equation for all the phases q. The porous drag contribution of all the phases are calculated after solving the energy equation since the knudsen diffusion coefficient (Dq k ) is temperature dependant. The calculations are done in the file porousSource.H and stored for all the phases at current time. The computing procedure is shown in table 4.4.

1. Begin porousSource.H 2. Loop over all phases q. • Calculate knudsen diffusion coefficient Dq k from 3.24 3. Calculate permeability K from 3.26 4. Loop over all phases q • Calculate the porous drag force multiplier Fmq,por from 3.31 without considering ~υq • Store Fmq,por as q.porousDragMultiplier() 5. Exit porousSource.H Table 4.4: Computing porous drag multiplier by porousSource.H Thus the terms required for a generic momemtum equation which can switch roles to represent a porous or a non-porous zone is implemented into multicomponentPorousFoam .

4.3.3

Energy transport

The implementation of equation 3.45 is done in Tfluid.H file and is as follows:

41


( fvm::ddt(Coefficient1,Tfluid)+fvm::div(Coefficient2,Tfluid) == fvm::laplacian(Coefficient3,Tfluid)+porosity*fvc::ddt(p));

(4.34)

The variables and functions needed by the above equation 4.35 is listed in table 4.5. Variable

Field

Calculation form

Cp kappa Coefficient1 Coefficient2 Coefficient3

volScalarField volScalarField volScalarField volScalarField surfaceScalarField

Pn αq C q Pq=1 n k=1 αq kq ρm Cp+(1−ǫ)(ρpor −Cpor ) ǫk Pn+ (1 − ǫ)kpor q=1 ρq φq Cq

k : kappa, Thermal conductivity subscript por: porous material

Invoking function fluid.Cp() fluid.kappa() — fluid.Coefficient2() —

ρm : Density of the mixture

Table 4.5: Variables and functions required for energy transport The algorithm to solve for temperature of the mixture is shown in table 4.6 and is as follows,

1. Calculate thermal conductivity kappa for the mixture 2. Calculate specific heat capacity of the mixture 3. Calculate Coefficient1 4. Calculate Coefficient2 by invoking the function Coefficient2() 5. Calculate Coefficient3 6. Set up the temperature transport equation for the mixture 4.35 7. Relax the temperature equation 8. Solve the temperature equation 9. Update viscosity of individual components 3.46 Table 4.6: Solution procedure for temperature of the mixture Only transport of temperature is solved in multicomponentPorousFoam . The densities of individual phases can be updated after solving the temperature equation but is out of the scope of this work.

4.3.4

Solution algorithm

The solution process of multicomponentPorousFoam is listed in table 4.7. Solving the temperature of the mixture results in updating temperature dependant properties in the flow. A condition at the start of setting up a case is that the porous

42

4. Numerical Simulation related fields have to be initialized in order to ditinguish between porous and nonporous zones which initiates appropriate consideration of source terms in the momentum equation as discussed earlier. The remaining steps involved can be seen in the solution algorithm.

a Set fields for phase fractions, porosity, porousPlug and tortuosity b runtime++ 1 Start Pimple loop 2 Update mixtue turbulence 3 Solve for phase fractions using fluid.solve() 4.4 4 Update mixture density 5 Solve for mixture temperature 4.25 – Update viscosity of individual components 3.46 – Update Knudsen diffusion coefficient of all components 3.24 – Update diffusion coefficients of two interacting components 3.47 6 Calculate porous drag multipliers 4.4 7 Set momentum equations for all components 4.33 8 Pressure corrector loop – Solve for pressure – Correct flux – Reconstruct velocities from corrected flux 9 Update Substantive derivative for all components 10 End Pimple loop c End time loop Table 4.7: multicomponentPorousFoam algorithm

4.4

Closure

A brief introduction and the essence of numerical simulation to solve fluid flow problems were discussed. Reasons for working with OpenFOAM were then stated. Characteristics of the existing solver multiphaseEulerFoam which solves multiphase flows by the Eulerian approach were then analyzed. The solution method used by the solver to solve phase fractions were then discussed and explained with appropriate schematics and equations. The MULES algorithm used to limit and solve a convective only transport system of a phase was also explained. Further, handling of pressure-velocity coupling by flux prediction and pressure correction was discussed and demonstrated with the help of discritized equations. The semi-implicit treatment of drag force by the solver was explained.

43

4. Numerical Simulation Development of multicomponentPorousFoam, an extended and modified solver for flow in porous media was then implemented. The modifications to the momentum equation and development of a generic equation for the same which represents both porous and non-porous zones with the help of certain fields; porosity, porousPlug and tortuosity is briefly explained and implemented. Volume-averaged energy equation to solve for temperature transport of the mixture is developed which is also generic and accommodates porous and non-porous zones. Solution Algorithm of multicomponentPorousFoam was then enumerated 4.7 .

44

5. Validation

5. Validation 5.1

Case1: Loschmidt tube

Validation of our solver for diffusion dominated flows involving mass transfer is undertaken in this section. Such flow scenarios are quite often validated by using diffusion tubes [27, 25, 23]. We use the Loschmidt tube for validation. Duncan and Toor conducted a few experiments to examine diffusion in an ideal ternary gas mixture containing Hydrogen, Nitrogen and Carbon-di-oxide [24]. It basically consists of two bulbs connected by a capillary as shown in figure 5.1. Among the three gases, one gas is made stationary and two gases are set in motion. This is achieved by adding almost equal quantities of one gas in both the bulbs (q2 ) resulting in a very small or in fact a negligible gradient of the gas rendering it almost stationary. The remaining two gases are added one in each bulb resulting in large gradients thus driving the two gases in motion when the stop cock is removed in order to achieve homogeniety.In our simulation, Loschmidt tube with Methane (CH4 ), Argon (Ar) and Hydrogen (H2 ) is simulated.

q1+ q2 q3-

q1 q2 q3+

Bulb 1

Bulb 2

Figure 5.1: Experimetal set up of loschmidt tube In such a scenario, the gas with a small gradient (q2 ) behaves in a certain way. Toor predicted the behaviour and named them; Osmotic diffusion, Reverse diffusion and Diffusion barrier [10]. It is important to elaborate these terms which occur in a ternary mixture.

45

5. Validation 1. Osmotic diffusion: At time t = a , though the gas has negligible driving force, it moves. 2. Reverse Diffusion: During a time interval a < t < b, the gas moves in a direction opposite to its existing trifle gradient thus defying Fick’s law which states that, gases diffuse in a direction normal to their concentration gradients. 3. Diffusion barrier: Between b < t < c, the gas does not diffuse even though a large gradient exists. The diffusion flux is zero despite a large driving force [24] The experimental results have been compared with the Maxwell-Stefan equations by Duncan and Toor [12] and is in good agreement. Hence the analytical solution of the Maxwell-Stefan equation is obtained by assuming a 1−D flow and verified with the results generated by the solver.

5.1.1

Case Setup

In order to simulate the Loschmidt tube, the whole process has to be abstracted to a simple model. The first abstraction done is the consideration of only the capillary and eliminating the two bulbs. Since the whole system is closed and entites of interest are the phase fractions, it is justified. The geometry used for simulation is shown in figure 5.2

Wall

Left tube

-

-

Figure 5.2: Simulation model of loschmidt tube The length of the capillary is reduced to 100µm. High gradients prevail initially at the division of left tube and right tube. In order to resolve them, edge lengths of 10−6 are required [27]. To save computation time, the length of the whole capillary is reduced. The experiment is conducted at atmospheric pressure of 101.3kP a and under isothermal conditions. Composition of the three gases in both left and right tube is given in table 5.1 It is observed that Ar behaves as gas q2 explained earlier and must undergo the stated osmotic diffusion, reverse diffusion and diffusion barrier effects. All the boundary patches are considered as walls ensuring a closed surface. The internal field is set to 101.3kP a. The pressure at the walls are set to zero gradient condition. Zero gradient means that there is no change of a quantity in normal direction. Velocity

46

5. Validation Component CH4 H2 Ar

Left tube 0.295 0.4 0.305

Right tube 0.405 0.3 0.295

Table 5.1: Initial composition

at the wall is set to slip condition for all components. The slip condition ensures that the tangential component of velocity is zero gradient as explained earlier and the normal velocity component near the wall is zero [5]. An analytical solution to the problem assuming a 1−D transient problem is developed in [25] (pages 110 − 115). A sci-lab program written based on the solution method which solves the phase fraction of a three component mixture from its initial composition for any time t in the left tube can be found in appendix A.1. A graph is generated for each component using the program to plot the variation of its phase fraction with time.

5.1.2

Results 0.31 b c

0.309 0.308 0.307 0.306

calculated alphaAr analytical alphaAr

0.305 a 0.304 0.303 0.302 0.301 0.3

0

2e-05

4e-05

6e-05

8e-05 0.0001 0.00012 0.00014

Figure 5.3: Phase fraction ofArgon in left tube: Analytical vs. Calculated a − b : Reverse diffusion. b − c :Diffusion barrier. > c :Normal diffusion Figure 5.3 shows the behaviour of Argon. From time a until time b, reverse diffusion takes place. Though the left tube has a higher concentration, there is an incoming flux to the left tube(0.305) from right tube(0.295) and hence an increase and decrease of phase fraction on left and right respectively. Between time b and c, a diffusion barrier exists resulting in stagnation of flux even though a considerable gradient exists. After c, Argon diffuses normally. This assures the physical reality of solution. The other two component gases begin to start homogenizing following the Fick’s law of diffusion and diffuse normally, down their concentration gradients. The animation of the intial stage of the other all the gases can be seen in appendix A.2

47

5. Validation 0.35 0.34 0.33 calculated alphaCH4 analytical alphaCH4

0.32 0.31 0.3 0.29

0

2e-05

4e-05

6e-05

8e-05 0.0001 0.00012 0.00014

Figure 5.4: Phase fraction of CH4 in left tube : Analytical vs. Calculated Methane behaves normally and obeys Fick’s law of diffusion. The initial composition of Methane in the left tube was comparatively lower. Hence phase fraction gradually increases until homogenization as shown in figure 5.4. 0.4 0.395 0.39 0.385 0.38 calculated alphaH2 analytical alphaH2

0.375 0.37 0.365 0.36 0.355 0.35

0

2e-05

4e-05

6e-05

8e-05 0.0001 0.00012 0.00014

Figure 5.5: Phase fraction of H2 in left tube: Analytical vs. Calculated There is not much of a surprise in case of Hydrogen too. It diffuses normally down its gradient. The initial composition of Hydrogen being high in left tube in comparison to right tube reduces in order to effect the homogenization process as seen in figure 5.5. The end time shown in the graph is not the time required to completely homogenize. There is no further change in the behaviour till the homogenization of both tubes.

48

5. Validation From figures 5.3, 5.4 and 5.5, it is observed that analytical results and the result of phase fraction of all three components calculated from the solver is in good agreement. Thus the solver is validated for diffusion dominated flows involving mass transfer.

5.2

Case2: Flow through a channel

Validation of volume averaged energy transport and porous drag source terms is based on calculation of heat transfer coefficients in thermally fully developed laminar flows. The dimensionless Nusselt number is finally calculated from the velocity temperature profile when fully developed. There is no influence of mass flow rate, gas properties and absolute temperature on Nusselt number when the temperature profile is fully developed. It then depends exclusively on the flow profile [31]. Hence it serves as a good criteria for validation since, in case of the energy equation and porous media. Nusselt number for a flow in a rectangular duct with infinite length in z-direction and with a height H is given by, ∂(Tfluid)

∂n

Tfluid

v

H y

x

Tw

Figure 5.6: Calculation of Nusselt number

2H Nu = Tw − T¯

∂T ∂n

(5.1) w

T¯ is the mean temperature over the cross section of the duct. Based on the velocity profile, it can be determined by [6],

T¯ =

RH

(υT )dy RH υdy 0

0

(5.2)

Figure 5.6 is a schematic depicting the determination of the required values of temperature and velocity of the fluid at a fixed location on the x-axis. The variation of the values along the y-direction and the normal gradient of temperature at the wall must be considered as shown in the figure. Thus the Nusselt number is calculated using equation 5.2 and equation 5.1.

49

5. Validation

5.2.1

Case Setup

A duct which is partially cooled by walls as shown in figure 5.7 is used for validation. Two zones 1 and 2 exist in the geometry. Zone 1 is surrounded by adiabatic walls and zone 2 by cooling walls. A homogeneous mixture of Oxygen and Nitrogen is then fed into the channel. The component fractions of Oxygen and Nitrogen are maintained at 0.21 and 0.79 respectively. The mixture enters through the inlet at 373.15K. The walls in zone 2 are set to 323.15K. The starting field of the two components are also set homogeneous. Tw=323.15

Adiabatic Walls

Tfluid=373.15K 2

1

0.04m 0.12m

Figure 5.7: Case set up for a channel flow The mass flow specified at inlet is 4 × 10−3 kgs−1 to generate a laminar flow profile. It results in a parabolic flow profile at the outlet. The velocities of individual components are then calculated proportional to their phase fraction and set in the case directory. Two variations for zone 2; non-porous and porous zones is set. For the non-porous case, the porosity is 1. In the second case, zone 2 is filled by a porous region with a porosity of 0.9. Further, the density of porous media is taken as 5900kgm−3 , specific heat capacity of 500 Jkg−1 K. Tortuosity factor of the porous media is 1.1. Diameter of pores is assumed to be 100µm. A symmetry plane is used for the case specified making the top plane as the plane of symmetry. When a steady state solution is sought, a solution symmetric to the plane of symmetry exists [14]. Since the fully developed temperature profile is the solution that we seek, the use of symmetry is justified. Dirichlet boundary condition [5] is used for pressure at the outlet and is 0. Velocity of both components are given at the inlet from the known flux. A no-slip condition for velocity is used at the walls. For temperature of the mixture (Tfluid ), Neumann boundary condition [5] with gradient as zero is set at outlet. At inlet, a temperature of 373.13K is fixed and for the walls in zone 2, 323.15K is prescribed. For adiabatic walls in zone 1, a zero gradient boundary condition is used. In order to find the Nusselt number from equation 5.1, the velocity and temperture profile in y-direction is probed at a location of x= 0.1m. A program in Sci-lab is used (appendix B) to calculate the Nusselt number using the equations 5.2 and 5.1. A comparison of the Nusselt number generated by the program and the literature values of Nusselt number for the two kinds of profiles generated gives the proximity of calculations by the solver to accurate values.

50

5. Validation

5.2.2

Validation of temperature transport

Velocity profile developed when zone 2 is non porous is shown in figure 5.8. The flow is fully developed and is as expected. The temperature profile since it is transported by flux develops a similar profile of velocity. There is no loss of heat in the adiabatic zone 1 and a boundary layer is developed in zone 2 where a temperature gradient exists.

Figure 5.8: Velocity profile of a non porous zone 2 Variation of temperature when zone 2 is non porous is shown in figure 5.9. Variation of velocity and temperature along y-axis at x= 0.1m is shown in figure 5.10. The average velocity developed is 0.16714 ms−1 . In order to calculate the gradient normal to the wall, a point on the y−axis has to be selected. Since the gradient of temperature normal to wall is required, the value of temperature sampled next to the wall at a y coordinate 4.20168 × 10−5 m is chosen. The temperature at this point is 323.661. The value of temperature gradient normal to the wall is thus 12161.802Km−1 . The mean temperature is 350.633 K. Thus the calculated Nusselt number is 8.85. The literature value of Nusselt number for a parabolic flow profile is 7.54 [30]. A relative error is observed between the values and the reasons are discussed subsequently.

5.2.3

Validation of porosity effects

The velocity profile for a porous zone is shown in figure 5.11. There is a marginal increase in maximum velocity. Temperature contour is shown in figure 5.12. It is observed that there is a change in the contour developed due to the presence of porous media. As the temperature

51

5. Validation

Figure 5.9: Temperature profile of a non porous zone 2 0.25

360 355

0.2

350 345

0.15

340

Tfluid

U 0.1

335 330

0.05

325 320

0

0.001

0.002

0.003

0.004

0.005

0

0

0.001

0.002

0.003

0.004

0.005

Figure 5.10: Temperature and velocity variation (non porous zone 2) at x= 0.1m gradient reduces with positive x-direction, the lesser it diffuses thus generating the contour shown. Variation of velocity and temperature along y-axis at x= 0.1m is shown in 5.10. The reduced values of temperature can be observed. The temperature profile follows the velocity profile. The average velocity developed in case of a porous zone 2 is 0.16718 ms− 1. Temperature is sampled at the same point considered for the previous case. The value of temperature here is 323.223K. This is because of the block profile case of a porous zone with a porosity value of 0.9. The value of temperature gradient normal to the wall is 1737.4003Km−1 . The mean temperature is 326.95619 K. Nusselt number for this case is 9.13. The literature value of Nusselt number for a block flow profile is 9.80 [31].

52

5. Validation

Figure 5.11: Velocity profile for a non porous zone 2

Figure 5.12: Temperature profile for a porous zone 2 328

0.25

327.5 327

0.2

326.5 326

0.15

325.5

Tfluid

U

325

0.1

324.5 324

0.05

323.5 323

0

0.001

0.002

0.003

0.004

0.005

0

0

0.001

0.002

0.003

0.004

0.005

Figure 5.13: Temperature and velocity variation (porous zone 2) at x= 0.1m

5.3

Closure

The Loschmidt tube considered as the validation case for Maxwell-Stefan source terms was set up and explained. Abstraction of experimental set up for CFD simu-

53

5. Validation lations were also shown. The behaviour patterns of a ternary mixture was explained. Since it is a closed geometry, the pressure value used here was 101.325KPa. The large value of pressure significantly reduced the time step to obtain a stable simulation. The time step observed for the simulation was 10−10 seconds. The analytical solution developed as a Scilab program for a transient 1 dimensional case was used to verify the results. The phenomenon observed for comparison was the variation of phase fraction of all three gases in the left tube of the capillary. By observing the variation of phase fraction, it is indeed equivalent to the observation of variation of flux of each phase. The simulation results are in very good agreement with the analytical results over the complete observation period. Additionally, the various kinds of diffusion behavior were also captured in the simulation. The second case of fully developed flow over parallel plates was used for the validation of temperature transport and porous drag effects. The case setup involves a sudden cooling of a mixture entering into a region cooled by the walls. Nusselt number was used as the validation criteria and compared with values available in the literature. In order to validate both the stated effects, the cooling zone was considered nonporous for the first case and as porous zone with a porosity of 0.9 for the second. For non-porous region, the Nusselt number calculated was 8.85 compared to the literature value of 7.54. For the porous case, the Nusselt number calculated was 9.13 compared to the literature value of 9.8. Since the solver deals with incompressible phases, instabilities were observed when total pressure was used due to the presence of large numbers which may increase cumulative errors. Hence only the relative pressure or pressure fluctuation values were used for this case. The terms related by ideal gas law were kept intact by adding a value 105 Pa to the pressure terms. This may be the cause for the observed deviation. The temperature profile developed for the two cases are realistic and rational.

54

6. Conclusions

6. Conclusions Scope of this work and its implications have been extensively discussed in the previous chapters. A few important conclusion and further potential improvements are stated here: • A circumspective attitude is generally taken towards the Eulerian approach to solve for transport of multicomponent gases. The reason being the complexity and computation costs expended for the approach. For flow scenarios which require slip velocities of individual components, the mixture approach which is the most chosen alternative to Eulerian approach cannot be used since it circumvents solving individual component velocities. Hence the Eulerian approach is chosen and models the mulitcomponent flow scenario accurately. • Maxwell stefan relations models diffusion dominated problems involving mass transfer. In multicomponent flows, inter-component and not intra-component momentum tranfer takes place which is established in section 3.1.2. The phenomena of osmotic diffusion, reverse diffusion and diffusion barrier is accurately modelled by Maxwell-Stefan relations which the Fick’s law fails to predict. • The intercomponent momentum transfer implemented to the modified solver multicomponentporousFoam simulates ternary diffusion in a Loschmidt tube accurately. The results are in perfect agreement with the analytical solution for the case solved as transient one dimensional flow. The graph of variation of phase fractions can be seen in figures 5.4 5.5 and 5.3. • Phase intensive momentum equations are not used in the solver. Hence flows involving sudden changes of phase fraction of a component must be avoided. The presence of the gradient term of phase fraction is responsible for this short coming which cause floating point exceptions when it encounters large gradients. The use of upwind interpolation scheme to solve for phase fractions restricts the use of large time steps. It is not a major problem for transient flows due to the presence of a small time step. When a steady state solution is sought, the demand for a small time step may be a drawback.

55

6. Conclusions • Volume averaged equations are implemented in the solver. The presence of porososity, porousPlug and tortuosity fileds provides a generic operating form for the solver. The porous drag terms and the effect of porosity can be switched on in a particular zone by setting the mentioned fields appropriately. Currently, isotropic and non deformable porous media have been modelled. Extensions to orthotropic media can be done by adding appropriate tensorial terms in the solver. The equations are inversed for this purpose. The results generated by the solver is in close agreement. • Transport of temperature of the mixture assuming local thermodynamic equilibrium between components accurately predicts the transfer of temperature in the pseudo homogeneous media. The difference between the contour of temperature in a non porous domain and porous domains can be seen in figures 5.9 and 5.12. • The deviation in Nusselt number is attributed to the use of only the relative pressure values for simulation. But it does not largely affect the simulation since all the terms which are related by ideal gas properties are properly accounted by adding a value of 105 Pa to the relative pressure solved for incompressible flows. The calculation of gradients at the wall in order calculate Nusselt number is also speculated to be the cause for the error. The probe location for calculation of gradient at a wall can vary continuously in the cell adjacent to the wall. The temperature also varies accordingly resulting in many possibilities to calculate temperature gradient at the wall.And even a small variation in its value affects the final solution quite largely. • Temperature dependance of viscosity and diffusion coefficients have been modelled in the solver. Density is considered as a scalar in the solver. This restricts its dependance on temperature. It is also predicted to be a source of error in the final solution. A large dependance of variables exist on density in the solver. Hence to convert it into a scalarField may prove detrimental since the aspect of considering individual components as incompressible may be violated and cause mass conservation problems. • Source terms which consists of pressure terms are modelled explicitly. Projecting pressure dependant source terms out of the momentum equation and solving for them in the poisson equation which is set for pressure improves the accuracy and may even reduce the small time step requirement. But the presence of velocity terms simultaneously poses a problem in doing so since pressure-velocity is solved using a segregated approach.Thus, using a coupled solution procedure would be suitable for handling the source terms which contain pressure and velocity terms. The article [16] which implements a coupled pressure based solution using a block matrix structure may be taken as a cue for establishing a coupled solution procedure. Further, Maxwell Stefan terms are added explicitly in the solver. A semi-implicit treatment of the terms may improve accuracy in flows involving large pressure gradients. • Reaction terms can be easily added to the phase continuity equation in the solver. The MULES algorithm which solves the equation contains explicit and implicit source terms as its arguments. Thus reaction can be modelled by passing the source term field to the function.

56

A. Loschmidt tube

A. Loschmidt tube A.1

Scilab Program to calculate 1-D analytical solution

function ternaryDiffusion () x1m = i n p u t ( ” Enter x1 i n l e f t h a l f ”) ; // x1 : x1p = i n p u t ( ” Enter x1 i n r i g h t h a l f ”) ; x2m = i n p u t ( ” Enter x2 i n l e f t h a l f ”) ; // x2 : x2p = i n p u t ( ” Enter x2 i n r i g h t h a l f ”) ; x3m = i n p u t ( ” Enter x3 i n l e f t h a l f ”) ; // x3 : x3p = i n p u t ( ” Enter x3 i n r i g h t h a l f ”) ; l = i n p u t ( ” Enter t h e l e n g t h o f t h e d i f f u s i o n d12 = i n p u t ( ” D i f f u s i o n 1−2 ”) ; d13 = i n p u t ( ” D i f f u s i o n 1−3 ”) ; d23 = i n p u t ( ” D i f f u s i o n 2−3 ”) ; time = i n p u t ( ” I n t e r e s t e d probe time ”) ; de l t a T=i n p u t ( ” d e l t a t ”) ; x1=(x1m+x1p ) / 2 ; x2=(x2m+x2p ) / 2 ; x3=(x3m+x3p ) / 2 ; S=x1∗ d23+x2∗ d13+x3∗ d12 ; D11=d13 ∗ ( ( x1∗ d23+(1−x1 )∗ d12 ) ) / S ; D12=(x1 ∗ d23 ∗( d13−d12 ) ) / S ; D21=(x2 ∗ d13 ∗( d23−d12 ) ) / S ; D22=(d23 ∗( x2∗ d13+(1−x2 )∗ d12 ) ) / S ; D=[D11 D12 , D21 D22 ] ; data = [ ] ;

57

Methane Argon Hydrogen tube ”) ;

A. Loschmidt tube

// Eigen Values a =1; b=−(D11+D22 ) ; c=D11∗D22−D21∗D12 ; D1=(−b+s q r t ( bˆ2−4∗a∗ c ) ) / 2 ; D2=(−b−s q r t ( bˆ2−4∗a∗ c ) ) / 2 ; d i s p (D1 , D2 ) ; // Modal matrix P p11 =1; p12 =1; p21=(D1−D11 ) / D12 ; p22=(D21 ) / ( D2−D22 ) ; P=[ p11 p12 ; p21 p22 ] ; d i s p (P ) ; Pinv=i n v (P ) ; r =60; t=0 p i l 2 =60; f 1 =(1/2) −((4/(% p i ˆ 2 ) ) . . ∗ exp(−(% p i ˆ2∗D1∗ t ) / ( 4 ∗ l ˆ 2 ) ) ) ; f 2 =(1/2) −((4/(% p i ˆ 2 ) ) . . ∗ exp(−(% p i ˆ2∗D2∗ t ) / ( 4 ∗ l ˆ 2 ) ) ) ; //End f =[ f 1 0 ; 0 f 2 ] ; // d i s p ( f ) ; x D i f f =[x1p−x1m ; x2p−x2m ] ; P i n v x D i f f=Pinv ∗ x D i f f ; f P i n v x D i f f=f ∗ P i n v x D i f f ; P f P i n v x D i f f=P∗ f P i n v x D i f f ; x i n i t i a l =[x1m ; x2m ] ; x= P f P i n v x D i f f+ x i n i t i a l ; xnormal=x ( 2 ) ; // F u n c t i o n s f o r t =0: de l t a T : time p i l 2 =60; f1 secondTerm = 0 ; f2 secondTerm = 0 ; f o r k =0:1:2

58

A. Loschmidt tube

m=k + 0 . 5 ; f1 secondTerm = f1 secondTerm . . + ( 1 / (mˆ 2 ) ) . . ∗ exp (−(mˆ2∗% p i ˆ2∗D1∗ t ) / ( l ˆ 2 ) ) ; f2 secondTerm = f2 secondTerm . . + ( 1 / (mˆ 2 ) ) . . ∗ exp (−(mˆ2∗% p i ˆ2∗D2∗ t ) / ( l ˆ 2 ) ) ; end f 1 =(1/2)−(1/% p i ˆ2)∗ f1 secondTerm ; f 2 =(1/2)−(1/% p i ˆ2)∗ f2 secondTerm ;

//End f =[ f 1 0 ; 0 f 2 ] ; // d i s p ( f ) ; x D i f f =[x1p−x1m ; x2p−x2m ] ; P i n v x D i f f=Pinv ∗ x D i f f ; f P i n v x D i f f=f ∗ P i n v x D i f f ; P f P i n v x D i f f=P∗ f P i n v x D i f f ; x i n i t i a l =[x1m ; x2m ] ; x= P f P i n v x D i f f+ x i n i t i a l ; // p l o t ( ( ( t ∗ d12 ) / ( l ˆ 2 ) ) , ( x ( 2 ) / xnormal ) , ’ ∗ ’ ) ; plot (t , x (2) , ’∗ ’); // data ( s i z e ( data , ’ r ’ ) + 1 , 1 ) = ( t ∗ d12 ) / ( l ˆ 2 ) ; // data ( s i z e ( data , ’ r ’ ) , 2 ) = x ( 2 ) / xnormal ; data ( s i z e ( data , ’ r ’)+1 ,1)= t ; data ( s i z e ( data , ’ r ’ ) , 2 ) = x ( 1 ) ; data ( s i z e ( data , ’ r ’ ) , 3 ) = x ( 2 ) ; data ( s i z e ( data , ’ r ’) ,4)=1 − x(1) −x ( 2 ) ;

end

d i s p ( data ) ; f p r i n t f M a t ( ’ data . dat ’ , data ) ; endfunction

59

A. Loschmidt tube

A.2

Animation of initial diffusion of gases

t=0.5 x 10-6 s

t=2.5 x 10-6 s

t=4 x 10-6 s

t=6 x 10-6 s

t=9.5 x 10-6 s

Figure A.1: Initial diffusion of Argon

60

A. Loschmidt tube

t=0.5 x 10-6 s

t=2.5 x 10-6 s

t=4 x 10-6 s

t=6 x 10-6 s

t=9.5 x 10-6 s

Figure A.2: Initial diffusion of Methane

61

A. Loschmidt tube

t=0.5 x 10-6 s

t=2.5 x 10-6 s

t=4 x 10-6 s

t=6 x 10-6 s

t=9.5 x 10-6 s

Figure A.3: Initial diffusion of Hydrogen

62

B. Scilab program to calculate Nusselt number

B. Scilab program to calculate Nusselt number a1 = f s c a n f M a t ( ” f i l e n a m e . xy ”) ; m=s i z e ( a1 ) ; sum ut =0; sum u =0; f o r i = 1 : 1 :m( 1 ) sum ut=sum ut+a1 ( i , 2 ) ∗ a1 ( i , 3 ) ; sum u=sum u+a1 ( i , 3 ) ; end d i s p ( sum ut ) ; T mean=sum ut /sum u ; d i s p ( T mean ) ; gradT=(a1 (2 ,2) − a1 ( 1 , 2 ) ) / ( a1 (2 ,1) − a1 ( 1 , 1 ) ) ; d i s p ( gradT ) ; Nu=(0.02∗ gradT ) / ( T mean − 3 2 3 . 1 5 ) ; d i s p (Nu ) ; d i s p ( a1 ( 1 , 2 ) , a1 ( 1 , 1 ) ) ; The above code reads a .xy file where the values of temperature and magnitude of velocity is stored as columns. The file can be generated by sampleDict utility in OpenFOAM where sampling for Tfluid and magU can be set for a particular time. Values of magU can be calculated by the utility foamCalc. Sampling must be done for the time when the temperature profile is fully developed. The code is specific for the case explained in section 5.2.1. In case of change in geometry or boundary conditions, height of the parallel plates and temperature at the wall in zone 2 must be changed appropriately.

63

BIBLIOGRAPHY

Bibliography [1] ANSYS FLUENT Theory Guide 12.0. [2] ANSYS FLUENT User’s Guide 12.0. [3] R.C.Reid et al. The properties of gases and liquids, 4 ed. McGraw-Hill, New York, 1988. [4] T.C.Bond et.al. Bounding the role of black Carbon in the climate system: A scientific assessment. American Geophysical Union, doi:10.1002/jgrd.50171, 2013. [5] OpenFOAM Programmer’s guide 2.1.1. [6] K.Stephan H.D.Baehr. Heat and mass transfer , first ed. Springer, 1998. [7] H.G.Weller. Derivation, modelling and solution of the conditionally averaged two-phase flow equations. [8] H.G.Weller. Alternative pressure-velocity algorithm. Private Communication, 1999. [9] H.G.Weller and C. Fureby G.Tabor, H.Jasak. A tensorial approach to computational continuum mechanics using object-oriented techniques. American Institute of Physics, 1998. [10] H.L.Toor. Diffusion in three component gas mixtures. A.I.Ch.E.J. 3, 198-207, 1957. [11] Hrvoje Jasak. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows. PhD thesis, Imperial college, University of London, 1996. [12] J.B.Duncan and H.L.Toor. An experimental study of three component gas diffusion. A.I.Ch.E.J. 3, 198-207, 1957. [13] J.C.Slattery. Momentum, Energy and Mass Transfer in continua. McGraw-Hill, New York, 1972. [14] Milovan Peric Joel H.Ferziger. Computational Methods for Fluid Dynamics, Second Edition. Springer, 1999. [15] C.L.Tien K. Vafai. Boundary and inertia effects on flow and heat transfer in porous media. int. J. Heat Mass Transfer 24, 1981.

64

BIBLIOGRAPHY [16] Hrvoje Jasak Steffen Schutz Karsten Urban Manfred Piesche Kathrin Kissling, Julia Springer. A Coupled pressure based solution algorithm based on the volume of fluid approach for two or more immicible fluids. V European conference on Computational Fluid Dynamics, ECCOMAS CFD, 2010. [17] H.A Kramers and J.Kistemaker. On the slip of a diffusing gas mixture along a wall. [18] V. Taivassalo M. Manninen and S. Kallio. On the mixture model for multiphase flow. VTT Publications 288,Technical Research Centre of Finland, 1996. [19] M.Kaviany. Laminar flow through a porous channel bounded by isothermal parallel plates. int. J. Heat Mass Transfer 28 (4), 1985. [20] M.Kaviany. Principles of Heat Transfer in Porous Media. Mechanical Engineering Series, Springer, 1991. [21] Suhas V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, 1980. [22] Edwin N. Lightfoot R.Byron Bird, Warren E. Stewart. Transport Phenomena, Second Edition. [23] R.Krishna. Problems and pitfalls in the use of Fick formulation for intraparticle diffusion. Chem. Eng. Sci. 48(5),(845-861), 1993. [24] J.A.Wesselingh R.Krishna. The Maxwell-Stefan approach to mass transfer. Elsevier Science Ltd. [25] R.Krishna Ross Taylor. MULTICOMPONENT MASS TRANSFER. Elsevier Science Ltd, 1996. [26] Henrik Rusche. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. [27] Manfred Piesche Stephan Goll. Multi-component gas transport in micro-porous domains:Mutltidimensional simulation at the macroscale. International Journal of Heat and Mass Transfer. [28] S.Whitaker. Advances in theory of fluid motion in porous media. Ind.Eng Chem. 61(12), 1969. [29] S.Whitaker. Role of species the species momentum equation in the analysis of the Stefan diffusion tube. Ind. Eng. Chem. Fund.30, 1991. [30] VDI-GesellschaftVerfahrenstchnik und Chemieingenieurwesen (Ed.). VDI Heat Atlas. second ed., Springer, 2010. [31] H.C.Perkins W.M.Kays. Forced Convection, internal flows in ducts, in J.P.Hartnett, W.M.Rohsenow(Eds.), Handbook of Heat transfer. McGraw-Hill, New York, p.124, 1973.

65

CFD Simulation of Multicomponent Gas Flow Through Porous Media

Recommend Documents