- Email: [email protected]
Multiphase flow in a capillary porous medium
Computational Materials Science 28 (2003) 704–713 www.elsevier.com/locate/commatsci
Multiphase flow in a capillary porous medium q T. Ricken *, R. de Boer Institute of Mechanics, University of Duisburg-Essen, Location Essen, 45117 Essen, Germany
Abstract Based on the theory of porous media, a calculation concept for the multiphase flow in a capillary porous medium will be presented. We will exclusively investigate the rise of liquids in porous bodies due to the capillarity phenomenon. The field equations used consist of the mechanical balance equations and the physical constraint conditions. The treatment of the capillary problem, based on thermomechanical investigations, yields the result that the capillarity force is a volume interaction force and depends on the free Helmholtz energy functions of the phases and the density gradient of the liquid. With respect to the aforementioned outcome, further constitutive relations for the ternary model are developed. The aim of this investigation is the numerical simulation of the behavior of liquid and gas phases in a rigid porous body at rest. Therefore, the needed weak formulations of the governing field equations, i.e. the balance equations of mass and the balance equations of momentum of the liquid and gas phases, will be given. The usefulness of the proposed theory will be demonstrated with an example. Ó 2003 Elsevier B.V. All rights reserved.
1. Introduction In the petroleum industry, soil mechanics, and building physics, the multiphase flow in capillary porous media is a point of high interest. On the one hand, one is interested in overcoming the partially strong capillary force in order, for example, to extract oil from reserves, while on the other hand one is interested in following the
q
A part of this manuscript has been published in Acta Mechanica 159 (2002) 173 and was submitted as a contribution to the Proceedings IUTAM Symposium on the Mechanics of Physicochemical and Electromechanical Interactions in Porous Media, Eindhoven University of Technology, The Netherlands, 18–23 May 2003, Kluwer Academic Publishers, 2003. * Corresponding author. Tel.: +49-201-183-2679; fax: +49201-183-2680. E-mail address: [email protected] (T. Ricken).
transport of liquid in capillary porous bodies. Due to the irregular distribution and different sizes and shapes of the pores a rigorous averaging is necessary. This is provided by the porous media theory allowing a continuum thermomechanical approach to the two-phase flow problem in a porous solid. The instationary and stationary motion of liquids in porous solids with small pores is complex and not all related problems have been solved. The main internal forces acting in these systems are the capillary and friction forces which are both influenced by the saturation, though unfortunately not in the same way. The capillary forces have been recognized as intermolecular forces which are created by cohesion and adhesion at interfaces. From the evaluation of the thermomechanical treatment we could clearly make out that these forces depend on the
0927-0256/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2003.08.032
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
free Helmholtz energy functions of the solid and gas phases and the density gradient of the liquid. The friction forces are taken into account using the well-known van Genuchten formulation. Both forces are introduced as interaction forces between the phases involved. Using these constitutive equations, the equations of motion of the liquid and gas phases are developed which form the framework for the presented ternary calculation concept. 2. The basic relations of the theory of porous media In the next section the volume fraction concept will be introduced and other basic relations will be taken from the mixture theory, with some modifications. 2.1. Volume fraction concept In the volume fraction concept, it is assumed that the porous solid always models a control space and that only the liquids and/or gases contained in the pores can leave the control space. Furthermore, it is assumed that the pores are statistically, homogeneously distributed and that an arbitrary volume element in the reference and the actual placement is composed of the volume elements of the real constituents. A porous medium occupying the control space of the porous solid BS , with the boundary oBS , consists of constituents ua , with real volumes va , where the index a denotes j individual constituents. The boundary oBS is a material surface for the solid phase and a non-material surface for the liquid and gas phases. The volume fraction concept contains the volume fractions na , which refer the volume elements dva of the individual constituents ua to the bulk volume element dv (the volume fraction concept can be founded by an average theory on the microscale, see [5]): j X dva na ðx; tÞ ¼ ; na ðx; tÞ ¼ 1; ð1Þ dv a¼1 where x is the position vector to the spatial point x in the actual placement and t is the time. The
705
volume fractions na in (1)1 satisfy the volume fraction condition (1)2 for j constituents ua . Moreover, the partial density of the constituent ua , namely qa , is related to the real density of the materials qaR involved via the volume fractions na : qa ¼ na qaR :
ð2Þ
The concept of volume fractions is an important part of the theory following in the next sections. Due to the volume fraction concept, all geometric and physical quantities, such as motion, deformation, and stress, are defined in the total control space, and, thus, they can be interpreted as the statistical average values of the real quantities. 2.2. Kinematics The ternary model, a rigid porous solid (a ¼ S) filled with a mixture of an incompressible liquid (a ¼ L) and a compressible gas (a ¼ G), will be treated as an immiscible mixture of all three constituents ua with particles Xa , where at any time t each spatial point of the current placement of the solid phase uS is simultaneously occupied by particles XS , XL and XG of the solid, liquid and gas phases. These particles proceed from different reference positions Xa at time t ¼ t0 . Thus, each constituent is assigned its own independent motion function: x ¼ va ðXa ; tÞ;
Xa ¼ v1 a ðx; tÞ:
ð3Þ
Eq. (3)1 represents the Lagrange description of motion. The function va is postulated to be unique, and uniquely invertible, at any time t. The existence of a function inverse to (3)1 leads to the Euler description of motion, see (3)2 . A mathematically necessary and sufficient condition for the existence of (3)2 is given, if the Jacobian Ja ¼ detFa differs from zero. Therein, Fa is the deformation gradient. Fa and its inverse F1 a are defined as Fa ¼ grada va ;
F1 a ¼ grad Xa :
ð4Þ
The differential operator ‘‘grada ’’ denotes a partial differentiation with respect to the reference position Xa of the constituent ua and the differential operator ‘‘grad’’ referring to the spatial point x. During the deformation process, Fa is restricted to detFa > 0.
706
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
With the Lagrange description of the motion (3)1 , the velocity and acceleration fields of the constituents ua are defined as material time derivatives of the motion function (3)1 : x0a
ov ðXa ; tÞ ; ¼ a ot
x00a
o2 va ðXa ; tÞ ¼ : ot2
ð5Þ
For scalar fields depending on x and t, the material time derivatives are defined as: 0
ð Þa ¼
oð Þ þ ½gradð Þ x0a : ot
ð6Þ
Readers interested in a review of the kinematics of porous media are referred to [5].
3. Balance equations As aforementioned, the investigated porous body consists of uS (solid), uL (liquid) and uG (gas). Due to the rigidity of the skeleton of the body, only the balance equations of the liquid and gas phases will be taken into account. Corresponding to that, the index b ¼ L, G only denotes the mobile liquid and gas phases in the following contribution. The given balance equations will be formulated in local forms for each individual constituent taking all interaction effects into consideration. The system is investigated under the condition of any mass exchange between the constituents and the restriction of an incompressible liquid (qLR ¼ const:) and a compressible gas. Under these assumptions, the local statements of the balance equations of mass are given for the constituents uL and uG by 0
ðnL ÞL þ nL div x0L ¼ 0; 0
0
nG ðqGR ÞG þ ðnG ÞG qGR þ qG div x0G ¼ 0
ð7Þ
volume force and ^pb describes the interaction of the constituents ub and are restricted to ^pL þ ^pG ¼ ^pS :
The balance of moment of momentum for nonpolar materials states that the material time derivative of the moment of momentum is equal to the moments of all external forces, where the moments are referred to a fixed point O. The evaluation of this balance equation yields CauchyÕs second equation of motion, namely the T symmetry of the stress tensor Tb ¼ ðTb Þ .
4. Constitutive modelling From the evaluation of the entropy inequality of the saturated porous body considering the saturation condition as a constraint, see [5], we obtain the constitutive relations for Tb : Tb ¼ nb kI þ TbE ¼ ðnb k þ pEb ÞI ¼ nb pbR I ¼ pb I;
wherein cG ¼ qGR 0 TRG denotes a parameter concerning the reference density qGR 0 , the temperature T and the specific gas constant RG . Moreover, we gain the definition of the interaction force ^pb with ^pb ¼ k grad nb þ ^pbE ; wherein
pL ¼ 0; div TL þ qL ðb x00L Þ þ ^
4.1. The partial Cauchy stress tensor ð8Þ
In (8)1;2 , ‘‘div’’ denotes the divergence operator, Tb is the partial Cauchy stress tensor, qb b is the
ð10Þ
with the interface pressure k, the partial extra (effective) stress TbE , and the partial extra (effective) pressure pEb of the mobile phases ub . The realistic pressure pbR cannot be determined for the incompressible liquid phase but in the case of the compressible gas phase, pGR is related to the real density qGR with the non-linear expression: GR q0 pGR ¼ cG ln ; ð11Þ qGR
and the local statements of the balance equations of momentum are given by
div TG þ nG qGR ðb x00G Þ þ ^ pG ¼ 0:
ð9Þ
^pbE
ð12Þ
denotes the extra supply term.
In a porous body saturated by two or more phases it is difficult to say how the pore pressure acts on the phases. In order to describe the pressure behavior of the mobile phases ub , the
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
so-called ‘‘carrier model’’ has been developed in [6]. In this model, the extra pressure pEb depends on the surface tension between the phases. Therefore, one has to distinguish between two different cases. In the first case, we look at a porous body which is saturated with liquid for the most part. As long as a dynamic balance exists in the system, the pore pressure and the pressure in the bubble differs in the extra pressure pEG . In this state, the liquid phase is called the carrier phase. On the contrary, if the porous body is filled mostly with gas, the value of the pressure difference in the bubbles becomes equal to the extra pressure pEL and the gas phase is called the carrier phase. With the expression k ¼ nS k þ nL pLR þ nG pGR (DaltonÕs law), the constrain pEL ¼ pEG has to be considered. Depending on the carrier phase, the formulations of the extra pressures pEb in (10) have been developed as follows, see [7]: Liquid as carrier: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pEG ¼ nL nG ^k G lnðnG =ðnF nW ÞÞ and ð13Þ pEL ¼ pEG ; Gas as carrier: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pEL ¼ nL nG ^k L lnðnL =nW Þ and pEG ¼ pEL :
ð14Þ
In (13) and (14), ^k L and ^k G denote material parameters depending on the surface tension between the constituents ub , nF ¼ nL þ nG describes the whole pore volume (porosity) and nW determines the point of change concerning the carrier phase. It should be mentioned that in the case of viscous fluids the partial extra stress TbE is additionally determined by motion. Moreover, the viscosity of the mobile phases influences the interaction forces ^ pbE , as a viscous friction. With regard to the consider partial fluid bodies the viscosity inside the fluids itself can be neglected. With the help of (13) and (14) in connection with the definition of the partial Cauchy stress tensor (10), we constitutively create a relation between the volume fraction nb and the realistic pressure of the phases pbR . In the case of the liquid as the carrier phase the function for nL reads as:
L
W
F
n ¼ ðn n Þ exp
ðpLR pGR Þ2 2 ð^k G nF Þ
707
! þ nF : ð15Þ
If the gas phase is the carrier, nL is determined by: ! 2 ðpLR pGR Þ L W n ¼ ðn Þ exp : ð16Þ ð^k L nF Þ2 Both relations have got a line of intersection at nW . The relations (15) and (16) are similar to the well-known capillary pressure–saturation relation by Brooks and Corey, see e.g. [4]. In this approach, the realistic liquid and gas pressures are coupled by the capillary pressure pC through pC ¼ pGR pLR . Indeed, this approach yields to a similar description but contains two important differences, this approach is not motivated by an thermomechanical investigation but on a phenomenological description, also, this approach does not contain a separate description for the pressure action between the mobile phases against each other on the one hand and the interaction forces between the mobile phases and the solid body (capillarity) on the other hand as given in this contribution. In fact, these two phenomena are put together in one constitutive relation using a capillary pressure– saturation relation. 4.2. Interaction forces The only interaction forces are the interaction forces ^pS , ^pL and ^pG , aforementioned. In an extended thermodynamic investigation, de Boer and Didwania [6] have developed constitutive equations for the interaction forces ^pL and ^pG , which read as: ! ^S ^G L LR L S ow G ow ^p ¼ p grad n q q grad qL oqL oqL SL x0L S G x0G ; ^pG ¼ pGR grad nG qG S L x0L SG x0G ;
ð17Þ ow^G grad qL oqL ð18Þ
708
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
where w^S and w^G denotes the free Helmholtz energy functions of the solid and gas phase and Sb and S b are friction functions. 4.2.1. Capillarity The capillary interaction forces ^tSL and ^tLG , with ^S ^tSL ¼ qS ow grad qL ¼ ^tSL grad qL ; oqL
ð19Þ
^G ^tLG ¼ qG ow grad qL ¼ ^tLG grad qL ; oqL
ð20Þ
grad qL ¼
oqL : oxi
ð22Þ
We consider only the capillary force in e3 direction (driving volume force). There are also capillary forces in both the e2 - and the e1 -directions. In particular, the force in e2 -direction is responsible for an increase in the static and dynamic friction, if the wall is rough. We assume here a smooth wall of the narrow tube and neglect the influence of the capillary forces in e2 - and e1 directions. Therefore, (22) simplifies to grad qL ¼ e3 ;
acting on the surface of the constituents. At the interactions liquid–solid and liquid–gas there are jumps in the densities and the existing thin films are under tensions, see e.g. [3]. Unfortunately, the parameters ^tSL and ^tLG are not directly measurable so that we develop a first approach in order to identify these quantities. Therefore, we investigate the rising meniscus in a thin tube, see Fig. 1. For the one-dimensional transport of the liquid in e3 direction with the coordinate x3 , the volume force is equal to zero, because grad qL vanishes due to the fact that qL is constant. Thus, the capillary driving volume force only exists in the thin films between the wall of the tube and the interface meniscus of the liquid and gas. Considering only the liquid phase with respect to (17), (19) and (20) the capillary driving volume force reads as: t ¼ ð^tSL þ ^tLG Þgrad qL :
ð21Þ
The gradient of the liquid density qL can be calculated in the base-system introduced in Fig. 1:
ð23Þ L
where we have made use of the fact that q is equal to minus or plus unity depending on the density ratio between the density of the liquid and the density of the other phase involved. Thus, considering (21) the capillary driving volume force t is given by t ¼ ð^tSL ^tLG Þe3 ;
ð24Þ
wherein we considered that the capillary forces in e3 -direction have opposite directions. We calculate the driving force from (24) by integrating over the thickness h of the thin film at the wall (see Fig. 1): Z sþ Z sþ SL ^t e3 ds ^tLG e3 ds: t¼ ð25Þ s
s
In e3 -direction at the wall the membrane forces are assumed to be constant. This assumption is in accordance with the assumption in the shell theory for membrane forces. Thus t ¼ ðtSL tLG Þe3 with tSL ¼ ^tSL h; tLG ¼ ^tLG h;
h ¼ sþ s
ð26Þ
r
e1
s+
ρG
e3 a3
is obtained. In general, the driving volume forces tSL ¼ tSL e3 and tLG ¼ tLG e3 are not accessible and measurable. However, the norms of the forces can be replaced equivalently by:
detail
t e2
s-
h
ϑ ρL
a1
meniscus
tS
x3 = s
t = tSL- t LG tS = tSa1
wall
Fig. 1. Force analysis at the meniscus.
tSL tLG ¼ tS cos #:
ð27Þ
It should be mentioned that tSL , tLG and tS are the norms of the forces per unit length in e3 direction and a1 -direction, which act only at the boundary of the meniscus and not over the whole
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
cross-section. In order to determine the total external force at the boundary of the meniscus we compare (27) with the Young equation which appoints the norm of the driving force with: t ¼ r cos #2pr;
ð28Þ
wherein r denotes the surface tension, # is the contact angle between the meniscus and the tube and r represents the radius of the tube. Motivated by this, for the average quantities ^tSL and ^tLG we proclaim that ^tSL ¼ nS rSL d cos #; ^tLG ¼ nG rLG d cos #:
ð29Þ
Herein, rSL and rSL with rSL rLG ¼ r denotes the surface tension between the solid and liquid phase and the liquid and gas phase, respectively. The material parameter d denotes a average value depending on the form and nature of the pores in the porous body. This new approach to the interaction forces allows for the capillary motion in porous solids. 4.2.2. Relative permeability The last terms on the right side of (17) and (18) denote the friction influence on each of the fluids in relation to their seepage velocities. Thereby, SL and SG concern the interaction between the solid phase and the constituents ub and they are restricted to a positive value. These material parameter functions are postulated in dependence on the volume fraction nL with well-known relations by van Genuchten [1], which read in a modified way as: m 2 SL1 ¼ k0L Se1=2 1 1 Se1=m ; ð30Þ 1 G 1=2 1=2 1=m 2m 1 Se : SG ¼ k 0 1 Se 2
In (30)1;2 , k0b ¼ kFb =½ðnL Þ cbR represents a specific soil parameter related to the Darcy permeability coefficient kFb and the specific weight cbR . Moreover, Se ¼ nL =nF denotes the saturation and m is a van Genuchten parameter. Note that SL and SG only describe the relative permeability concerning to the mutational volume fractions. Do not mistake these for the capillary pressure–saturation relation mentioned before.
709
The parameter S b in (30)1;2 characterizes the friction influence of liquid and gas velocities on each other. However, with regard to the investigated fluids these friction rates will be neglected in the following contribution. With these formulations, the set of constitutive equations is complete and we proceed to the numerical modelling.
5. Numerical treatment Within the framework of a standard Galerkin procedure, weak formulations are formed. Therefore, one can either choose the saturation–density formulation or the pressure–pressure formulation, depending on the problem being solved. In the saturation–density formulation we gain the realistic density of the gas qGR , the volume fraction of the liquid nL and the velocities x0b of ub as primary variables and in the pressure–pressure formulation we receive the pressures of both phases pLR and pGR as unknowns. The most important difference between the two formulations can be found in the numerical cost. Because of the reduction from 6 to 2 unknowns the pressure–pressure formulation is very effective. However, this has to be paid for by neglecting dynamic effects. Moreover, the numerical benchmarks like stability and convergence are different concerning the two formulations. Currently, investigations to this points are treated by the authors. In the following sections, we will present the weak formulations of the two versions. The structure of the linearized weak formulations will be given according to the general form: _ Dh_ ¼ R þ RoB ; DG Dh þ DG
ð31Þ
wherein DG multiplied by the vector of the increase of process variables, Dh represents the share _ multiplied by the vector of time of stiffness, DG derivation of increase of the process variables, Dh_ denotes the vector of is the share of damping, R known quantities and RoB represents the vector concerning boundary conditions on the surface oB of the porous body B.
710
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
2
5.1. Saturation–density formulation In the saturation–density formulation we attain the set of unknowns RSD with RSD ¼ RSD ðx; tÞ ¼ fnL ; qGR ; x0L ; x0G g:
ð32Þ
In order to determine these quantities, the weak formulations of the balance of mass and balance of momentum of the liquid and gas phases can be used. They read as follows: Balance of mass for the liquid phase multiplied by the weight function dnL : Z 0 ML ¼ fðnL ÞS þ nL div x0L gdnL dv ¼ 0: ð33Þ B
Balance of mass for the gas phase multiplied by the weight function dnG : Z 0 0 MG ¼ fnG ðqGR ÞG þ ðnG ÞG qGR B
þ nG qGR div x0G gdqGR dv ¼ 0:
ML 6 MG DG ¼ 6 4 IL IG
ð34Þ
Balance of momentum for the liquid phase multiplied by the weight function dx0L : Z IL ¼ fðnL pLR ÞI grad dx0L þ qLR nL ðb x00L Þdx0L B
þ pLR grad nL dx0L qLR ð^tSL þ ^tLG Þgrad nL dx0L Z SL x0L dx0L g dv ¼ fnL pLR n dx0L g da: oBL
ð35Þ Balance of momentum for the gas phase multiplied by the weight function dx0G : Z IG ¼ fðnG pGR ÞI grad dx0G þ qGR nG ðb x00G Þdx0G
2
ML 6 MG _ ¼6 DG 4 0 0
0 MG IL IG
ML 0 IL IG
0 MG 0 0
0 0 IL 0
3 32 0 dnL 6 GR 7 MG 7 76 dq 0 7; 5 IL 4 dxL 5 IG dx0G 32 3 0 dnL 6 GR 7 0 7 76 dq 0 7; 5 0 4 dxL 5 IG dx0G
ð38Þ
2
RoB
3 0 6 0 7 7 ¼6 4 RIL 5: RIG
ð37Þ
ð39Þ
In (39), the boundary conditions RIL ¼ nL pLR n and RIG ¼ nG pGR n denote the partial pressure on the surface of the porous body and result from the surface share in (35) and (36). 5.2. Pressure–pressure formulation The second possibility to reformulate the governing set of equations is the pressure–pressure formulation. Therefore, we modify the set of unknowns to RPP ¼ RPP ðx; tÞ ¼ fpLR ; pGR g:
ð40Þ
In order to reduce the set of equations, the velocities of the constituents x0L and x0G have to be explicitly determined using the balance of momentum in the Darcy formulation: x0L ¼ SL1 fnL grad pLR þ qL b ð^tSL þ ^tLG Þgrad qL g; ð41Þ
B
þ pGR grad nG dx0G qLR^tLG grad nL dx0G Z SG x0G dx0G g dv ¼ fnG pGR n dx0G g da: oBG
ð36Þ In Eqs. (35) and (36), Neumann boundary conditions are formulated on the surface oBL and oBG. Herein, n is the outward oriented unit surface normal on the respective boundary. _ and the vectors Dh The matrix of DG and DG and Dh_ as well as the vector concerning the boundary conditions RoB are given as:
x0G ¼ SG1 fnG grad pGR þ qG b ^tLG grad qL g:
ð42Þ
It must be noted that in this formulation one has to neglect the acceleration x00b . After some modifications the weak formulations of the balances of mass with the weight function dpLR and dpGR reads as: Z 0 fðnL ÞS dpLR grad dpLR nL x0L g dv ML ¼ B Z ¼ fnL x0L dpLR ng da; ð43Þ oBL
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
MG ¼ ¼
Z
6. Numerical example
0
fðqGR nG ÞS dpGR grad dpGR qG x0G g dv
ZB
fqGR nG x0G dpGR ng da:
The linearized weak equations of the saturation–density and pressure–pressure formulation are implemented in the well-known finite element program FEAP by Taylor. The finite element discretization was realized with the Galerkin method. In the example, see Fig. 2(a), we simulate the moistening of a brick standing with one of the small sides in a water reservoir. The left side is closed by an isolation for the liquid and gas phase but the other sides are open. Fig. 2(b) shows the finite element discretization with 403 isoparametric rectangular elements. Due to the expected high saturation gradient at the lower right corner of the brick, the mesh becomes refined at this area. In an experiment, K€ unzel [2] measured the distribution of the volume fraction after 72 days, see Fig. 2(c). The solution of the numerical simulation are depicted in the following. In Fig. 3, the calculated distribution of the liquid versus time (color chart) and the vectors of the capillary force are given. The result of the calculation after 72 days in Fig. 3(c) shows a good approach to the measurement.
ð44Þ
oBG
Thus, for the linearized weak formulation we gain: LR ML ML dp DG ¼ ; ð45Þ MG MG dpGR _ ¼ DG
R
oB
ML MG
ML MG
dpLR ; dpGR
ð46Þ
RML ¼ : RMG
ð47Þ
For the Neumann boundary conditions we have got on the Neumann boundaries oBL a liquid volume efflux RML ¼ nL x0L n and on oBG a gas mass efflux RMG ¼ qGR nG x0G n. In the following example, we use the pressure– pressure formulation for the simulation with linear shape functions concerning the fluid pressures pLR and pGR .
modell
711
measurement
discretization 50
[cm]
49 0.05
impermeable 0.4
40
0.3
30
Height[cm]
open
g 0.2
0.10
0.15
20
0.1
10
0
0
volume fraction nL
0.20
0.25
0
suction
(a)
0
0.1 [cm]
(b) Fig. 2. Saturation of a brick.
5
11 10
Width [cm]
(c)
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
0.5
0.5
0.5 L
L
n [-] 0.255874 0.233806 0.211737 0.199324 0.18829 0.177256 0.164842 0.153808 0.142774 0.13174 0.120706 0.109672 0.0986374 0.0872236 0.0751898 0.0600179 0.0549228 0.053676 0.0521647
y
0.3
0.2
0.1
0.4
0.255874 0.233806 0.211737 0.199324 0.18829 0.177256 0.164842 0.153808 0.142774 0.13174 0.120706 0.109672 0.0986374 0.0872236 0.0751898 0.0600179 0.0549228 0.053676 0.0521647
0.3
y
0.4
L
n [-]
0.2
0.1
0
n [-] 0.4
0.3
0.2
0.1
0
0 0
0.2
x
0.255874 0.233806 0.211737 0.199324 0.18829 0.177256 0.164842 0.153808 0.142774 0.13174 0.120706 0.109672 0.0986374 0.0872236 0.0751898 0.0600179 0.0549228 0.053676 0.0521647
y
712
0
(a)
0
0.2
x
0.2
x
(b)
(c)
Fig. 3. Liquid saturation and capillarity force versus time. Volume fraction nL after 1 min (a), after 8 days (b), and after 72 days (c).
0.5
0.5
L
L
p [Pa] 0.4 2.58346 2.5819 2.56332 2.47346 2.29679 2.12011 1.94343 1.76676 1.59008 1.41341 1.23673 1.06006 0.883379 0.706704 0.530028 0.353352 0.176676 0
y
0.3
0.2
0.1
0.3
0.2
0.1
0
0
0
0.2
x
p [Pa] 0.4
2.58346 2.5819 2.56332 2.47346 2.29679 2.12011 1.94343 1.76676 1.59008 1.41341 1.23673 1.06006 0.883379 0.706704 0.530028 0.353352 0.176676 0
y
0.4
L
p [Pa]
2.58346 2.5819 2.56332 2.47346 2.29679 2.12011 1.94343 1.76676 1.59008 1.41341 1.23673 1.06006 0.883379 0.706704 0.530028 0.353352 0.176676 0
0.3
y
0.5
0.2
0.1
0 0
0.2
x
(a)
0
0.2
x
(b)
(c)
Fig. 4. Liquid pressure and liquid velocity direction versus time. Liquid pressure pL after 1 min (a), after 8 days (b), and after 72 days (c).
The distribution of the liquid pressure (color chart) and the velocity direction of the liquid phase (unit vectors) are presented in Fig. 4. One can recognize that the moving saturation front causes
a pressure swelling in the liquid phase, see Fig. 4(b) and (c). In the case of dynamic effects, like the quick capillary saturation of tissues with saturation
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
times of several seconds, instead of the pressure– pressure formulation the saturation–pressure formulation has to be used. Furthermore, the authors are interested in wave propagations in multiphase saturated porous materials. Therefore, only the saturation–pressure formulation can be applied. Readers who are interested in an example with dynamic effects are referred to [7].
7. Conclusion A closed calculation concept for the simulation of multiphase flow in capillary active porous media has been developed. After a brief introduction to the theory of porous media, the used field equations, the constitutive equations and the numerical treatment have been presented. Herein, an essential contribution is given by implementation of the capillary attraction into the interaction forces on the basis of an thermodynamical investigation. To this subject, there have been many empirical approaches but only a limited number of continuum mechanical treatments are known. The presented contribution by the authors in the framework of a ‘‘smeared’’ continua has led to
713
reasonable field equations, in particular, to proper constitutive equations for the interaction volume forces which cause the capillary flow. References [1] M.T. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Science Society of America Journal 44 (1980) 892–898. [2] H.M. K€ unzel, Verfahren zur ein- und zweidimensionalen Berechnung des gekoppelten W€arme- und Feuchtetransports in Bauteilen mit einfachen Kennwerten, Dissertation, Lehrstuhl Konstruktive Bauphysik, Fakult€at Bauingenieurund Vermessungswesen der Universit€at Stuttgart, 1994. [3] H.T. Davis, Statistical Mechanics of Phases, Interfaces and Thin Films, VCH Publishers, New York, Weinheim, Cambridge, 1996. [4] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems, Springer-Verlag, Berlin, Heidelberg, New York, 1997. [5] R. de Boer, Theory of Porous Media: Highlights in the Historical Development and Current State, Springer-Verlag, Berlin, Heidelberg, New York, 2000. [6] R. de Boer, A.K. Didwania, Capillarity in porous bodies: contribution of the Vienna school and recent findings, Acta Mechanica 159 (2002) 173–188. [7] T. Ricken, Kapillarit€at in por€ osen medien––theoretische untersuchung und numerische simulation, Dissertation, Shaker Verlag, Aachen, 2002.
Multiphase flow in a capillary porous medium q T. Ricken *, R. de Boer Institute of Mechanics, University of Duisburg-Essen, Location Essen, 45117 Essen, Germany
Abstract Based on the theory of porous media, a calculation concept for the multiphase flow in a capillary porous medium will be presented. We will exclusively investigate the rise of liquids in porous bodies due to the capillarity phenomenon. The field equations used consist of the mechanical balance equations and the physical constraint conditions. The treatment of the capillary problem, based on thermomechanical investigations, yields the result that the capillarity force is a volume interaction force and depends on the free Helmholtz energy functions of the phases and the density gradient of the liquid. With respect to the aforementioned outcome, further constitutive relations for the ternary model are developed. The aim of this investigation is the numerical simulation of the behavior of liquid and gas phases in a rigid porous body at rest. Therefore, the needed weak formulations of the governing field equations, i.e. the balance equations of mass and the balance equations of momentum of the liquid and gas phases, will be given. The usefulness of the proposed theory will be demonstrated with an example. Ó 2003 Elsevier B.V. All rights reserved.
1. Introduction In the petroleum industry, soil mechanics, and building physics, the multiphase flow in capillary porous media is a point of high interest. On the one hand, one is interested in overcoming the partially strong capillary force in order, for example, to extract oil from reserves, while on the other hand one is interested in following the
q
A part of this manuscript has been published in Acta Mechanica 159 (2002) 173 and was submitted as a contribution to the Proceedings IUTAM Symposium on the Mechanics of Physicochemical and Electromechanical Interactions in Porous Media, Eindhoven University of Technology, The Netherlands, 18–23 May 2003, Kluwer Academic Publishers, 2003. * Corresponding author. Tel.: +49-201-183-2679; fax: +49201-183-2680. E-mail address: [email protected] (T. Ricken).
transport of liquid in capillary porous bodies. Due to the irregular distribution and different sizes and shapes of the pores a rigorous averaging is necessary. This is provided by the porous media theory allowing a continuum thermomechanical approach to the two-phase flow problem in a porous solid. The instationary and stationary motion of liquids in porous solids with small pores is complex and not all related problems have been solved. The main internal forces acting in these systems are the capillary and friction forces which are both influenced by the saturation, though unfortunately not in the same way. The capillary forces have been recognized as intermolecular forces which are created by cohesion and adhesion at interfaces. From the evaluation of the thermomechanical treatment we could clearly make out that these forces depend on the
0927-0256/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2003.08.032
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
free Helmholtz energy functions of the solid and gas phases and the density gradient of the liquid. The friction forces are taken into account using the well-known van Genuchten formulation. Both forces are introduced as interaction forces between the phases involved. Using these constitutive equations, the equations of motion of the liquid and gas phases are developed which form the framework for the presented ternary calculation concept. 2. The basic relations of the theory of porous media In the next section the volume fraction concept will be introduced and other basic relations will be taken from the mixture theory, with some modifications. 2.1. Volume fraction concept In the volume fraction concept, it is assumed that the porous solid always models a control space and that only the liquids and/or gases contained in the pores can leave the control space. Furthermore, it is assumed that the pores are statistically, homogeneously distributed and that an arbitrary volume element in the reference and the actual placement is composed of the volume elements of the real constituents. A porous medium occupying the control space of the porous solid BS , with the boundary oBS , consists of constituents ua , with real volumes va , where the index a denotes j individual constituents. The boundary oBS is a material surface for the solid phase and a non-material surface for the liquid and gas phases. The volume fraction concept contains the volume fractions na , which refer the volume elements dva of the individual constituents ua to the bulk volume element dv (the volume fraction concept can be founded by an average theory on the microscale, see [5]): j X dva na ðx; tÞ ¼ ; na ðx; tÞ ¼ 1; ð1Þ dv a¼1 where x is the position vector to the spatial point x in the actual placement and t is the time. The
705
volume fractions na in (1)1 satisfy the volume fraction condition (1)2 for j constituents ua . Moreover, the partial density of the constituent ua , namely qa , is related to the real density of the materials qaR involved via the volume fractions na : qa ¼ na qaR :
ð2Þ
The concept of volume fractions is an important part of the theory following in the next sections. Due to the volume fraction concept, all geometric and physical quantities, such as motion, deformation, and stress, are defined in the total control space, and, thus, they can be interpreted as the statistical average values of the real quantities. 2.2. Kinematics The ternary model, a rigid porous solid (a ¼ S) filled with a mixture of an incompressible liquid (a ¼ L) and a compressible gas (a ¼ G), will be treated as an immiscible mixture of all three constituents ua with particles Xa , where at any time t each spatial point of the current placement of the solid phase uS is simultaneously occupied by particles XS , XL and XG of the solid, liquid and gas phases. These particles proceed from different reference positions Xa at time t ¼ t0 . Thus, each constituent is assigned its own independent motion function: x ¼ va ðXa ; tÞ;
Xa ¼ v1 a ðx; tÞ:
ð3Þ
Eq. (3)1 represents the Lagrange description of motion. The function va is postulated to be unique, and uniquely invertible, at any time t. The existence of a function inverse to (3)1 leads to the Euler description of motion, see (3)2 . A mathematically necessary and sufficient condition for the existence of (3)2 is given, if the Jacobian Ja ¼ detFa differs from zero. Therein, Fa is the deformation gradient. Fa and its inverse F1 a are defined as Fa ¼ grada va ;
F1 a ¼ grad Xa :
ð4Þ
The differential operator ‘‘grada ’’ denotes a partial differentiation with respect to the reference position Xa of the constituent ua and the differential operator ‘‘grad’’ referring to the spatial point x. During the deformation process, Fa is restricted to detFa > 0.
706
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
With the Lagrange description of the motion (3)1 , the velocity and acceleration fields of the constituents ua are defined as material time derivatives of the motion function (3)1 : x0a
ov ðXa ; tÞ ; ¼ a ot
x00a
o2 va ðXa ; tÞ ¼ : ot2
ð5Þ
For scalar fields depending on x and t, the material time derivatives are defined as: 0
ð Þa ¼
oð Þ þ ½gradð Þ x0a : ot
ð6Þ
Readers interested in a review of the kinematics of porous media are referred to [5].
3. Balance equations As aforementioned, the investigated porous body consists of uS (solid), uL (liquid) and uG (gas). Due to the rigidity of the skeleton of the body, only the balance equations of the liquid and gas phases will be taken into account. Corresponding to that, the index b ¼ L, G only denotes the mobile liquid and gas phases in the following contribution. The given balance equations will be formulated in local forms for each individual constituent taking all interaction effects into consideration. The system is investigated under the condition of any mass exchange between the constituents and the restriction of an incompressible liquid (qLR ¼ const:) and a compressible gas. Under these assumptions, the local statements of the balance equations of mass are given for the constituents uL and uG by 0
ðnL ÞL þ nL div x0L ¼ 0; 0
0
nG ðqGR ÞG þ ðnG ÞG qGR þ qG div x0G ¼ 0
ð7Þ
volume force and ^pb describes the interaction of the constituents ub and are restricted to ^pL þ ^pG ¼ ^pS :
The balance of moment of momentum for nonpolar materials states that the material time derivative of the moment of momentum is equal to the moments of all external forces, where the moments are referred to a fixed point O. The evaluation of this balance equation yields CauchyÕs second equation of motion, namely the T symmetry of the stress tensor Tb ¼ ðTb Þ .
4. Constitutive modelling From the evaluation of the entropy inequality of the saturated porous body considering the saturation condition as a constraint, see [5], we obtain the constitutive relations for Tb : Tb ¼ nb kI þ TbE ¼ ðnb k þ pEb ÞI ¼ nb pbR I ¼ pb I;
wherein cG ¼ qGR 0 TRG denotes a parameter concerning the reference density qGR 0 , the temperature T and the specific gas constant RG . Moreover, we gain the definition of the interaction force ^pb with ^pb ¼ k grad nb þ ^pbE ; wherein
pL ¼ 0; div TL þ qL ðb x00L Þ þ ^
4.1. The partial Cauchy stress tensor ð8Þ
In (8)1;2 , ‘‘div’’ denotes the divergence operator, Tb is the partial Cauchy stress tensor, qb b is the
ð10Þ
with the interface pressure k, the partial extra (effective) stress TbE , and the partial extra (effective) pressure pEb of the mobile phases ub . The realistic pressure pbR cannot be determined for the incompressible liquid phase but in the case of the compressible gas phase, pGR is related to the real density qGR with the non-linear expression: GR q0 pGR ¼ cG ln ; ð11Þ qGR
and the local statements of the balance equations of momentum are given by
div TG þ nG qGR ðb x00G Þ þ ^ pG ¼ 0:
ð9Þ
^pbE
ð12Þ
denotes the extra supply term.
In a porous body saturated by two or more phases it is difficult to say how the pore pressure acts on the phases. In order to describe the pressure behavior of the mobile phases ub , the
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
so-called ‘‘carrier model’’ has been developed in [6]. In this model, the extra pressure pEb depends on the surface tension between the phases. Therefore, one has to distinguish between two different cases. In the first case, we look at a porous body which is saturated with liquid for the most part. As long as a dynamic balance exists in the system, the pore pressure and the pressure in the bubble differs in the extra pressure pEG . In this state, the liquid phase is called the carrier phase. On the contrary, if the porous body is filled mostly with gas, the value of the pressure difference in the bubbles becomes equal to the extra pressure pEL and the gas phase is called the carrier phase. With the expression k ¼ nS k þ nL pLR þ nG pGR (DaltonÕs law), the constrain pEL ¼ pEG has to be considered. Depending on the carrier phase, the formulations of the extra pressures pEb in (10) have been developed as follows, see [7]: Liquid as carrier: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pEG ¼ nL nG ^k G lnðnG =ðnF nW ÞÞ and ð13Þ pEL ¼ pEG ; Gas as carrier: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pEL ¼ nL nG ^k L lnðnL =nW Þ and pEG ¼ pEL :
ð14Þ
In (13) and (14), ^k L and ^k G denote material parameters depending on the surface tension between the constituents ub , nF ¼ nL þ nG describes the whole pore volume (porosity) and nW determines the point of change concerning the carrier phase. It should be mentioned that in the case of viscous fluids the partial extra stress TbE is additionally determined by motion. Moreover, the viscosity of the mobile phases influences the interaction forces ^ pbE , as a viscous friction. With regard to the consider partial fluid bodies the viscosity inside the fluids itself can be neglected. With the help of (13) and (14) in connection with the definition of the partial Cauchy stress tensor (10), we constitutively create a relation between the volume fraction nb and the realistic pressure of the phases pbR . In the case of the liquid as the carrier phase the function for nL reads as:
L
W
F
n ¼ ðn n Þ exp
ðpLR pGR Þ2 2 ð^k G nF Þ
707
! þ nF : ð15Þ
If the gas phase is the carrier, nL is determined by: ! 2 ðpLR pGR Þ L W n ¼ ðn Þ exp : ð16Þ ð^k L nF Þ2 Both relations have got a line of intersection at nW . The relations (15) and (16) are similar to the well-known capillary pressure–saturation relation by Brooks and Corey, see e.g. [4]. In this approach, the realistic liquid and gas pressures are coupled by the capillary pressure pC through pC ¼ pGR pLR . Indeed, this approach yields to a similar description but contains two important differences, this approach is not motivated by an thermomechanical investigation but on a phenomenological description, also, this approach does not contain a separate description for the pressure action between the mobile phases against each other on the one hand and the interaction forces between the mobile phases and the solid body (capillarity) on the other hand as given in this contribution. In fact, these two phenomena are put together in one constitutive relation using a capillary pressure– saturation relation. 4.2. Interaction forces The only interaction forces are the interaction forces ^pS , ^pL and ^pG , aforementioned. In an extended thermodynamic investigation, de Boer and Didwania [6] have developed constitutive equations for the interaction forces ^pL and ^pG , which read as: ! ^S ^G L LR L S ow G ow ^p ¼ p grad n q q grad qL oqL oqL SL x0L S G x0G ; ^pG ¼ pGR grad nG qG S L x0L SG x0G ;
ð17Þ ow^G grad qL oqL ð18Þ
708
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
where w^S and w^G denotes the free Helmholtz energy functions of the solid and gas phase and Sb and S b are friction functions. 4.2.1. Capillarity The capillary interaction forces ^tSL and ^tLG , with ^S ^tSL ¼ qS ow grad qL ¼ ^tSL grad qL ; oqL
ð19Þ
^G ^tLG ¼ qG ow grad qL ¼ ^tLG grad qL ; oqL
ð20Þ
grad qL ¼
oqL : oxi
ð22Þ
We consider only the capillary force in e3 direction (driving volume force). There are also capillary forces in both the e2 - and the e1 -directions. In particular, the force in e2 -direction is responsible for an increase in the static and dynamic friction, if the wall is rough. We assume here a smooth wall of the narrow tube and neglect the influence of the capillary forces in e2 - and e1 directions. Therefore, (22) simplifies to grad qL ¼ e3 ;
acting on the surface of the constituents. At the interactions liquid–solid and liquid–gas there are jumps in the densities and the existing thin films are under tensions, see e.g. [3]. Unfortunately, the parameters ^tSL and ^tLG are not directly measurable so that we develop a first approach in order to identify these quantities. Therefore, we investigate the rising meniscus in a thin tube, see Fig. 1. For the one-dimensional transport of the liquid in e3 direction with the coordinate x3 , the volume force is equal to zero, because grad qL vanishes due to the fact that qL is constant. Thus, the capillary driving volume force only exists in the thin films between the wall of the tube and the interface meniscus of the liquid and gas. Considering only the liquid phase with respect to (17), (19) and (20) the capillary driving volume force reads as: t ¼ ð^tSL þ ^tLG Þgrad qL :
ð21Þ
The gradient of the liquid density qL can be calculated in the base-system introduced in Fig. 1:
ð23Þ L
where we have made use of the fact that q is equal to minus or plus unity depending on the density ratio between the density of the liquid and the density of the other phase involved. Thus, considering (21) the capillary driving volume force t is given by t ¼ ð^tSL ^tLG Þe3 ;
ð24Þ
wherein we considered that the capillary forces in e3 -direction have opposite directions. We calculate the driving force from (24) by integrating over the thickness h of the thin film at the wall (see Fig. 1): Z sþ Z sþ SL ^t e3 ds ^tLG e3 ds: t¼ ð25Þ s
s
In e3 -direction at the wall the membrane forces are assumed to be constant. This assumption is in accordance with the assumption in the shell theory for membrane forces. Thus t ¼ ðtSL tLG Þe3 with tSL ¼ ^tSL h; tLG ¼ ^tLG h;
h ¼ sþ s
ð26Þ
r
e1
s+
ρG
e3 a3
is obtained. In general, the driving volume forces tSL ¼ tSL e3 and tLG ¼ tLG e3 are not accessible and measurable. However, the norms of the forces can be replaced equivalently by:
detail
t e2
s-
h
ϑ ρL
a1
meniscus
tS
x3 = s
t = tSL- t LG tS = tSa1
wall
Fig. 1. Force analysis at the meniscus.
tSL tLG ¼ tS cos #:
ð27Þ
It should be mentioned that tSL , tLG and tS are the norms of the forces per unit length in e3 direction and a1 -direction, which act only at the boundary of the meniscus and not over the whole
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
cross-section. In order to determine the total external force at the boundary of the meniscus we compare (27) with the Young equation which appoints the norm of the driving force with: t ¼ r cos #2pr;
ð28Þ
wherein r denotes the surface tension, # is the contact angle between the meniscus and the tube and r represents the radius of the tube. Motivated by this, for the average quantities ^tSL and ^tLG we proclaim that ^tSL ¼ nS rSL d cos #; ^tLG ¼ nG rLG d cos #:
ð29Þ
Herein, rSL and rSL with rSL rLG ¼ r denotes the surface tension between the solid and liquid phase and the liquid and gas phase, respectively. The material parameter d denotes a average value depending on the form and nature of the pores in the porous body. This new approach to the interaction forces allows for the capillary motion in porous solids. 4.2.2. Relative permeability The last terms on the right side of (17) and (18) denote the friction influence on each of the fluids in relation to their seepage velocities. Thereby, SL and SG concern the interaction between the solid phase and the constituents ub and they are restricted to a positive value. These material parameter functions are postulated in dependence on the volume fraction nL with well-known relations by van Genuchten [1], which read in a modified way as: m 2 SL1 ¼ k0L Se1=2 1 1 Se1=m ; ð30Þ 1 G 1=2 1=2 1=m 2m 1 Se : SG ¼ k 0 1 Se 2
In (30)1;2 , k0b ¼ kFb =½ðnL Þ cbR represents a specific soil parameter related to the Darcy permeability coefficient kFb and the specific weight cbR . Moreover, Se ¼ nL =nF denotes the saturation and m is a van Genuchten parameter. Note that SL and SG only describe the relative permeability concerning to the mutational volume fractions. Do not mistake these for the capillary pressure–saturation relation mentioned before.
709
The parameter S b in (30)1;2 characterizes the friction influence of liquid and gas velocities on each other. However, with regard to the investigated fluids these friction rates will be neglected in the following contribution. With these formulations, the set of constitutive equations is complete and we proceed to the numerical modelling.
5. Numerical treatment Within the framework of a standard Galerkin procedure, weak formulations are formed. Therefore, one can either choose the saturation–density formulation or the pressure–pressure formulation, depending on the problem being solved. In the saturation–density formulation we gain the realistic density of the gas qGR , the volume fraction of the liquid nL and the velocities x0b of ub as primary variables and in the pressure–pressure formulation we receive the pressures of both phases pLR and pGR as unknowns. The most important difference between the two formulations can be found in the numerical cost. Because of the reduction from 6 to 2 unknowns the pressure–pressure formulation is very effective. However, this has to be paid for by neglecting dynamic effects. Moreover, the numerical benchmarks like stability and convergence are different concerning the two formulations. Currently, investigations to this points are treated by the authors. In the following sections, we will present the weak formulations of the two versions. The structure of the linearized weak formulations will be given according to the general form: _ Dh_ ¼ R þ RoB ; DG Dh þ DG
ð31Þ
wherein DG multiplied by the vector of the increase of process variables, Dh represents the share _ multiplied by the vector of time of stiffness, DG derivation of increase of the process variables, Dh_ denotes the vector of is the share of damping, R known quantities and RoB represents the vector concerning boundary conditions on the surface oB of the porous body B.
710
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
2
5.1. Saturation–density formulation In the saturation–density formulation we attain the set of unknowns RSD with RSD ¼ RSD ðx; tÞ ¼ fnL ; qGR ; x0L ; x0G g:
ð32Þ
In order to determine these quantities, the weak formulations of the balance of mass and balance of momentum of the liquid and gas phases can be used. They read as follows: Balance of mass for the liquid phase multiplied by the weight function dnL : Z 0 ML ¼ fðnL ÞS þ nL div x0L gdnL dv ¼ 0: ð33Þ B
Balance of mass for the gas phase multiplied by the weight function dnG : Z 0 0 MG ¼ fnG ðqGR ÞG þ ðnG ÞG qGR B
þ nG qGR div x0G gdqGR dv ¼ 0:
ML 6 MG DG ¼ 6 4 IL IG
ð34Þ
Balance of momentum for the liquid phase multiplied by the weight function dx0L : Z IL ¼ fðnL pLR ÞI grad dx0L þ qLR nL ðb x00L Þdx0L B
þ pLR grad nL dx0L qLR ð^tSL þ ^tLG Þgrad nL dx0L Z SL x0L dx0L g dv ¼ fnL pLR n dx0L g da: oBL
ð35Þ Balance of momentum for the gas phase multiplied by the weight function dx0G : Z IG ¼ fðnG pGR ÞI grad dx0G þ qGR nG ðb x00G Þdx0G
2
ML 6 MG _ ¼6 DG 4 0 0
0 MG IL IG
ML 0 IL IG
0 MG 0 0
0 0 IL 0
3 32 0 dnL 6 GR 7 MG 7 76 dq 0 7; 5 IL 4 dxL 5 IG dx0G 32 3 0 dnL 6 GR 7 0 7 76 dq 0 7; 5 0 4 dxL 5 IG dx0G
ð38Þ
2
RoB
3 0 6 0 7 7 ¼6 4 RIL 5: RIG
ð37Þ
ð39Þ
In (39), the boundary conditions RIL ¼ nL pLR n and RIG ¼ nG pGR n denote the partial pressure on the surface of the porous body and result from the surface share in (35) and (36). 5.2. Pressure–pressure formulation The second possibility to reformulate the governing set of equations is the pressure–pressure formulation. Therefore, we modify the set of unknowns to RPP ¼ RPP ðx; tÞ ¼ fpLR ; pGR g:
ð40Þ
In order to reduce the set of equations, the velocities of the constituents x0L and x0G have to be explicitly determined using the balance of momentum in the Darcy formulation: x0L ¼ SL1 fnL grad pLR þ qL b ð^tSL þ ^tLG Þgrad qL g; ð41Þ
B
þ pGR grad nG dx0G qLR^tLG grad nL dx0G Z SG x0G dx0G g dv ¼ fnG pGR n dx0G g da: oBG
ð36Þ In Eqs. (35) and (36), Neumann boundary conditions are formulated on the surface oBL and oBG. Herein, n is the outward oriented unit surface normal on the respective boundary. _ and the vectors Dh The matrix of DG and DG and Dh_ as well as the vector concerning the boundary conditions RoB are given as:
x0G ¼ SG1 fnG grad pGR þ qG b ^tLG grad qL g:
ð42Þ
It must be noted that in this formulation one has to neglect the acceleration x00b . After some modifications the weak formulations of the balances of mass with the weight function dpLR and dpGR reads as: Z 0 fðnL ÞS dpLR grad dpLR nL x0L g dv ML ¼ B Z ¼ fnL x0L dpLR ng da; ð43Þ oBL
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
MG ¼ ¼
Z
6. Numerical example
0
fðqGR nG ÞS dpGR grad dpGR qG x0G g dv
ZB
fqGR nG x0G dpGR ng da:
The linearized weak equations of the saturation–density and pressure–pressure formulation are implemented in the well-known finite element program FEAP by Taylor. The finite element discretization was realized with the Galerkin method. In the example, see Fig. 2(a), we simulate the moistening of a brick standing with one of the small sides in a water reservoir. The left side is closed by an isolation for the liquid and gas phase but the other sides are open. Fig. 2(b) shows the finite element discretization with 403 isoparametric rectangular elements. Due to the expected high saturation gradient at the lower right corner of the brick, the mesh becomes refined at this area. In an experiment, K€ unzel [2] measured the distribution of the volume fraction after 72 days, see Fig. 2(c). The solution of the numerical simulation are depicted in the following. In Fig. 3, the calculated distribution of the liquid versus time (color chart) and the vectors of the capillary force are given. The result of the calculation after 72 days in Fig. 3(c) shows a good approach to the measurement.
ð44Þ
oBG
Thus, for the linearized weak formulation we gain: LR ML ML dp DG ¼ ; ð45Þ MG MG dpGR _ ¼ DG
R
oB
ML MG
ML MG
dpLR ; dpGR
ð46Þ
RML ¼ : RMG
ð47Þ
For the Neumann boundary conditions we have got on the Neumann boundaries oBL a liquid volume efflux RML ¼ nL x0L n and on oBG a gas mass efflux RMG ¼ qGR nG x0G n. In the following example, we use the pressure– pressure formulation for the simulation with linear shape functions concerning the fluid pressures pLR and pGR .
modell
711
measurement
discretization 50
[cm]
49 0.05
impermeable 0.4
40
0.3
30
Height[cm]
open
g 0.2
0.10
0.15
20
0.1
10
0
0
volume fraction nL
0.20
0.25
0
suction
(a)
0
0.1 [cm]
(b) Fig. 2. Saturation of a brick.
5
11 10
Width [cm]
(c)
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
0.5
0.5
0.5 L
L
n [-] 0.255874 0.233806 0.211737 0.199324 0.18829 0.177256 0.164842 0.153808 0.142774 0.13174 0.120706 0.109672 0.0986374 0.0872236 0.0751898 0.0600179 0.0549228 0.053676 0.0521647
y
0.3
0.2
0.1
0.4
0.255874 0.233806 0.211737 0.199324 0.18829 0.177256 0.164842 0.153808 0.142774 0.13174 0.120706 0.109672 0.0986374 0.0872236 0.0751898 0.0600179 0.0549228 0.053676 0.0521647
0.3
y
0.4
L
n [-]
0.2
0.1
0
n [-] 0.4
0.3
0.2
0.1
0
0 0
0.2
x
0.255874 0.233806 0.211737 0.199324 0.18829 0.177256 0.164842 0.153808 0.142774 0.13174 0.120706 0.109672 0.0986374 0.0872236 0.0751898 0.0600179 0.0549228 0.053676 0.0521647
y
712
0
(a)
0
0.2
x
0.2
x
(b)
(c)
Fig. 3. Liquid saturation and capillarity force versus time. Volume fraction nL after 1 min (a), after 8 days (b), and after 72 days (c).
0.5
0.5
L
L
p [Pa] 0.4 2.58346 2.5819 2.56332 2.47346 2.29679 2.12011 1.94343 1.76676 1.59008 1.41341 1.23673 1.06006 0.883379 0.706704 0.530028 0.353352 0.176676 0
y
0.3
0.2
0.1
0.3
0.2
0.1
0
0
0
0.2
x
p [Pa] 0.4
2.58346 2.5819 2.56332 2.47346 2.29679 2.12011 1.94343 1.76676 1.59008 1.41341 1.23673 1.06006 0.883379 0.706704 0.530028 0.353352 0.176676 0
y
0.4
L
p [Pa]
2.58346 2.5819 2.56332 2.47346 2.29679 2.12011 1.94343 1.76676 1.59008 1.41341 1.23673 1.06006 0.883379 0.706704 0.530028 0.353352 0.176676 0
0.3
y
0.5
0.2
0.1
0 0
0.2
x
(a)
0
0.2
x
(b)
(c)
Fig. 4. Liquid pressure and liquid velocity direction versus time. Liquid pressure pL after 1 min (a), after 8 days (b), and after 72 days (c).
The distribution of the liquid pressure (color chart) and the velocity direction of the liquid phase (unit vectors) are presented in Fig. 4. One can recognize that the moving saturation front causes
a pressure swelling in the liquid phase, see Fig. 4(b) and (c). In the case of dynamic effects, like the quick capillary saturation of tissues with saturation
T. Ricken, R. de Boer / Computational Materials Science 28 (2003) 704–713
times of several seconds, instead of the pressure– pressure formulation the saturation–pressure formulation has to be used. Furthermore, the authors are interested in wave propagations in multiphase saturated porous materials. Therefore, only the saturation–pressure formulation can be applied. Readers who are interested in an example with dynamic effects are referred to [7].
7. Conclusion A closed calculation concept for the simulation of multiphase flow in capillary active porous media has been developed. After a brief introduction to the theory of porous media, the used field equations, the constitutive equations and the numerical treatment have been presented. Herein, an essential contribution is given by implementation of the capillary attraction into the interaction forces on the basis of an thermodynamical investigation. To this subject, there have been many empirical approaches but only a limited number of continuum mechanical treatments are known. The presented contribution by the authors in the framework of a ‘‘smeared’’ continua has led to
713
reasonable field equations, in particular, to proper constitutive equations for the interaction volume forces which cause the capillary flow. References [1] M.T. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Science Society of America Journal 44 (1980) 892–898. [2] H.M. K€ unzel, Verfahren zur ein- und zweidimensionalen Berechnung des gekoppelten W€arme- und Feuchtetransports in Bauteilen mit einfachen Kennwerten, Dissertation, Lehrstuhl Konstruktive Bauphysik, Fakult€at Bauingenieurund Vermessungswesen der Universit€at Stuttgart, 1994. [3] H.T. Davis, Statistical Mechanics of Phases, Interfaces and Thin Films, VCH Publishers, New York, Weinheim, Cambridge, 1996. [4] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems, Springer-Verlag, Berlin, Heidelberg, New York, 1997. [5] R. de Boer, Theory of Porous Media: Highlights in the Historical Development and Current State, Springer-Verlag, Berlin, Heidelberg, New York, 2000. [6] R. de Boer, A.K. Didwania, Capillarity in porous bodies: contribution of the Vienna school and recent findings, Acta Mechanica 159 (2002) 173–188. [7] T. Ricken, Kapillarit€at in por€ osen medien––theoretische untersuchung und numerische simulation, Dissertation, Shaker Verlag, Aachen, 2002.