- Email: [email protected]
Finite element approximation for single-phase multicomponent flows
MECHANICS RESEARCH COMMUNICATIONS
Mechanics Research Communications 32 (2005) 53–64 www.elsevier.com/locate/mechrescom
Finite element approximation for single-phase multicomponent flows Fl avia Franceschini 1, Sergio Frey
*
Laboratory of Applied and Computational Fluid Mechanics (LAMAC), Mechanical Engineering Department, Universidade Federal do Rio Grande do Sul, Rua Sarmento Leite 425, 90050-170 Porto Alegre (RS), Brazil Available online 17 May 2004
Abstract With the objective of performing computational simulations of mixture problems, this work employed a mechanical modeling of single-phase incompressible multicomponent flows able to deal with geometric and material non-linearities. This model is based on conservative laws of mass and momentum in continuum mechanics, assuming the mixture as a superposition of a number of continuous media, each one with a variable volume fraction field. Using appropriate hypothesis, the model was formulated as a set of non-linear partial differential equations and associated boundary conditions. This set was approximated by a stabilized finite element method, based on a Galerkin/least-squares scheme, in order to circumvent Babuska–Brezzi condition and to generate stable approximations even in highly advective situations. Some preliminary numerical results have been obtained for non-Newtonian axial injections in backward facing step incompressible flows. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Continuum mixture theory; Multicomponent flows; Finite elements in fluids
1. Introduction This work has a main goal of simulating mixture flow problems––also called multicomponent flows––via stabilized finite element approximations. This kind of phenomenon finds application in many subjects in engineering. In petroleum engineering, modeling and simulation of mixture volumes that arise between batch transfers of petroleum fractions in multiproduct pipelines has been studied by Rachid et al. (2002). In environmental engineering, problems of liquid effluent dispersion in lakes and rivers are common. The prediction of the locals affected by the pollutant and the environmental impact caused are the focus of researches on those multicomponent systems (Curran, 1981). In food industry, problems of in-line mixture of food colorants and additives, clean-in-place devices, mixture tanks as homogenizers and bioreactors are
*
Corresponding author. Fax: +55-51-33163355. E-mail address: [email protected] (S. Frey). 1 This author is a graduate student at Mechanical Engineering Graduate Program (PROMEC/UFRGS). 0093-6413/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2004.04.003
54
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
also examples of multicomponent systems, in which modeling and computational simulation could be of great interest (Geankoplis, 1995). Initially in this paper, a mechanical modeling for multicomponent incompressible flows in a single phase was reviewed (Atkin and Craine, 1976; Sampaio and Williams, 1979; Slattery, 1999). The mechanical model describes mathematically a material body B consisting of N species, which might be undergoing an arbitrary number of homogeneous chemical reactions, which might experiment a condition of stress distribution and answer to this condition as a generalized Newtonian fluid. The resulting boundary value problem was approximated by a stabilized finite element method. The Galerkin method in fluids suffers from two major difficulties. First, the need to satisfy Babuska–Brezzi condition (Ciarlet, 1978) in order to achieve a compatible combination of velocity and pressure subspaces. Further, the inherent instability of central difference schemes in approximating advective dominated equations (Brooks and Hughes, 1982)––in the case considered herein characterized by high Reynolds and mass Peclet numbers. The stabilization technique employed herein was based on a Galerkin/least-squares scheme (Franca and Frey, 1992; Hughes et al., 1986). This method has the ability to circumvent Babuska–Brezzi condition and to generate stable approximations for highly advective flows preserving good accuracy properties. This is achieved by adding residual-based terms to the classic Galerkin formulation, retaining its weighted residual structure and not damaging its consistency.
2. Continuum mechanics conservation and constitutive equations A general mixture theory may be particularized into a mechanical model for incompressible liquid mixtures. This model takes a mixture of N components as the superposition of the N continuous single material bodies BðaÞ , which are mapped into Euclidean space e by the pair placement and volume fraction (vðaÞ ; uðaÞ ), as follows: vðaÞ : BðaÞ Rþ ! e uðaÞ : BðaÞ Rþ ! ½0; 1
ð1Þ
A configuration for a material particle PðaÞ of a body BðaÞ is given by x ¼ vðaÞ ðPðaÞ Þ
ð2Þ
The family of deformations suffered for the material body, parameterized by the time t, from a reference configuration XðaÞ , defines its motion (Slattery, 1999) x ¼ vXð aÞ XðaÞ ; t ð3Þ The velocity and acceleration fields are defined for each component by the first and second time derivatives of the placement. In the Eulerian point of view, these fields are described as functions of a fixed position in space uðaÞ ¼ uðaÞ ðx; tÞ aðaÞ ¼ aðaÞ ðx; tÞ
ð4Þ
The single-component bodies Bð1Þ ; Bð2Þ ; . . . ; BðN Þ are given motions whose trajectories may overlap, then it is required that (Sampaio and Williams, 1979) N X a¼1
uðaÞ ðxÞ ¼
N X dVðaÞ ¼1 dV a¼1
ð5Þ
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
55
~ðaÞ denoting the density of species a in its pure state from the definition of volume fraction it Now, for q follows that ~ðaÞ ðxÞuðaÞ ðxÞ ¼ qðaÞ ðxÞ q
ð6Þ
The mass density of species a with respect to mixture’s mass density in a position x is called mass fraction of a, and denoted by xðaÞ xðaÞ ðxÞ ¼
qðaÞ ðxÞ dmðaÞ ðxÞ=dV ¼ dMðxÞ=dV qðxÞ
ð7Þ
where q represents the density field for the multicomponent body. From the assumption of the mixture as the superposition of the N continuous single material bodies BðaÞ , one assumes the field q as qðxÞ
N X
qðaÞ ðxÞ
ð8Þ
a¼1
and the mass averaged velocity or bulk velocity as: u
N N X 1X qðaÞ uðaÞ ¼ xðaÞ uðaÞ q a¼1 a¼1
ð9Þ
Consequently, from Eq. (7), one achieves that N X
xðaÞ ðxÞ ¼ 1
ð10Þ
a¼1
The relative velocity for the species a with respect to bulk velocity is (uðaÞ u). The N relative velocities are related to each other, as follows from Eqs. (8) and (9), by N X
qðaÞ ðuðaÞ uÞ ¼ 0
ð11Þ
a
Thus, Eq. (11) shows that the definition of bulk velocity is so as to ensure that the sum of the mass flux of all diffusive movements is null. The mass flux of species a with respect to bulk velocity, called relative mass flux, denoted by jðaÞ , is defined as ~ðaÞ uðaÞ ðuðaÞ uÞ ¼ qðaÞ ðuðaÞ uÞ ¼ qxðaÞ ðuðaÞ uÞ jðaÞ q
ð12Þ
and it is finally concluded from Eq. (11) that N X
jðaÞ ¼ 0
ð13Þ
a¼1
As in engineering applications there is a major interest in obtaining results for velocity and pressure fields for the mixture as a whole (called bulk velocity and bulk pressure), a mechanical modeling developed in terms of these bulk variables and also in terms of each component volume or mass fraction field would seem the most proper. Thus, the model exploited in this work is consisted of two equations––in order to solve the bulk velocity and pressure fields––plus N 1 equations for a mixture of N components, that allow the solution of N mass fraction fields using Eq. (10) in sequence. This system of equations is derived from conservation principles approached in sequence, for the multicomponent body B that occupies the region of volume X of boundary C.
56
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
2.1. Principle of mass conservation for a single-component body BðaÞ The time rate of change of species a in a region X of B equals the rate in which this species is produced or consumed by chemical reactions. In the case of a single-phase, this principle is written as Z Z d qðaÞ dX ¼ rðaÞ dX ð14Þ dt X X where rðaÞ denotes the rate of production or consumption of a by homogeneous chemical reactions per unit volume. Applying Reynolds transport theorem (Truesdell and Toupin, 1960) to Eq. (14), results Z dðaÞ qðaÞ ð15Þ þ qðaÞ div uðaÞ rðaÞ dX ¼ 0 dt X With a localization argument (Gurtin, 1981), the differential mass balance for species a is derived as dðaÞ qðaÞ þ qðaÞ div uðaÞ ¼ rðaÞ dt
ð16Þ
In terms of the relative mass flux defined by Eq. (12) and the mass fraction (Eq. (7)) for an incompressible species, Eq. (16) results in q
dðuÞ xðaÞ þ div jðaÞ ¼ rðaÞ dt
ð17Þ
where dðuÞ ðÞ=dt represents the multicomponent material derivative (Slattery, 1999). Applying, in Eq. (16), a summation over all the N species of the mixture, the overall differential mass balance for the multicomponent body B is obtained, in terms of bulk velocity u and mixture density q dðuÞ q þ qdiv u ¼ 0 dt
ð18Þ
For an isochoric motion (Slattery, 1999), where q is assumed constant Eq. (18) becomes div u ¼ 0
ð19Þ
2.2. Principle of momentum conservation for a single-component body BðaÞ The rate of change of momentum of species a in a material volume X relative to an inertial frame of reference equals the resultant force upon X. In this case, body forces are represented by f ðaÞ per unit mass of a, and contact forces by the partial stress vector tðaÞ ðnÞ, defined over C. In addition, the effects of interactions between mixture constituents must be considered (Atkin and Craine, 1976). The momentum supplied to a due to chemical reactions with other constituents is denoted by JðaÞ , and the momentum supplied to species a due to other interactions, such as the relative motion between the constituents because of a diffusive force, is denoted by the force pðaÞ per unit mass of a. The postulate is mathematically written as Z Z Z Z d q uðaÞ dX ¼ rðaÞ JðaÞ dX þ qðaÞ ðf ðaÞ þ pðaÞ Þ dX þ tðaÞ ðnÞ dC ð20Þ dt X ðaÞ X X C Applying Reynolds transport theorem (Truesdell and Toupin, 1960) to Eq. (20), results Z Z dðaÞ uðaÞ qðaÞ tðaÞ ðnÞ dC f ðaÞ þ rðaÞ JðaÞ pðaÞ dX ¼ dt X C
ð21Þ
Following Cauchy’s hypothesis for a single-component body, the existence of a partial stress tensor TðaÞ associated with species a is stated, such as (Atkin and Craine, 1976)
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
tðaÞ ðnÞ ¼ TðaÞ n
57
ð22Þ
Green’s transformation (Slattery, 1999) and localization theorem (Gurtin, 1981) are employed in Eq. (21) to transform the integral formulation in the differential form of the momentum conservation for species a qðaÞ
dðaÞ uðaÞ ¼ div TðaÞ þ pðaÞ rðaÞ ðuðaÞ JðaÞ Þ þ qðaÞ f ðaÞ dt
ð23Þ
To postulate the momentum balance for the mixture as a whole, the following hypotheses have been assumed, according to Atkin and Craine (1976): f¼
N 1X q f ðaÞ q a¼1 ðaÞ
T ¼
N X
ð24Þ
TðaÞ
a¼1 N X
ðpðaÞ þ rðaÞ JðaÞ Þ ¼ 0
a¼1
where f is the total external force per unit mass acting on the mixture, T is the total stress tensor, and the last relation in Eq. (24) accounts for the fact that the sum of all momentum interactions by chemical reactions and relative motions is null. A summation over the N species in Eq. (23) results that X N
qðaÞ
dðaÞ uðaÞ ¼ divT þ qf; dt
ð25Þ
which is approximated by q
dðuÞ u ¼ div T þ qf dt
ð26Þ
where T stands for the total stress plus the apparent stress tensor that arises from diffusion (Atkin and Craine, 1976). Eq. (26) then represents the First Cauchy’s law written for a multicomponent material body. As T is taken to be symmetric, then the moment of momentum conservation is also satisfied (Gurtin, 1981; Atkin and Craine, 1976). 2.3. Constitutive equation for the relative mass flux A classic discussion on the diffusion in multicomponent systems may be found in Slattery (1999) and references therein on the context of the Chapman–Enskog solution for Boltzmann equation to get to an expression for jðaÞ . This expression describes the dependence of jðaÞ on thermal effects, pressure, external forces and concentration gradients. In binary mixtures, Curtiss equation is reduced to the following form: ðoÞ
ðpÞ
ðfÞ
ðTÞ
jðaÞ ¼ jðaÞ þ jðaÞ þ jðaÞ þ jðaÞ ðoÞ
ð27Þ ðpÞ
where jðaÞ , represents the ordinary diffusion due to mass concentration gradients, jðaÞ the mass flux due to ðfÞ ðTÞ pressure effects, jðaÞ the mass flux due to external forces and jðaÞ the mass flux due to thermal effects. According to Slattery (1999), these last three terms may be neglected in many cases of interest; they account for only a small portion in the mass diffusion mechanism. Fundamentally, the term of ordinary diffusion is ðoÞ a function of each species’ concentration gradient. In addition, jðaÞ may be function of parameters such as species physicochemical properties, and also function of mixture properties and characteristics (for a
58
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
detailed discussion see for instance Bird (1960)). In this work, the term of ordinary diffusion is taken as in Fick Law (Bird, 1960), as follows: ðoÞ
jðaÞ ¼ qDðabÞ rxðaÞ
ð28Þ ðoÞ
where the term DðabÞ accounts for the dependence of jðaÞ on parameters as cited above. 2.4. Constitutive equations for T The First Cauchy Law (Eq. (26)) is a conservation law that describes the dynamical conditions of any material body. The specific material behavior is given in the form that the stress tensor T depends on the motion and deformation suffered by the body. This dependence is described by constitutive equations. For a single-component body, one single constitutive equation is necessary to describe its behavior, while for a multicomponent body, it may be necessary to have more some more information in order to describe the mixture’s material behavior in any concentration of its constituents. The model presented in this work, a generalized Newtonian behavior is assumed, which means a linear relation between T and the rate of deformation tensor D. Generalized Newtonian models are largely applied to describe fluid in viscometric flows, where there is the absence of normal stresses or elastic effects (Bird et al., 1987). In those models, a constitutive equation for T is given by T ¼ pI þ 2gðcÞD
ð29Þ
The function gðcÞ is the apparent viscosity, which depends on the model employed to describe the fluid; and c is the magnitude of tensor D. The apparent viscosity for a generalized Newtonian fluid is defined as gap ¼
ð1=2 tr S2 Þ1=2 2 1=2
ð2 tr D Þ
¼
r c
ð30Þ
where the scalar r represents the shear magnitude. In this study, the employed generalized Newtonian model is the constitutive equation known as Carreau model. It is usually used to describe the behavior of blood and many polymer solutions. The relation below gives the Carreau apparent viscosity ! 2 ðn1Þ=2 g c ¼ gc ~ g þ 1 þ W~c ð31Þ where its parameters have been rewritten in a dimensionless form as suggested by Tran-Cahn and TranCong (2002) g c ¼ g0 g1 ;
~ g¼
g1 ; gc
~c ¼
c ; u1 =L
W¼k
u 1
L
;
g ðcÞ ¼
gðcÞ gðcÞ
ð32Þ
The parameter gc is used to characterize the Reynolds number of the flow, with the indices 0 and 1 denoting the limiting viscosity on zero and infinitive rate of deformation. The parameters u1 and L denote a velocity and a length of reference, k is a time parameter and W denotes the Weissenberg number (Bird et al., 1987). In most practical cases, the power law coefficient n in this model is less than one. In those cases, the fluid is called shear–thinning and its apparent viscosity decreases as the shear stress increases.
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
59
3. Finite element approximation A dimensionless boundary-value problem for a generalized Newtonian fluid mixture with and jðaÞ given by Fick law (Eq. (28)) is given by þ ½ruu ¼ rp þ 2Re1 div½g ðcÞDðuÞ þ Fr2 f div u ¼ 0 oxðaÞ þ rxðaÞ u Pe1 m divðrxðaÞ Þ RðaÞ ¼ 0 ot ou ot
in X ð0; t1 Þ in X ð0; t1 Þ in X ð0; t1 Þ
ð33Þ
subjected to the following boundary conditions: u ¼ ug xðaÞ ¼ xðaÞg Tn ¼ th jðaÞ n ¼ jðaÞh
on Cg ð0; t1 Þ on Cg ð0; t1 Þ
ð34Þ
on Ch ð0; t1 Þ on Ch ð0; t1 Þ
ð35Þ
where u is the bulk velocity, p the bulk pressure, g ðcÞ is the non-Newtonian dynamic viscosity normalized with gc , T the stress tensor defined by Eq. (29), DðuÞ the rate of deformation tensor, xðaÞ the mass fraction field and jðaÞ the relative mass flux and RðaÞ the mass generation rate for species a. The Cg and Ch are the portions of boundary C with outward unit normal n on which Dirichlet and Neumann conditions are prescribed. The dimensionless Reynolds, Froude and mass Peclet numbers are given by Re ¼
qu1 L ; gc
u1 Fr ¼ pffiffiffiffiffi ; fL
Pem ¼
u1 L DðabÞ
ð36Þ
where q and gc are mixture’s density and Carreau viscosity defined by Eq. (32), respectively, L and u1 the characteristic length and velocity, f a characteristic body force, and DðabÞ the binary mass diffusivity coefficient. 3.1. A stabilized finite element formulation The variational formulation of Eqs. (33) was approximated by a finite element method using a stabilization scheme which was based on Galerkin/least-squares formulations developed in (Hughes et al., 1986; Franca and Frey, 1992 and references therein). The approximation is performed over the partition Ch of the closed domain X, consisting of convex quadrilateral Q1 elements in Rnsd¼2 In order to approximate the velocity, pressure and mass fraction fields, it was employed the usual spaces of functions (Ciarlet, 1978) for fluid dynamics
n o nsd
nsd Vh ¼ v 2 H01 ðXÞ vjK 2 Rk ðKÞ ; K 2 XK ð37Þ
Ph ¼ p 2 C 0 ðXÞ \ L20 ðXÞ pjK 2 Rl ðKÞ; K 2 XK
ð38Þ
Wh ¼ w 2 H01 ðXÞ wjK 2 Rk ðKÞ; K 2 XK
ð39Þ
n o
Vgh ¼ vð; tÞ 2 H 1 ðXÞnsd ; t 2 ½0; t1 vjK 2 Rk ðKÞnsd ; K 2 XK ; vð; tÞ ¼ ug on Cg
ð40Þ
n o
Whg ¼ wð; tÞ 2 H 1 ðXÞ; t 2 ½0; t1 wjK 2 Rk ðKÞ; K 2 XK ; wð; tÞ ¼ xðaÞg on Cg
ð41Þ
where Rk and Rl denote the polynomial spaces of degrees k and l, respectively.
60
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
Based on the finite element subspaces defined by Eqs. (37)–(41), a stabilized formulation for the model g h h h presented in Eqs. (33)–(36), may be written: Find the triple u ; p ; xðaÞ 2 Vh ; Ph Whg such as B uh ; ph ; xhðaÞ ; v; q; w ¼ F ðv; q; wÞ; ðv; q; wÞ 2 Vh Ph Wh ð42Þ where B u; p; xðaÞ ; v; q; w ¼
ou ; v þ ð½ruu; vÞ þ ð2Re1 g ðcÞDðuÞ; DðvÞÞ ðdiv v; pÞ ðdiv u; qÞ ot oxðaÞ ; w þ ðrxðaÞ u; wÞ þ ðPe1 þ m rxðaÞ ; rwÞ ot X ou þ þ ½ruu þ rp 2Re1 divðg ðcÞDðuÞÞ; sv ðReK Þð½rvu ot K2Ch 1 2Re divðg ðcÞDðuÞÞ rqÞ K X oxðaÞ 1 2 m 1 2 þ þ rxðaÞ u Pem r xðaÞ ; sw ðPeK Þðrw u Pem r wÞ ot K K2Ch ð43Þ
and F ðv; q; wÞ ¼ ðFr2 f; vÞ þ ðth ; vÞC þ ðRðaÞ ; wÞ þ ðjðaÞh ; wÞC þ 1
2Re divðg DðvÞÞ rqÞÞK þ
X
X
ðFr2 f; sv ðReK Þð½rvu
K2Ch 2 ðRðaÞ ; sw ðPemK Þðrw u Pe1 m r wÞÞK
ð44Þ
K2Ch
with the symbol ð; Þ representing the L2 -inner product in X and C, and ð; ÞK the L2 -inner product in the element domain K. The terms inside the sums over the element domains are controlled by the stabilization parameters sv ðReK Þ and sw ðPemK Þ, which according to Franca and Frey (1992) and Franca et al. (1992) were defined as follows, where hK is a mesh length parameter: 8 mik jujp hK > > < for i ¼ v ¼ Re K hK 4m ð45Þ si ðXi Þ ¼ nðXi Þ; with Xi ¼ i m juj h > 2jujp p K k > for i ¼ w : PemK ¼ 2DðabÞ
Xi ; 1;
nðXi Þ ¼ mik ¼ min
( jujp ¼ Ckv
X K2Ch
0 6 Xi < 1 Xi P 1
1 ; 2Cki 3
ð46Þ
ð47Þ
1=p p ju j ; j j¼1 maxj¼1;nsd juj j; Pnsd
2
16p < 1 p¼1 2
h2K kdiv DðvÞk0;K 6 kDðvÞk0
v 2 Vh
ð48Þ
ð49Þ
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
Ckv
X
h2K kr2 vk20;K 6 krvk20
v 2 Wh
61
ð50Þ
K2Ch
Remark 1. Constants Ck derives from inverse estimate (Ciarlet, 1978). As indicated in Franca et al. (1992), their values correspond to Ck ¼ 1 for a bilinear element and Ck ¼ 1=12 for a biquadratic one. New designs of stabilization parameters are found in Tezduyar and Osawa (2000) and Harari et al. (2001). Tezduyar and Osawa (2000) describe the computation of stability parameters based on element-level matrices and vectors, both for advection–diffusion and Navier–Stokes time-dependent equations. Harari et al. (2001) analyze the influence of the flow direction in flow computations, and suggest new definitions of stabilization parameters, incorporating the flow direction and showing good performance. Remark 2. For analysis goals, we must study a linearized version of the formulation introduced by Eqs. (33)–(36). Based on Oseen improvement for the Stokes solution (Landau and Lifchitz, 1971), the balance equations of this system reduce to the following equations subjected to appropriated boundary conditions: ðruÞva þ rp 2Re1 div DðuÞ Fr2 f ¼ 0 div u ¼ 0 rxðaÞ va Pe1 m div rxðaÞ RðaÞ ¼ 0
in X in X in X
ð51Þ
where va is a given advection field, linearizing and uncoupling both momentum and mass transport equations. Setting the stabilization parameters sv and sw to zero, the stabilized formulation (42)–(50) of system (51) reduces to the classical Galerkin approximation for Oseen version of system (33)–(36). The instability of Galerkin approximation of linearized model (51) for advective dominated flows is due to the nsd lack of coercivity of its trilinear form. Choosing ðv; q; wÞ ¼ ðu; p; xðaÞ Þ and selecting v 2 H01 ðXÞ and w 2 H01 ðXÞ, it follows that 2
2
Bðv; q; w; v; q; wÞ P 2kmDðvÞk0 þ DðabÞ krwk0
ðv; q; wÞ 2 Vh Ph Wh
ð52Þ
Therefore, when m and DðabÞ simultaneously tend to zero, the linearized Galerkin approximation will be polluted by spurious oscillations, generating physically unreal solutions. 3.2. Matrix problem Discretization of Eqs. (42)–(50) is carried out by expanding the trial functions u, p, xðaÞ and v, q, w in terms of their finite element basis or shape functions. This leads to the following set of semi-discrete equations: ½M þ Mvs a þ NðuÞ þ Nvs ðuÞ þ K þ Kvs u þ G þ Gvs p ¼ F þ Fs q Ms a þ Nvs ðuÞ þ ½GT þ Kqs u þ Gqs p ¼ Es ð53Þ w w w M þ Ms x_ ðaÞ þ NðuÞ þ Ns ðuÞ xðaÞ þ K þ Ks xðaÞ ¼ F þ Fs where u, p and xðaÞ are the degrees of freedom for uh , ph and xhðaÞ , respectively, a is the vector of degrees of freedom for the transient term ouh =ot, x_ ðaÞ is the vector of degrees of freedom for the transient term oxhðaÞ =ot. The matrices ½M, ½K, ½G arise from the transient, viscous and pressure terms, respectively, and ½M, ½K arise from the transient and diffusive terms. NðuÞ and NðuÞ arise from the advective transport terms. The matrices on the right-hand side are the font terms, and all the other matrices correspond to the stabilized terms. The time integration algorithm is based on Tezduyar et al. (1992), and is a predictor/ corrector scheme for the pressure–velocity problem.
62
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
4. Numerical results In this section, some numerical results for the finite element approximation (Eqs. (42)–(50)) of Eqs. (33)– (36) are presented. The problem analyzed is the axial injection of a fluid (species one) into second fluid (species zero) flowing over a backward facing step––as may be found in Fig. 1. As a motivation for this study, there is the characterization of the interface contamination formed between batch transfers of two different fluids and its spread along the flow. Petroleum transport through long pipelines is an example of application for that (Rachid et al., 2002). The bulk properties were considered constant, the constitutive equations for a Carreau fluid were employed (Eq. (31)), and the mass diffusion was assumed to be governed by Fick Law (Eq. (28)) with a constant mass diffusion coefficient. The imposed boundary conditions are non-slippery walls, constant mass fraction for species one equals to 1.0 at entrance region, free traction and free mass flux for species one at outflow boundary. For initial condition, it is assumed, at instant t0 , a mixture consisting only of species one is injected into the other species. The bulk Reynolds number employed was equal to 100 and the physical properties were assumed independent of species’ concentration. Three orders of magnitude for mass Peclet number have been tested, Pem ¼ 50, 500 and 5000 in order to exploit a large range of interest. The mean velocity at the entrance of the channel u0 was set to 1 m/s and the characteristic length H h was set equal to 0.5 m. For a finite element mesh consisting of 4600 equal-order bilinear elements, for Pem ¼ 5000, the maximum of the mass grid Peclet number occurs in the regions of maximum velocity (u ¼ 1:5u0 ), at the center of channel. There, a mass grid Peclet number has orders of magnitude calculated as PemK ¼
jujp hK 2:12E þ 02 2Dab
ð54Þ
The following pictures show the mass fraction evolution over the domain, with the time evolution based on the non-dimensional time parameter t , t ¼
tu0 ðH hÞ
ð55Þ
Figs. 2 and 3 show the feature of the stabilized method defined by Eqs. (42)–(50) to approximate diffusive and advective-dominated problems, generating results with a physically realistic meaning and not yielding spurious oscillations. The hydrodynamic results obtained were similar to those generated by the Streamline Upwind/Petrov–Galerkin formulation implemented in FLOTRAN code of Ansys Incorporation at Supercomputer Center of Federal University of Rio Grande do Sul (Franceschini, 2002). The interface contamination or mixture front is the region where mass fractions of both species can be found. The mixture front formed in the diffusive case (Pem ¼ 50) is longer and smoother than that in the advective case (Pem ¼ 5000), as it is noticed comparing Figs. 2(a) and 3(a). That is easily explained using Fick Law: for higher mass diffusion coefficients, the relation between mass flux and mass fraction space derivatives happens because of. This result would make one expect that the amount of contaminated fluid
Fig. 1. Problem statement.
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
63
Fig. 2. Result for Pem ¼ 50. Time steps: (a) t ¼ 2; (b) t ¼ 30 and (c) t ¼ 60.
Fig. 3. Result for Pem ¼ 5000. Time steps: (a) t ¼ 2, (b) t ¼ 150 and (c) t ¼ 290.
generated after a process like this would be much smaller for higher Pem . But comparing Figs. 2(b) and (c) and 3(b) and (c), and analyzing the time parameter, one finds out the opposite. Fig. 2(b) and (c) shows that the contamination is largely spread into the flow, but when t ¼ 60, there is only a small fraction of species zero still dispersed. After that species one will easily occupy the whole flow region. In the advective case, the contaminated interface does not spread along the flow, but even after a much longer period (t ¼ 150 and t ¼ 290) there is still a considerable amount of species zero in the flow. That happens because the mixture is retained in zones of recirculation and momentum boundary layer, forming a mass boundary layer near the channel walls. In those regions, where velocities tend to vanish (at momentum boundary layer) or do not take the same direction as the main flow (at the recirculation zone), the mixing of species in the flow direction strongly depends on the diffusion mechanism, which is very inefficient in those high Peclet number cases. Real liquid mixtures usually present very low diffusive characteristics (Geankoplis, 1995) tending to behave more like the flow illustrated in Fig. 3. In cases when there is the rejection of a contaminated volume, for instance in batch transfers in long pipelines of the petroleum industry, the mass boundary layer generated in advective-dominated cases might hinder the procedure for obvious economic reasons.
Acknowledgements The authors F. Franceschini and S. Frey gratefully acknowledge the financial support provided by the agencies CAPES and CNPq (grant 350747/93-8), respectively.
64
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
References Atkin, R.J., Craine, R.E., 1976. Continuum theories of mixtures: basic theory and historical development.Q. J. Mech. Appl. Math. 29 (Pt. 2). Bird, R.B., 1960. Transport Phenomena. John Wiley & Sons, New York. Bird, R.B., Armstrong, R.C., Hassager, O., 1987. In: Dynamics of Polymeric Liquids, vol. 1. John Wiley & Sons, USA. Brooks, A.N., Hughes, T.J.R., 1982. Streamline upwind/Petrov–Galerkin formulations for convective dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comp. Meth. Appl. Mech. Eng. 32, 199–259. Ciarlet, P.G., 1978. The Finite Element Method for Elliptic Problems. North Holland, Amsterdam. Curran, J.C., 1981. A finite-element model of pollution in the Clyde estuary: formulation, validation and utilization. Appl. Math. Model. 5 (3), 137–142. Franca, L.P., Frey, S., 1992. Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Comp. Meth. Appl. Mech. Eng. 99, 209–233. Franca, L.P., Frey, S., Hughes, T.J.R., 1992. Stabilized finite element methods: I. Application to the advective–diffusive model. Comp. Meth. Appl. Mech. Eng. 95, 253–276. Franceschini, F., 2002. Modelagem mec^anica e aproximacß~ao por metodos estabilizados de escoamentos multicomponentes (Mechanical modeling and approximation via stabilized methods for multicomponent flows). Dissertation at PROMEC/UFRGS, Porto Alegre, Brasil, 2002. Geankoplis, C.J., 1995. Procesos de transporte y operaciones unitarias, second ed. Compa~ nıa Editorial Continental S. A., Mexico. Gurtin, M.E., 1981. An Introduction to Continuum Mechanics. Academic Press, New York, USA. Harari, I., Franca, L.P., Oliveira, S., 2001. Streamline design of stability parameter for advective–diffusion problems. J. Comput. Phys. 171, 115–131. Hughes, T.J.R., Franca, L.P., Balestra, M., 1986. A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuska–Brezzi condition: A stable Petrov–Galerkin formulation of the Stokes problem accommodating equalorder interpolations. Comp. Meth. Appl. Mech. Eng. 59, 85–99. Landau, L., Lifchitz, E., 1971. Mecanique des Fluides. Mir, Moscou. Rachid, F.B.F., de Araujo, J.H.C., Baptista, R.M., 2002. Predicting mixing volumes in serial transport in pipelines. J. Fluids Eng. Trans. ASME 124, 528–534. Sampaio, R., Williams, W.O., 1979. Thermodynamics of diffusing mixtures. J. Mech., 18. Slattery, J.C., 1999. Advanced Transport Phenomena. Cambridge University Press, USA. Tezduyar, T.E., Mittal, S., Ray, S.E., Shih, R., 1992. Incompressible flow computations with stabilized bilinear and linear equal-orderinterpolation velocity–pressure elements. Int. J. Numer. Meth. Fluids 95, 221–242. Tezduyar, T.E., Osawa, Y., 2000. Finite element stabilization parameters computed from element matrices and vectors. Comp. Meth. Appl. Mech. Eng. 190, 411–430. Tran-Cahn, D., Tran-Cong, T., 2002. BEM-NN computation of generalized Newtonian flows. Eng. Anal. Bound. Elem. 26, 15–28. Truesdell, C., Toupin, R.A., 1960. The classical field theories. In: Flugge, S. (Ed.), Handbuch Der Physic, vol. 3/1. Springer-Verlag, Berlin.
Mechanics Research Communications 32 (2005) 53–64 www.elsevier.com/locate/mechrescom
Finite element approximation for single-phase multicomponent flows Fl avia Franceschini 1, Sergio Frey
*
Laboratory of Applied and Computational Fluid Mechanics (LAMAC), Mechanical Engineering Department, Universidade Federal do Rio Grande do Sul, Rua Sarmento Leite 425, 90050-170 Porto Alegre (RS), Brazil Available online 17 May 2004
Abstract With the objective of performing computational simulations of mixture problems, this work employed a mechanical modeling of single-phase incompressible multicomponent flows able to deal with geometric and material non-linearities. This model is based on conservative laws of mass and momentum in continuum mechanics, assuming the mixture as a superposition of a number of continuous media, each one with a variable volume fraction field. Using appropriate hypothesis, the model was formulated as a set of non-linear partial differential equations and associated boundary conditions. This set was approximated by a stabilized finite element method, based on a Galerkin/least-squares scheme, in order to circumvent Babuska–Brezzi condition and to generate stable approximations even in highly advective situations. Some preliminary numerical results have been obtained for non-Newtonian axial injections in backward facing step incompressible flows. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Continuum mixture theory; Multicomponent flows; Finite elements in fluids
1. Introduction This work has a main goal of simulating mixture flow problems––also called multicomponent flows––via stabilized finite element approximations. This kind of phenomenon finds application in many subjects in engineering. In petroleum engineering, modeling and simulation of mixture volumes that arise between batch transfers of petroleum fractions in multiproduct pipelines has been studied by Rachid et al. (2002). In environmental engineering, problems of liquid effluent dispersion in lakes and rivers are common. The prediction of the locals affected by the pollutant and the environmental impact caused are the focus of researches on those multicomponent systems (Curran, 1981). In food industry, problems of in-line mixture of food colorants and additives, clean-in-place devices, mixture tanks as homogenizers and bioreactors are
*
Corresponding author. Fax: +55-51-33163355. E-mail address: [email protected] (S. Frey). 1 This author is a graduate student at Mechanical Engineering Graduate Program (PROMEC/UFRGS). 0093-6413/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2004.04.003
54
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
also examples of multicomponent systems, in which modeling and computational simulation could be of great interest (Geankoplis, 1995). Initially in this paper, a mechanical modeling for multicomponent incompressible flows in a single phase was reviewed (Atkin and Craine, 1976; Sampaio and Williams, 1979; Slattery, 1999). The mechanical model describes mathematically a material body B consisting of N species, which might be undergoing an arbitrary number of homogeneous chemical reactions, which might experiment a condition of stress distribution and answer to this condition as a generalized Newtonian fluid. The resulting boundary value problem was approximated by a stabilized finite element method. The Galerkin method in fluids suffers from two major difficulties. First, the need to satisfy Babuska–Brezzi condition (Ciarlet, 1978) in order to achieve a compatible combination of velocity and pressure subspaces. Further, the inherent instability of central difference schemes in approximating advective dominated equations (Brooks and Hughes, 1982)––in the case considered herein characterized by high Reynolds and mass Peclet numbers. The stabilization technique employed herein was based on a Galerkin/least-squares scheme (Franca and Frey, 1992; Hughes et al., 1986). This method has the ability to circumvent Babuska–Brezzi condition and to generate stable approximations for highly advective flows preserving good accuracy properties. This is achieved by adding residual-based terms to the classic Galerkin formulation, retaining its weighted residual structure and not damaging its consistency.
2. Continuum mechanics conservation and constitutive equations A general mixture theory may be particularized into a mechanical model for incompressible liquid mixtures. This model takes a mixture of N components as the superposition of the N continuous single material bodies BðaÞ , which are mapped into Euclidean space e by the pair placement and volume fraction (vðaÞ ; uðaÞ ), as follows: vðaÞ : BðaÞ Rþ ! e uðaÞ : BðaÞ Rþ ! ½0; 1
ð1Þ
A configuration for a material particle PðaÞ of a body BðaÞ is given by x ¼ vðaÞ ðPðaÞ Þ
ð2Þ
The family of deformations suffered for the material body, parameterized by the time t, from a reference configuration XðaÞ , defines its motion (Slattery, 1999) x ¼ vXð aÞ XðaÞ ; t ð3Þ The velocity and acceleration fields are defined for each component by the first and second time derivatives of the placement. In the Eulerian point of view, these fields are described as functions of a fixed position in space uðaÞ ¼ uðaÞ ðx; tÞ aðaÞ ¼ aðaÞ ðx; tÞ
ð4Þ
The single-component bodies Bð1Þ ; Bð2Þ ; . . . ; BðN Þ are given motions whose trajectories may overlap, then it is required that (Sampaio and Williams, 1979) N X a¼1
uðaÞ ðxÞ ¼
N X dVðaÞ ¼1 dV a¼1
ð5Þ
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
55
~ðaÞ denoting the density of species a in its pure state from the definition of volume fraction it Now, for q follows that ~ðaÞ ðxÞuðaÞ ðxÞ ¼ qðaÞ ðxÞ q
ð6Þ
The mass density of species a with respect to mixture’s mass density in a position x is called mass fraction of a, and denoted by xðaÞ xðaÞ ðxÞ ¼
qðaÞ ðxÞ dmðaÞ ðxÞ=dV ¼ dMðxÞ=dV qðxÞ
ð7Þ
where q represents the density field for the multicomponent body. From the assumption of the mixture as the superposition of the N continuous single material bodies BðaÞ , one assumes the field q as qðxÞ
N X
qðaÞ ðxÞ
ð8Þ
a¼1
and the mass averaged velocity or bulk velocity as: u
N N X 1X qðaÞ uðaÞ ¼ xðaÞ uðaÞ q a¼1 a¼1
ð9Þ
Consequently, from Eq. (7), one achieves that N X
xðaÞ ðxÞ ¼ 1
ð10Þ
a¼1
The relative velocity for the species a with respect to bulk velocity is (uðaÞ u). The N relative velocities are related to each other, as follows from Eqs. (8) and (9), by N X
qðaÞ ðuðaÞ uÞ ¼ 0
ð11Þ
a
Thus, Eq. (11) shows that the definition of bulk velocity is so as to ensure that the sum of the mass flux of all diffusive movements is null. The mass flux of species a with respect to bulk velocity, called relative mass flux, denoted by jðaÞ , is defined as ~ðaÞ uðaÞ ðuðaÞ uÞ ¼ qðaÞ ðuðaÞ uÞ ¼ qxðaÞ ðuðaÞ uÞ jðaÞ q
ð12Þ
and it is finally concluded from Eq. (11) that N X
jðaÞ ¼ 0
ð13Þ
a¼1
As in engineering applications there is a major interest in obtaining results for velocity and pressure fields for the mixture as a whole (called bulk velocity and bulk pressure), a mechanical modeling developed in terms of these bulk variables and also in terms of each component volume or mass fraction field would seem the most proper. Thus, the model exploited in this work is consisted of two equations––in order to solve the bulk velocity and pressure fields––plus N 1 equations for a mixture of N components, that allow the solution of N mass fraction fields using Eq. (10) in sequence. This system of equations is derived from conservation principles approached in sequence, for the multicomponent body B that occupies the region of volume X of boundary C.
56
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
2.1. Principle of mass conservation for a single-component body BðaÞ The time rate of change of species a in a region X of B equals the rate in which this species is produced or consumed by chemical reactions. In the case of a single-phase, this principle is written as Z Z d qðaÞ dX ¼ rðaÞ dX ð14Þ dt X X where rðaÞ denotes the rate of production or consumption of a by homogeneous chemical reactions per unit volume. Applying Reynolds transport theorem (Truesdell and Toupin, 1960) to Eq. (14), results Z dðaÞ qðaÞ ð15Þ þ qðaÞ div uðaÞ rðaÞ dX ¼ 0 dt X With a localization argument (Gurtin, 1981), the differential mass balance for species a is derived as dðaÞ qðaÞ þ qðaÞ div uðaÞ ¼ rðaÞ dt
ð16Þ
In terms of the relative mass flux defined by Eq. (12) and the mass fraction (Eq. (7)) for an incompressible species, Eq. (16) results in q
dðuÞ xðaÞ þ div jðaÞ ¼ rðaÞ dt
ð17Þ
where dðuÞ ðÞ=dt represents the multicomponent material derivative (Slattery, 1999). Applying, in Eq. (16), a summation over all the N species of the mixture, the overall differential mass balance for the multicomponent body B is obtained, in terms of bulk velocity u and mixture density q dðuÞ q þ qdiv u ¼ 0 dt
ð18Þ
For an isochoric motion (Slattery, 1999), where q is assumed constant Eq. (18) becomes div u ¼ 0
ð19Þ
2.2. Principle of momentum conservation for a single-component body BðaÞ The rate of change of momentum of species a in a material volume X relative to an inertial frame of reference equals the resultant force upon X. In this case, body forces are represented by f ðaÞ per unit mass of a, and contact forces by the partial stress vector tðaÞ ðnÞ, defined over C. In addition, the effects of interactions between mixture constituents must be considered (Atkin and Craine, 1976). The momentum supplied to a due to chemical reactions with other constituents is denoted by JðaÞ , and the momentum supplied to species a due to other interactions, such as the relative motion between the constituents because of a diffusive force, is denoted by the force pðaÞ per unit mass of a. The postulate is mathematically written as Z Z Z Z d q uðaÞ dX ¼ rðaÞ JðaÞ dX þ qðaÞ ðf ðaÞ þ pðaÞ Þ dX þ tðaÞ ðnÞ dC ð20Þ dt X ðaÞ X X C Applying Reynolds transport theorem (Truesdell and Toupin, 1960) to Eq. (20), results Z Z dðaÞ uðaÞ qðaÞ tðaÞ ðnÞ dC f ðaÞ þ rðaÞ JðaÞ pðaÞ dX ¼ dt X C
ð21Þ
Following Cauchy’s hypothesis for a single-component body, the existence of a partial stress tensor TðaÞ associated with species a is stated, such as (Atkin and Craine, 1976)
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
tðaÞ ðnÞ ¼ TðaÞ n
57
ð22Þ
Green’s transformation (Slattery, 1999) and localization theorem (Gurtin, 1981) are employed in Eq. (21) to transform the integral formulation in the differential form of the momentum conservation for species a qðaÞ
dðaÞ uðaÞ ¼ div TðaÞ þ pðaÞ rðaÞ ðuðaÞ JðaÞ Þ þ qðaÞ f ðaÞ dt
ð23Þ
To postulate the momentum balance for the mixture as a whole, the following hypotheses have been assumed, according to Atkin and Craine (1976): f¼
N 1X q f ðaÞ q a¼1 ðaÞ
T ¼
N X
ð24Þ
TðaÞ
a¼1 N X
ðpðaÞ þ rðaÞ JðaÞ Þ ¼ 0
a¼1
where f is the total external force per unit mass acting on the mixture, T is the total stress tensor, and the last relation in Eq. (24) accounts for the fact that the sum of all momentum interactions by chemical reactions and relative motions is null. A summation over the N species in Eq. (23) results that X N
qðaÞ
dðaÞ uðaÞ ¼ divT þ qf; dt
ð25Þ
which is approximated by q
dðuÞ u ¼ div T þ qf dt
ð26Þ
where T stands for the total stress plus the apparent stress tensor that arises from diffusion (Atkin and Craine, 1976). Eq. (26) then represents the First Cauchy’s law written for a multicomponent material body. As T is taken to be symmetric, then the moment of momentum conservation is also satisfied (Gurtin, 1981; Atkin and Craine, 1976). 2.3. Constitutive equation for the relative mass flux A classic discussion on the diffusion in multicomponent systems may be found in Slattery (1999) and references therein on the context of the Chapman–Enskog solution for Boltzmann equation to get to an expression for jðaÞ . This expression describes the dependence of jðaÞ on thermal effects, pressure, external forces and concentration gradients. In binary mixtures, Curtiss equation is reduced to the following form: ðoÞ
ðpÞ
ðfÞ
ðTÞ
jðaÞ ¼ jðaÞ þ jðaÞ þ jðaÞ þ jðaÞ ðoÞ
ð27Þ ðpÞ
where jðaÞ , represents the ordinary diffusion due to mass concentration gradients, jðaÞ the mass flux due to ðfÞ ðTÞ pressure effects, jðaÞ the mass flux due to external forces and jðaÞ the mass flux due to thermal effects. According to Slattery (1999), these last three terms may be neglected in many cases of interest; they account for only a small portion in the mass diffusion mechanism. Fundamentally, the term of ordinary diffusion is ðoÞ a function of each species’ concentration gradient. In addition, jðaÞ may be function of parameters such as species physicochemical properties, and also function of mixture properties and characteristics (for a
58
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
detailed discussion see for instance Bird (1960)). In this work, the term of ordinary diffusion is taken as in Fick Law (Bird, 1960), as follows: ðoÞ
jðaÞ ¼ qDðabÞ rxðaÞ
ð28Þ ðoÞ
where the term DðabÞ accounts for the dependence of jðaÞ on parameters as cited above. 2.4. Constitutive equations for T The First Cauchy Law (Eq. (26)) is a conservation law that describes the dynamical conditions of any material body. The specific material behavior is given in the form that the stress tensor T depends on the motion and deformation suffered by the body. This dependence is described by constitutive equations. For a single-component body, one single constitutive equation is necessary to describe its behavior, while for a multicomponent body, it may be necessary to have more some more information in order to describe the mixture’s material behavior in any concentration of its constituents. The model presented in this work, a generalized Newtonian behavior is assumed, which means a linear relation between T and the rate of deformation tensor D. Generalized Newtonian models are largely applied to describe fluid in viscometric flows, where there is the absence of normal stresses or elastic effects (Bird et al., 1987). In those models, a constitutive equation for T is given by T ¼ pI þ 2gðcÞD
ð29Þ
The function gðcÞ is the apparent viscosity, which depends on the model employed to describe the fluid; and c is the magnitude of tensor D. The apparent viscosity for a generalized Newtonian fluid is defined as gap ¼
ð1=2 tr S2 Þ1=2 2 1=2
ð2 tr D Þ
¼
r c
ð30Þ
where the scalar r represents the shear magnitude. In this study, the employed generalized Newtonian model is the constitutive equation known as Carreau model. It is usually used to describe the behavior of blood and many polymer solutions. The relation below gives the Carreau apparent viscosity ! 2 ðn1Þ=2 g c ¼ gc ~ g þ 1 þ W~c ð31Þ where its parameters have been rewritten in a dimensionless form as suggested by Tran-Cahn and TranCong (2002) g c ¼ g0 g1 ;
~ g¼
g1 ; gc
~c ¼
c ; u1 =L
W¼k
u 1
L
;
g ðcÞ ¼
gðcÞ gðcÞ
ð32Þ
The parameter gc is used to characterize the Reynolds number of the flow, with the indices 0 and 1 denoting the limiting viscosity on zero and infinitive rate of deformation. The parameters u1 and L denote a velocity and a length of reference, k is a time parameter and W denotes the Weissenberg number (Bird et al., 1987). In most practical cases, the power law coefficient n in this model is less than one. In those cases, the fluid is called shear–thinning and its apparent viscosity decreases as the shear stress increases.
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
59
3. Finite element approximation A dimensionless boundary-value problem for a generalized Newtonian fluid mixture with and jðaÞ given by Fick law (Eq. (28)) is given by þ ½ruu ¼ rp þ 2Re1 div½g ðcÞDðuÞ þ Fr2 f div u ¼ 0 oxðaÞ þ rxðaÞ u Pe1 m divðrxðaÞ Þ RðaÞ ¼ 0 ot ou ot
in X ð0; t1 Þ in X ð0; t1 Þ in X ð0; t1 Þ
ð33Þ
subjected to the following boundary conditions: u ¼ ug xðaÞ ¼ xðaÞg Tn ¼ th jðaÞ n ¼ jðaÞh
on Cg ð0; t1 Þ on Cg ð0; t1 Þ
ð34Þ
on Ch ð0; t1 Þ on Ch ð0; t1 Þ
ð35Þ
where u is the bulk velocity, p the bulk pressure, g ðcÞ is the non-Newtonian dynamic viscosity normalized with gc , T the stress tensor defined by Eq. (29), DðuÞ the rate of deformation tensor, xðaÞ the mass fraction field and jðaÞ the relative mass flux and RðaÞ the mass generation rate for species a. The Cg and Ch are the portions of boundary C with outward unit normal n on which Dirichlet and Neumann conditions are prescribed. The dimensionless Reynolds, Froude and mass Peclet numbers are given by Re ¼
qu1 L ; gc
u1 Fr ¼ pffiffiffiffiffi ; fL
Pem ¼
u1 L DðabÞ
ð36Þ
where q and gc are mixture’s density and Carreau viscosity defined by Eq. (32), respectively, L and u1 the characteristic length and velocity, f a characteristic body force, and DðabÞ the binary mass diffusivity coefficient. 3.1. A stabilized finite element formulation The variational formulation of Eqs. (33) was approximated by a finite element method using a stabilization scheme which was based on Galerkin/least-squares formulations developed in (Hughes et al., 1986; Franca and Frey, 1992 and references therein). The approximation is performed over the partition Ch of the closed domain X, consisting of convex quadrilateral Q1 elements in Rnsd¼2 In order to approximate the velocity, pressure and mass fraction fields, it was employed the usual spaces of functions (Ciarlet, 1978) for fluid dynamics
n o nsd
nsd Vh ¼ v 2 H01 ðXÞ vjK 2 Rk ðKÞ ; K 2 XK ð37Þ
Ph ¼ p 2 C 0 ðXÞ \ L20 ðXÞ pjK 2 Rl ðKÞ; K 2 XK
ð38Þ
Wh ¼ w 2 H01 ðXÞ wjK 2 Rk ðKÞ; K 2 XK
ð39Þ
n o
Vgh ¼ vð; tÞ 2 H 1 ðXÞnsd ; t 2 ½0; t1 vjK 2 Rk ðKÞnsd ; K 2 XK ; vð; tÞ ¼ ug on Cg
ð40Þ
n o
Whg ¼ wð; tÞ 2 H 1 ðXÞ; t 2 ½0; t1 wjK 2 Rk ðKÞ; K 2 XK ; wð; tÞ ¼ xðaÞg on Cg
ð41Þ
where Rk and Rl denote the polynomial spaces of degrees k and l, respectively.
60
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
Based on the finite element subspaces defined by Eqs. (37)–(41), a stabilized formulation for the model g h h h presented in Eqs. (33)–(36), may be written: Find the triple u ; p ; xðaÞ 2 Vh ; Ph Whg such as B uh ; ph ; xhðaÞ ; v; q; w ¼ F ðv; q; wÞ; ðv; q; wÞ 2 Vh Ph Wh ð42Þ where B u; p; xðaÞ ; v; q; w ¼
ou ; v þ ð½ruu; vÞ þ ð2Re1 g ðcÞDðuÞ; DðvÞÞ ðdiv v; pÞ ðdiv u; qÞ ot oxðaÞ ; w þ ðrxðaÞ u; wÞ þ ðPe1 þ m rxðaÞ ; rwÞ ot X ou þ þ ½ruu þ rp 2Re1 divðg ðcÞDðuÞÞ; sv ðReK Þð½rvu ot K2Ch 1 2Re divðg ðcÞDðuÞÞ rqÞ K X oxðaÞ 1 2 m 1 2 þ þ rxðaÞ u Pem r xðaÞ ; sw ðPeK Þðrw u Pem r wÞ ot K K2Ch ð43Þ
and F ðv; q; wÞ ¼ ðFr2 f; vÞ þ ðth ; vÞC þ ðRðaÞ ; wÞ þ ðjðaÞh ; wÞC þ 1
2Re divðg DðvÞÞ rqÞÞK þ
X
X
ðFr2 f; sv ðReK Þð½rvu
K2Ch 2 ðRðaÞ ; sw ðPemK Þðrw u Pe1 m r wÞÞK
ð44Þ
K2Ch
with the symbol ð; Þ representing the L2 -inner product in X and C, and ð; ÞK the L2 -inner product in the element domain K. The terms inside the sums over the element domains are controlled by the stabilization parameters sv ðReK Þ and sw ðPemK Þ, which according to Franca and Frey (1992) and Franca et al. (1992) were defined as follows, where hK is a mesh length parameter: 8 mik jujp hK > > < for i ¼ v ¼ Re K hK 4m ð45Þ si ðXi Þ ¼ nðXi Þ; with Xi ¼ i m juj h > 2jujp p K k > for i ¼ w : PemK ¼ 2DðabÞ
Xi ; 1;
nðXi Þ ¼ mik ¼ min
( jujp ¼ Ckv
X K2Ch
0 6 Xi < 1 Xi P 1
1 ; 2Cki 3
ð46Þ
ð47Þ
1=p p ju j ; j j¼1 maxj¼1;nsd juj j; Pnsd
2
16p < 1 p¼1 2
h2K kdiv DðvÞk0;K 6 kDðvÞk0
v 2 Vh
ð48Þ
ð49Þ
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
Ckv
X
h2K kr2 vk20;K 6 krvk20
v 2 Wh
61
ð50Þ
K2Ch
Remark 1. Constants Ck derives from inverse estimate (Ciarlet, 1978). As indicated in Franca et al. (1992), their values correspond to Ck ¼ 1 for a bilinear element and Ck ¼ 1=12 for a biquadratic one. New designs of stabilization parameters are found in Tezduyar and Osawa (2000) and Harari et al. (2001). Tezduyar and Osawa (2000) describe the computation of stability parameters based on element-level matrices and vectors, both for advection–diffusion and Navier–Stokes time-dependent equations. Harari et al. (2001) analyze the influence of the flow direction in flow computations, and suggest new definitions of stabilization parameters, incorporating the flow direction and showing good performance. Remark 2. For analysis goals, we must study a linearized version of the formulation introduced by Eqs. (33)–(36). Based on Oseen improvement for the Stokes solution (Landau and Lifchitz, 1971), the balance equations of this system reduce to the following equations subjected to appropriated boundary conditions: ðruÞva þ rp 2Re1 div DðuÞ Fr2 f ¼ 0 div u ¼ 0 rxðaÞ va Pe1 m div rxðaÞ RðaÞ ¼ 0
in X in X in X
ð51Þ
where va is a given advection field, linearizing and uncoupling both momentum and mass transport equations. Setting the stabilization parameters sv and sw to zero, the stabilized formulation (42)–(50) of system (51) reduces to the classical Galerkin approximation for Oseen version of system (33)–(36). The instability of Galerkin approximation of linearized model (51) for advective dominated flows is due to the nsd lack of coercivity of its trilinear form. Choosing ðv; q; wÞ ¼ ðu; p; xðaÞ Þ and selecting v 2 H01 ðXÞ and w 2 H01 ðXÞ, it follows that 2
2
Bðv; q; w; v; q; wÞ P 2kmDðvÞk0 þ DðabÞ krwk0
ðv; q; wÞ 2 Vh Ph Wh
ð52Þ
Therefore, when m and DðabÞ simultaneously tend to zero, the linearized Galerkin approximation will be polluted by spurious oscillations, generating physically unreal solutions. 3.2. Matrix problem Discretization of Eqs. (42)–(50) is carried out by expanding the trial functions u, p, xðaÞ and v, q, w in terms of their finite element basis or shape functions. This leads to the following set of semi-discrete equations: ½M þ Mvs a þ NðuÞ þ Nvs ðuÞ þ K þ Kvs u þ G þ Gvs p ¼ F þ Fs q Ms a þ Nvs ðuÞ þ ½GT þ Kqs u þ Gqs p ¼ Es ð53Þ w w w M þ Ms x_ ðaÞ þ NðuÞ þ Ns ðuÞ xðaÞ þ K þ Ks xðaÞ ¼ F þ Fs where u, p and xðaÞ are the degrees of freedom for uh , ph and xhðaÞ , respectively, a is the vector of degrees of freedom for the transient term ouh =ot, x_ ðaÞ is the vector of degrees of freedom for the transient term oxhðaÞ =ot. The matrices ½M, ½K, ½G arise from the transient, viscous and pressure terms, respectively, and ½M, ½K arise from the transient and diffusive terms. NðuÞ and NðuÞ arise from the advective transport terms. The matrices on the right-hand side are the font terms, and all the other matrices correspond to the stabilized terms. The time integration algorithm is based on Tezduyar et al. (1992), and is a predictor/ corrector scheme for the pressure–velocity problem.
62
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
4. Numerical results In this section, some numerical results for the finite element approximation (Eqs. (42)–(50)) of Eqs. (33)– (36) are presented. The problem analyzed is the axial injection of a fluid (species one) into second fluid (species zero) flowing over a backward facing step––as may be found in Fig. 1. As a motivation for this study, there is the characterization of the interface contamination formed between batch transfers of two different fluids and its spread along the flow. Petroleum transport through long pipelines is an example of application for that (Rachid et al., 2002). The bulk properties were considered constant, the constitutive equations for a Carreau fluid were employed (Eq. (31)), and the mass diffusion was assumed to be governed by Fick Law (Eq. (28)) with a constant mass diffusion coefficient. The imposed boundary conditions are non-slippery walls, constant mass fraction for species one equals to 1.0 at entrance region, free traction and free mass flux for species one at outflow boundary. For initial condition, it is assumed, at instant t0 , a mixture consisting only of species one is injected into the other species. The bulk Reynolds number employed was equal to 100 and the physical properties were assumed independent of species’ concentration. Three orders of magnitude for mass Peclet number have been tested, Pem ¼ 50, 500 and 5000 in order to exploit a large range of interest. The mean velocity at the entrance of the channel u0 was set to 1 m/s and the characteristic length H h was set equal to 0.5 m. For a finite element mesh consisting of 4600 equal-order bilinear elements, for Pem ¼ 5000, the maximum of the mass grid Peclet number occurs in the regions of maximum velocity (u ¼ 1:5u0 ), at the center of channel. There, a mass grid Peclet number has orders of magnitude calculated as PemK ¼
jujp hK 2:12E þ 02 2Dab
ð54Þ
The following pictures show the mass fraction evolution over the domain, with the time evolution based on the non-dimensional time parameter t , t ¼
tu0 ðH hÞ
ð55Þ
Figs. 2 and 3 show the feature of the stabilized method defined by Eqs. (42)–(50) to approximate diffusive and advective-dominated problems, generating results with a physically realistic meaning and not yielding spurious oscillations. The hydrodynamic results obtained were similar to those generated by the Streamline Upwind/Petrov–Galerkin formulation implemented in FLOTRAN code of Ansys Incorporation at Supercomputer Center of Federal University of Rio Grande do Sul (Franceschini, 2002). The interface contamination or mixture front is the region where mass fractions of both species can be found. The mixture front formed in the diffusive case (Pem ¼ 50) is longer and smoother than that in the advective case (Pem ¼ 5000), as it is noticed comparing Figs. 2(a) and 3(a). That is easily explained using Fick Law: for higher mass diffusion coefficients, the relation between mass flux and mass fraction space derivatives happens because of. This result would make one expect that the amount of contaminated fluid
Fig. 1. Problem statement.
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
63
Fig. 2. Result for Pem ¼ 50. Time steps: (a) t ¼ 2; (b) t ¼ 30 and (c) t ¼ 60.
Fig. 3. Result for Pem ¼ 5000. Time steps: (a) t ¼ 2, (b) t ¼ 150 and (c) t ¼ 290.
generated after a process like this would be much smaller for higher Pem . But comparing Figs. 2(b) and (c) and 3(b) and (c), and analyzing the time parameter, one finds out the opposite. Fig. 2(b) and (c) shows that the contamination is largely spread into the flow, but when t ¼ 60, there is only a small fraction of species zero still dispersed. After that species one will easily occupy the whole flow region. In the advective case, the contaminated interface does not spread along the flow, but even after a much longer period (t ¼ 150 and t ¼ 290) there is still a considerable amount of species zero in the flow. That happens because the mixture is retained in zones of recirculation and momentum boundary layer, forming a mass boundary layer near the channel walls. In those regions, where velocities tend to vanish (at momentum boundary layer) or do not take the same direction as the main flow (at the recirculation zone), the mixing of species in the flow direction strongly depends on the diffusion mechanism, which is very inefficient in those high Peclet number cases. Real liquid mixtures usually present very low diffusive characteristics (Geankoplis, 1995) tending to behave more like the flow illustrated in Fig. 3. In cases when there is the rejection of a contaminated volume, for instance in batch transfers in long pipelines of the petroleum industry, the mass boundary layer generated in advective-dominated cases might hinder the procedure for obvious economic reasons.
Acknowledgements The authors F. Franceschini and S. Frey gratefully acknowledge the financial support provided by the agencies CAPES and CNPq (grant 350747/93-8), respectively.
64
F. Franceschini, S. Frey / Mechanics Research Communications 32 (2005) 53–64
References Atkin, R.J., Craine, R.E., 1976. Continuum theories of mixtures: basic theory and historical development.Q. J. Mech. Appl. Math. 29 (Pt. 2). Bird, R.B., 1960. Transport Phenomena. John Wiley & Sons, New York. Bird, R.B., Armstrong, R.C., Hassager, O., 1987. In: Dynamics of Polymeric Liquids, vol. 1. John Wiley & Sons, USA. Brooks, A.N., Hughes, T.J.R., 1982. Streamline upwind/Petrov–Galerkin formulations for convective dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comp. Meth. Appl. Mech. Eng. 32, 199–259. Ciarlet, P.G., 1978. The Finite Element Method for Elliptic Problems. North Holland, Amsterdam. Curran, J.C., 1981. A finite-element model of pollution in the Clyde estuary: formulation, validation and utilization. Appl. Math. Model. 5 (3), 137–142. Franca, L.P., Frey, S., 1992. Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Comp. Meth. Appl. Mech. Eng. 99, 209–233. Franca, L.P., Frey, S., Hughes, T.J.R., 1992. Stabilized finite element methods: I. Application to the advective–diffusive model. Comp. Meth. Appl. Mech. Eng. 95, 253–276. Franceschini, F., 2002. Modelagem mec^anica e aproximacß~ao por metodos estabilizados de escoamentos multicomponentes (Mechanical modeling and approximation via stabilized methods for multicomponent flows). Dissertation at PROMEC/UFRGS, Porto Alegre, Brasil, 2002. Geankoplis, C.J., 1995. Procesos de transporte y operaciones unitarias, second ed. Compa~ nıa Editorial Continental S. A., Mexico. Gurtin, M.E., 1981. An Introduction to Continuum Mechanics. Academic Press, New York, USA. Harari, I., Franca, L.P., Oliveira, S., 2001. Streamline design of stability parameter for advective–diffusion problems. J. Comput. Phys. 171, 115–131. Hughes, T.J.R., Franca, L.P., Balestra, M., 1986. A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuska–Brezzi condition: A stable Petrov–Galerkin formulation of the Stokes problem accommodating equalorder interpolations. Comp. Meth. Appl. Mech. Eng. 59, 85–99. Landau, L., Lifchitz, E., 1971. Mecanique des Fluides. Mir, Moscou. Rachid, F.B.F., de Araujo, J.H.C., Baptista, R.M., 2002. Predicting mixing volumes in serial transport in pipelines. J. Fluids Eng. Trans. ASME 124, 528–534. Sampaio, R., Williams, W.O., 1979. Thermodynamics of diffusing mixtures. J. Mech., 18. Slattery, J.C., 1999. Advanced Transport Phenomena. Cambridge University Press, USA. Tezduyar, T.E., Mittal, S., Ray, S.E., Shih, R., 1992. Incompressible flow computations with stabilized bilinear and linear equal-orderinterpolation velocity–pressure elements. Int. J. Numer. Meth. Fluids 95, 221–242. Tezduyar, T.E., Osawa, Y., 2000. Finite element stabilization parameters computed from element matrices and vectors. Comp. Meth. Appl. Mech. Eng. 190, 411–430. Tran-Cahn, D., Tran-Cong, T., 2002. BEM-NN computation of generalized Newtonian flows. Eng. Anal. Bound. Elem. 26, 15–28. Truesdell, C., Toupin, R.A., 1960. The classical field theories. In: Flugge, S. (Ed.), Handbuch Der Physic, vol. 3/1. Springer-Verlag, Berlin.