- Email: [email protected]
Coupling of free surface and groundwater flows
Computers & Fluids 32 (2003) 73–83 www.elsevier.com/locate/compfluid
Coupling of free surface and groundwater flows Edie Miglio a
a,*
, Alfio Quarteroni
a,b
, Fausto Saleri
c
Dipartimento di Matematica ‘‘F. Brioschi’’, Politecnico di Milano, P.zza Leonardo da Vinci 32, 20133 Milano, Italy b D epart ement de Math ematiques, Ecole Polytechnique F ed erale de Lausanne, CH-1015, Lausanne, Switzerland c Dipartimento di Matematica ‘‘F. Enriques’’, Universita degli Studi di Milano, Via Saldini 50, 20133 Milano, Italy
Abstract In this paper we present some preliminary results about the coupling of shallow water equations for free surface flows and Darcy equation for groundwater flows. A suitable set of interface conditions is discussed: the Beavers and Joseph formula for the bottom stress is used. An iterative algorithm to solve the coupled problem is proposed and some numerical results are presented. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Free surface flows; Shallow water equations; Groundwater flows; Darcy’s law; Porous media; Interface conditions
1. Introduction In this paper we address the problem of the coupling between the surface and the subsurface b p of the motion of a fluid. In particular, we will consider the problem in which the lower portion X b computational domain is formed by a porous medium while in the upper part X f we consider a free surface fluid. The problem in the two regions is governed by different partial differential equations, typically the free surface Navier Stokes equations and the Darcy’s law, giving rise to a coupled problem of heterogeneous type. For the closure of the problem, we need suitable interface conditions at C ¼ b p that relate the different mathematical unknowns from either subdomain across the b f \ oX oX interface. More precisely we demand for the continuity of the normal component of the velocity b p. and the equality of the fluid pressure and the piezometric head in X
*
Corresponding author. E-mail address: [email protected] (E. Miglio).
0045-7930/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 1 0 2 - 5
74
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
This coupling between surface and groundwater flows, furtherly coupled with numerical models of transport-diffusion reaction of chemical agents, has interesting applications, e.g. it could help in assessing the short- and medium-term effects of a pollutant dispersed in a fluid. The plan of the paper is as follows. In Section 2 we describe the mathematical models for the b p . Section 3 is devoted to the b f and of the fluid in the ground X motion of the free surface fluid in X discussion of the interface conditions between the two models; Section 4 deals with the development of an iterative scheme to decouple the solution of the coupled problem and the weak formulation of the problem is derived. Finally in the last section some preliminary numerical results are presented.
2. Problem setting b is divided into two parts: the one occupied by the porous medium denoted by The domain X b f ; the interface between the fluid and b X p and that occupied by the free surface fluid denoted by X the porous media is denoted by C. The reference level correspond with the lower point of the porous domain and hence the free fluid is bounded above by the free surface (Cs ) and below by the bottom described by the function hðx; yÞ (see Fig. 1). In the following, to avoid confusions, b p , respectively. b f and X the subscripts ‘f’ and ‘p’ will denote quantities relative to X As for the porous media we introduce the following quantity, called piezometric head, pp ; U¼zþ qf g where z is the elevation from the reference level i.e., the potential energy per unit weight of water, b p the motion of the fluid is described by the pp is the pressure and qf the density of the fluid. In X Darcy law and the continuity equation [2,17]
b. Fig. 1. Schematic representation of the domain X
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
8 <
oU e q¼0 þr ot : e q ¼ KrU S0
b p; 8t > 0 and 8ðx; y; zÞ 2 X b p; 8t > 0 and 8ðx; y; zÞ 2 X
75
ð1Þ
e is the 3D where K is the permeability tensor, S0 is the specific mass storativity coefficient and r gradient operator. Once the specific storage vector q is known it is possible to calculate the average velocity Vp ðx; y; z; tÞ ¼ ðup ; wp ÞT in the porous media using the relation Vp ¼ q=n, where n is the porosity coefficient (up and wp are the horizontal and vertical component of the velocity vector). b p =C ¼ CpD [ CpN : on CpD we assign the piezometric head while on CpN we e p ¼ oX Finally C assign the normal component of the velocity. The domain in which we consider the free surface fluid is characterized as follows: let X be a bounded domain of R2 which represents the undisturbed free surface of the fluid. Moreover, let z ¼ hðx; yÞ and z ¼ gðx; y; tÞ be the functions describing the bathymetry and the free surface with b f ¼ fðx; y; zÞjðx; yÞ 2 X; respect to the reference level z ¼ 0 (see Fig. 1). The domain of interest is X b f is a normal domain with respect to the z-axis, thus each integral z 2 ð h; gÞg. Hence the domain X R R Rg bf ¼ b ð Þ dz dX: This consideration will be used in Section over X f can be written as X^ f ð Þ d X X h 4.1 for deriving the weak formulation of the problem. b f the motion of the free surface fluid is described by the 3D non-hydrostatic shallow water In X equations (briefly 3D-NH-SW) [4], in which the density is constant and the total pressure is written as the sum of an hydrostatic part and an hydrodynamic correction, i.e., p ¼ qgðg zÞ þ q. The model reads: 8 Duf o ouf > b f; > mv
þ g rg þ rq ¼ f 8t > 0 and 8ðx; y; zÞ 2 X > > > Dt oz oz > > > > Dwf o owf oq > b f; < mv ¼0 8t > 0 and 8ðx; y; zÞ 2 X
þ oz oz Dt oz ð2Þ > owf > b > ¼0 8t > 0 and 8ðx; y; zÞ 2 X f ; r uf þ > > ozg > Z > > og > > e : þr uf dz ¼ Q 8t > 0 and 8ðx; yÞ 2 X; ot h where g is the gravity acceleration, f is the external force vector (for instance the Coriolis force), uf ðx; y; z; tÞ and wf ðx; y; z; tÞ represent the horizontal and vertical components of the velocity T vector (Vf ¼ ðuf ; wf Þ will denote the 3D velocity vector), qðx; y; z; tÞ is the hydrodynamic pressure, gðx; y; tÞ is the elevation from the reference level and mv is the vertical viscosity coefficient. Finally r is the 2D horizontal gradient operator and D=Dt ¼ uf r þ wf o=oz is the total derivative. Let us remind that the last equation of system (2) is used to calculate the free surface position. It is obtained integrating along the vertical direction the continuity equation and using the kinematic boundary condition on the free surface og þ uf rg ¼ wf ot
8t > 0 and 8ðx; yÞ 2 X; z ¼ gðx; y; tÞ;
ð3Þ
as well as the following condition on the bottom surface z ¼ hðx; yÞ e Vf n f ¼ Q
8t > 0 and ðx; yÞ 2 X; z ¼ hðx; yÞ:
ð4Þ
76
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
e is taken equal to zero, i.e., a no-flux condition Usually the bottom surface is impermeable, thus Q through the bottom is imposed, yielding a null right-hand side in the fourth equation of (2). e to be different from zero accounting for the presence However, in view of the coupling we allow Q of sinks and/or sources. In this paper we neglect the effect of the wind on the free surface i.e., ouf ¼ 0 8t > 0 and 8ðx; yÞ 2 X; z ¼ gðx; y; tÞ oz and moreover on a fictitious water–water boundary (Copen that we suppose to be vertical) we assign the normal component of the velocity. In the following nf ¼ ðnxyf ; nzf ÞT and np ¼ ðnxyp ; nzp ÞT will denote the outward unit normal vectors b p , respectively; nxy and nz are the horizontal and vertical components of these normal b f and o X to o X vectors. On the interface C we have nf ¼ np . 3. Interface conditions In 1967 Beavers and Joseph [3], observed that the interfacial velocities of the exterior fluid and of the fluid in the porous medium could be related by an ad hoc matching condition which admits a discontinuity in the tangential component. Using this idea to model the friction effect on the bottom we propose the following interface conditions on C 8 Vf nf ¼ Vp np ; > > > pffiffiffi < ouf aBJ 3 ð5Þ ¼ pffiffiffiffiffiffiffiffi ðuf up Þ; > oz > tr K > : pf ¼ qf gU ¼ qf gH þ pp ; where pf is the pressure exerted by the fluid on the bottom and H ¼ g h is the total height of the b f. fluid in X The parameter aBJ is an empiric dimensionless coefficient whose value depends on the properties of the porous medium. Its range of values is 0.01–5 [10]. The first condition in (5) enforces the continuity of the normal velocities on C. The set of boundary conditions we have advocated is alternative to the following one proposed in [12] for the coupling of the Stokes equations with the Darcy’s law (still to be held on C) 8 Vf nf ¼ Vp np ; > > > pffiffiffi > < ouf aBJ 3 ¼ pffiffiffiffiffiffiffiffi ðuf up Þ; ð6Þ oz tr K > > > > : p ow ¼ p : f p oz It must be noticed that, contrarily to what done by Ene and Sanchez-Palencia [5], Jones [9] and Ochoa-Tapia and Whitaker [11], in both sets of interface conditions the pressure is allowed to be discontinuous across the interface as advocated also by Bayada and Chambat [1] and J€ager and Mikelic [7,8].
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
77
4. Subdomain iterations between surface and subsurface flows To decouple the computation of the physical variables in the free surface fluid domain from that in the porous medium we adopt an iterative scheme which is inspired to the Dirichlet– Neumann method in heterogeneous domain decomposition methods [13]. To this aim the first and the second condition of (5) are assigned to the free surface fluid: in e appearing on the particular the first condition of (5) is enforced through the sink/source term Q e right-hand side of the fourth equation of (2) i.e., Q ¼ Vp np , while the second condition of (5) substitutes the Chezy or Manning formulas accounting for the bottom friction. Finally the third condition of (5) is used as a Dirichlet condition for the porous media. It is worthwhile to notice that these choices correspond to the specification of the tangential component of the normal stress and the normal component of the velocity field for the free surface fluid and a Dirichlet condition on the piezometric head for the porous media. More precisely, the iterative method reads as follows: 8t > 0, given k0 , find for j P 0 the sojþ1 ðtÞ, qjþ1 ðtÞ, Ujþ1 ðtÞg of lution fðVjþ1 f ðtÞ, g 8 oUjþ1 > > e qjþ1 S ¼ r > 0 > < ot e jþ1 qjþ1 ¼ KrU > jþ1 > qgU ¼ kj > > : þ b:c:
b p; in X b p; in X on C; b p =C on o X
ð7Þ
and 8 ! jþ1 jþ1 > Du o ou > f > > mv f ¼ g rgjþ1 rqjþ1 þ > > oz Dt oz > > ! > > jþ1 jþ1 > > Dwf oqjþ1 o ow > > mv f ¼
þ > > oz Dt oz oz > > > > jþ1 < ow þ f ¼0 r ujþ1 f > > Zozgjþ1 > h i jþ1 > og > jþ1 jþ1 > > þ r u dz ¼ V n p > p f > C ot > h pffiffiffi > > jþ1 > ou mv aBJ 3 > > >
ujþ1 mv f ¼ pffiffiffiffiffiffiffiffi ðujþ1 > f p Þ > oz tr K > > : þ b:c:
b f; in X b f; in X b f; in X
ð8Þ
in X; on C; b f =C; on o X
finally setting kjþ1 ¼ hpfjþ1 þ ð1 hÞkj on C where h 2 ½1=2; 1 is a relaxation parameter. As for the time discretization the semi implicit Euler scheme proposed in [4] can be adopted: in this context when solving the problem at time tnþ1 we have g ¼ gðtn Þ and similarly k0 ¼ pf jC ðtn Þ. To summarize, the iterative procedure to solve the coupled problem (7) and (8) reads as follows: for every t, let k0 be the an initial guess for the interface condition on C; then iterate on j till kkjþ1 kj kL2 ðCÞ =kkjþ1 kL2 ðCÞ 6 :
78
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
(i) solve problem (7), (ii) solve problem (8) using Uj computed at step (i), (iii) compute the new value kjþ1 and go back to step (i).
4.1. Weak formulation and finite element approximation In this section we derive the weak formulation of the two problems. The porous domain is decomposed using a mesh (Tp ) of prisms and each prism of this mesh will be denoted by P. We adopt the lowest order Raviart–Thomas finite elements (RT0 ) [14] for the space discretization while the time advancing is achieved using the implicit Euler scheme [15]. Neumann condition. For simplicity we assume CpD ¼ ; and on CpN n n o we take a homogeneous 2 b b p Þjs 2 RT0 ðP Þ; Let Qph and Wph ¼ s h 2 HC~ p ;div ð X hjP o ¼ uh 2 L ð X p ÞjuhjP 2 P0 ðP Þ; 8P 2 Tp 8P 2 Tp
where
b pÞ ¼ HC~ p ;div ð X
3 o b b p Þ and s n ¼ 0 on C ep e s 2 L2 ð X s 2 L ð X p Þ r
2
and P0 is the space of constant functions. In the following j is the index of the iterations by subdomains. The weak form of the problem for the porous media reads: for all t > 0, find ðqh ; Uh Þ 2 Wph Qph such that 8Z Z Z > jþ1 e j
1 jþ1 b b > > < ^ K qh s h d X p ^ Uh r sh d X p þ k sh np dC ¼ 0 8sh 2 Wph ; C Xp X Z p Z ð9Þ jþ1 oUh > jþ1 b > b e > d X þ u d X ¼ 0 8u 2 Q : r q u S p p ph h h h : 0 ^ h ot ^p Xp X The velocity in the porous media can be calculated a-posteriori when the discharge vector is jþ1 known using the relation Vjþ1 ph ¼ qh =n. As for the 3D-NH-SW model the domain is discretized using a multilayer approach [6] which induces a grid Tf of prism; each prism P can be seen as the cartesian product of a triangle T 2 Txy (where Txy is a subdivision into triangles of X) with a vertical interval dz. The time advancing is achieved using a fractional step scheme with a Lagrangian treatment of the convective terms [4]. To derive the weak form of the 3D-NH-SW model we assume Copen ¼ ; and we introduce the following spaces n o b f ; @z Þ; r ^ b f Þj^ bf ; u 2 H 1ðX u 2 L2 ð X u nxy ¼ 0 on o X Vf ¼ ^ ) ^ o w 2 b f Þ; b f Þjw bf ; ^ 2 L ðX ^ ¼ 0 on o X w 2 L ðX oz
( Wf ¼
2
b f ; @z Þ is the following anisotropic Sobolev space [16] where H 1 ð X 2 ow 2 1 b 2 b 2 b ^ ðx; y; zÞ 2 ðL ð X f ÞÞ H ð X f ; @z Þ ¼ w 2 ðL ð X f ÞÞ : oz
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
79
The discrete n spaces are o uh 2 Vf j^ uhjP 2 RT0 ðT Þ P1 ðdzÞ 8P 2 Tf ; Vfh ¼ ^ n o ^ h 2 Wf jw ^ hjP 2 P0 ðT Þ P1 ðdzÞ 8P 2 Tf ; Wfh ¼ w n o ghjT 2 P0 ðT Þ 8T 2 Txy ; Pfh ¼ g^h 2 L2 ðXÞj^ n o b f Þj^ qhjP 2 P0 ðP Þ 8P 2 Tf ; Qfh ¼ q^h 2 L2 ð X where P1 is the space of linear functions in the z variable. The weak form reads: for all t > 0, find ðufh ; wfh ; gh ; qh Þ 2 Vfh Wfh Pfh Qfh such that 8Z Z > Dujþ1 oujþ1 o^ uh b fh > b > ^ uh d X f þ dXf mv fh > > Dt oz oz ^f ^f > X X > p ffiffi ffi > Z > > mv aBJ 3 jþ1 > >
pffiffiffiffiffiffiffiffi ðufh ujþ1 uh dX > ph Þnz ^ > > tr K X > Z Z > > > jþ1 bf
bf > uh d X qjþ1 uh d X > h r^ > ^ ggh r ^ > ^f Xf X > Z > > > > b f 8^ > ¼ f^ uh d X uh 2 Vf ; > > ^ > Xf > Z Dt oz oz ^f ^f > X X Z > > > ^h b ow > > ^ h 2 Wf ;
qjþ1 d X f ¼ 0 8w > h > oz > ^f X > > Z Z > > owjþ1 > jþ1 fh b b f ¼ 0 8^ > r ufh q^h d X f þ qh 2 Qf ; q^ d X > > > oz h ^f ^f X X > > Z Z Z > gh > ogjþ1 > h > ^ ^h dX > g dX þ r ujþ1 h fh dz g > > ot > X X h Z > > > > : ¼ ðVjþ1 gh dX 8^ gh 2 Pf ; ph np Þ^ X
jþ1 þ ð1 hÞkj . and finally k ¼ hpnh The weak formulation of the two problems shows that the interface conditions we have proposed are easily incorporated in the models we have used. jþ1
5. Numerical results In this section we present some preliminary numerical results concerning the coupling of free surface and groundwater flows. In this tests we have considered the situation presented in Fig. 2 in which the interface and the bottom of the porous media are flat. In the first simulation we consider a water channel 1 m wide and 10 m long with a water depth of 1 m covering a porous media whose height is 9 m. At the inlet of the channel a discharge of 10 4 m3 /s is prescribed while at outlet a non-reflective boundary condition is enforced. The porous
80
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
b in the case of a flat interface. Fig. 2. Schematic representation of the domain X
media is homogeneous (K ¼ dI, where d ¼ 10 5 and I is the identity tensor) and aBJ ¼ 1. The time step is set to Dt ¼ 600 s. The mesh for the porous media is composed by 15 320 prisms while the mesh for the free surface flows is composed by 13 788 prisms. In Fig. 3 the velocity profile in the center of the channel is shown: the set of boundary conditions (5) using the Beavers and Joseph model for the bottom friction gives rise to a parabolic velocity profile as expected. In the second test we consider the situation shown in Fig. 4 in which a water basin of 20 m side and 10 m depth covers a porous media whose height is 9 m. In the basin the free surface is initially flat and the velocity is null while in the porous media the initial potential head is 0. The mesh for the porous media is composed by 23 706 prisms while the mesh for the free surface flows is
Fig. 3. Vertical profile of the velocity field in the center of the basin (longitudinal section).
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
81
Fig. 4. Setting for the second test case.
Fig. 5. Vertical section of the velocity profile at the beginning of the simulation in the free surface fluid.
composed by 26 340 prisms. Due to the weight of the water a filtration in the porous media takes place: in Figs. 5 and 6 the vertical section of the velocity profile is shown at the beginning of the simulation and after two days, while Fig. 7 shows the piezometric head and the velocity field after t ¼ 48 h (the scales of the velocity vectors in the figure for the free surface fluid and for the porous media are different). In Fig. 8 the solid and the dashed lines show the water volume in the basin as
Fig. 6. Vertical section of the velocity profile after two days in the free surface fluid.
82
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
Fig. 7. Piezometric head and velocity field after two days.
Fig. 8. Volume balance.
a function of time and the net flux of water through the bottom, respectively: it can be appreciated that the total mass of water is conserved. In the numerical results presented in this paper the iteration by subdomain is performed at each time step. Nevertheless the algorithm we have presented can be used also in those situations where the time scales in the porous media and in the free surface fluid are very different (for instance in presence of porous media with very low permeability). In this case the iteration by subdomain is performed only every N time steps of the free surface fluid and N is such that Dtp ¼ N Dtf , where Dtp and Dtf are the time steps for the porous media and for the free surface fluid, respectively. Remark 1. The numerical results show that the iterative scheme we have proposed is convergent: in particular, we have taken a tolerance ¼ 10 3 and the convergence is obtained in about five iterations. From the theoretical point of view the convergence analysis of this procedure will rely on the formulation of a Steklov–Poincare problem for the unknown U at the interface C and will be the subject of future work.
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
83
6. Conclusion In this paper we have addressed the problem of coupling free surface and groundwater flows. A set of boundary conditions based on physical considerations is proposed; to model the bottom stress in the free surface fluid the Beavers and Joseph formula is employed. A suitable iterative scheme to decouple the solution of the full problem is proposed. The numerical tests show that the iterative scheme is convergent and the effectiveness of the interface conditions we have used.
Acknowledgement We thank Eng. T. Ponta for having performed some simulations.
References [1] Bayada G, Chambat M. On the interface conditions for a thin film flow past a porous medium. SIAM J Math Anal 1995;26:1113–29. [2] Bear J. Hydraulics of groundwater. New York: McGraw-Hill; 1979. [3] Beavers GS, Joseph DD. Boundary conditions at a naturally permeable wall. J Fluid Mech 1967;30:197–207. [4] Causin P, Miglio E, Saleri F. Algebraic factorizations for 3D non-hydrostatic free surface flows. Accepted for publication in Comput Visual Sci. quations et phenomenes de surface pour l’ecoulement dans un modele de milieu [5] Ene HI, Sanchez-Palencia E. E poreaux. J Mecanique 1975;14:73–108. [6] Fontana L, Miglio E, Quarteroni A, Saleri F. A finite element method for 3D hydrostatic water flows. Comput Visual Sci 1999;2(2/3):85–93. [7] J€ager W, Mikelic A. On the boundary conditions at the contact interface between a porous medium and a free fluid. Ann Scuola Norm Sup Pisa Cl Sci 1996;23:403–65. [8] J€ager W, Mikelic A. On the interface boundary condition of Beavers, Joseph and Saffman. SIAM J Appl Math 2000;60:1111–27. [9] Jones IP. Low Reynolds number flow past a porous spherical shell. Proc Cambridge Philos Soc 1973;73:231–8. [10] Nield DA, Bejan A. Convection in porous media. New York: Springer; 1999. [11] Ochoa-Tapia JA, Whitaker S. Momentum transfer at the boundary between a porous medium and a homogeneous fluid. I. Theoretical development. Int J Heat Mass Transfer 1995;38:2635–46. [12] Payne LE, Straughan B. Analysis of the boundary condition at the interface between a viscous fluid and a porous medium and related modelling questions. J Math Pures Appl 1998;77:317–54. [13] Quarteroni A, Valli A. Domain decomposition method for partial differential equations. Oxford: Oxford Science Publications; 1999. [14] Raviart PA, Thomas JM. A mixed finite element method for 2nd order elliptic problems. In: Dold A, Eckmann B, editors. Mathematical aspects of finite elements methods. Lecture Notes 606. Berlin, Heidelberg, New York: Springer; 1977. [15] Thomee V. Galerkin finite elements methods for parabolic problems. Berlin, Heidelberg: Springer; 1997. [16] Temam R. Sur la stabilite et la convergence de la methode des pas fractionnaires. Ann Mat Pura et Appl IV Ser 1968;79:191–379. [17] Wood WL. Introduction to numerical methods for water resources. Oxford: Oxford Science Publications; 1993.
Coupling of free surface and groundwater flows Edie Miglio a
a,*
, Alfio Quarteroni
a,b
, Fausto Saleri
c
Dipartimento di Matematica ‘‘F. Brioschi’’, Politecnico di Milano, P.zza Leonardo da Vinci 32, 20133 Milano, Italy b D epart ement de Math ematiques, Ecole Polytechnique F ed erale de Lausanne, CH-1015, Lausanne, Switzerland c Dipartimento di Matematica ‘‘F. Enriques’’, Universita degli Studi di Milano, Via Saldini 50, 20133 Milano, Italy
Abstract In this paper we present some preliminary results about the coupling of shallow water equations for free surface flows and Darcy equation for groundwater flows. A suitable set of interface conditions is discussed: the Beavers and Joseph formula for the bottom stress is used. An iterative algorithm to solve the coupled problem is proposed and some numerical results are presented. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Free surface flows; Shallow water equations; Groundwater flows; Darcy’s law; Porous media; Interface conditions
1. Introduction In this paper we address the problem of the coupling between the surface and the subsurface b p of the motion of a fluid. In particular, we will consider the problem in which the lower portion X b computational domain is formed by a porous medium while in the upper part X f we consider a free surface fluid. The problem in the two regions is governed by different partial differential equations, typically the free surface Navier Stokes equations and the Darcy’s law, giving rise to a coupled problem of heterogeneous type. For the closure of the problem, we need suitable interface conditions at C ¼ b p that relate the different mathematical unknowns from either subdomain across the b f \ oX oX interface. More precisely we demand for the continuity of the normal component of the velocity b p. and the equality of the fluid pressure and the piezometric head in X
*
Corresponding author. E-mail address: [email protected] (E. Miglio).
0045-7930/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 1 0 2 - 5
74
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
This coupling between surface and groundwater flows, furtherly coupled with numerical models of transport-diffusion reaction of chemical agents, has interesting applications, e.g. it could help in assessing the short- and medium-term effects of a pollutant dispersed in a fluid. The plan of the paper is as follows. In Section 2 we describe the mathematical models for the b p . Section 3 is devoted to the b f and of the fluid in the ground X motion of the free surface fluid in X discussion of the interface conditions between the two models; Section 4 deals with the development of an iterative scheme to decouple the solution of the coupled problem and the weak formulation of the problem is derived. Finally in the last section some preliminary numerical results are presented.
2. Problem setting b is divided into two parts: the one occupied by the porous medium denoted by The domain X b f ; the interface between the fluid and b X p and that occupied by the free surface fluid denoted by X the porous media is denoted by C. The reference level correspond with the lower point of the porous domain and hence the free fluid is bounded above by the free surface (Cs ) and below by the bottom described by the function hðx; yÞ (see Fig. 1). In the following, to avoid confusions, b p , respectively. b f and X the subscripts ‘f’ and ‘p’ will denote quantities relative to X As for the porous media we introduce the following quantity, called piezometric head, pp ; U¼zþ qf g where z is the elevation from the reference level i.e., the potential energy per unit weight of water, b p the motion of the fluid is described by the pp is the pressure and qf the density of the fluid. In X Darcy law and the continuity equation [2,17]
b. Fig. 1. Schematic representation of the domain X
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
8 <
oU e q¼0 þr ot : e q ¼ KrU S0
b p; 8t > 0 and 8ðx; y; zÞ 2 X b p; 8t > 0 and 8ðx; y; zÞ 2 X
75
ð1Þ
e is the 3D where K is the permeability tensor, S0 is the specific mass storativity coefficient and r gradient operator. Once the specific storage vector q is known it is possible to calculate the average velocity Vp ðx; y; z; tÞ ¼ ðup ; wp ÞT in the porous media using the relation Vp ¼ q=n, where n is the porosity coefficient (up and wp are the horizontal and vertical component of the velocity vector). b p =C ¼ CpD [ CpN : on CpD we assign the piezometric head while on CpN we e p ¼ oX Finally C assign the normal component of the velocity. The domain in which we consider the free surface fluid is characterized as follows: let X be a bounded domain of R2 which represents the undisturbed free surface of the fluid. Moreover, let z ¼ hðx; yÞ and z ¼ gðx; y; tÞ be the functions describing the bathymetry and the free surface with b f ¼ fðx; y; zÞjðx; yÞ 2 X; respect to the reference level z ¼ 0 (see Fig. 1). The domain of interest is X b f is a normal domain with respect to the z-axis, thus each integral z 2 ð h; gÞg. Hence the domain X R R Rg bf ¼ b ð Þ dz dX: This consideration will be used in Section over X f can be written as X^ f ð Þ d X X h 4.1 for deriving the weak formulation of the problem. b f the motion of the free surface fluid is described by the 3D non-hydrostatic shallow water In X equations (briefly 3D-NH-SW) [4], in which the density is constant and the total pressure is written as the sum of an hydrostatic part and an hydrodynamic correction, i.e., p ¼ qgðg zÞ þ q. The model reads: 8 Duf o ouf > b f; > mv
þ g rg þ rq ¼ f 8t > 0 and 8ðx; y; zÞ 2 X > > > Dt oz oz > > > > Dwf o owf oq > b f; < mv ¼0 8t > 0 and 8ðx; y; zÞ 2 X
þ oz oz Dt oz ð2Þ > owf > b > ¼0 8t > 0 and 8ðx; y; zÞ 2 X f ; r uf þ > > ozg > Z > > og > > e : þr uf dz ¼ Q 8t > 0 and 8ðx; yÞ 2 X; ot h where g is the gravity acceleration, f is the external force vector (for instance the Coriolis force), uf ðx; y; z; tÞ and wf ðx; y; z; tÞ represent the horizontal and vertical components of the velocity T vector (Vf ¼ ðuf ; wf Þ will denote the 3D velocity vector), qðx; y; z; tÞ is the hydrodynamic pressure, gðx; y; tÞ is the elevation from the reference level and mv is the vertical viscosity coefficient. Finally r is the 2D horizontal gradient operator and D=Dt ¼ uf r þ wf o=oz is the total derivative. Let us remind that the last equation of system (2) is used to calculate the free surface position. It is obtained integrating along the vertical direction the continuity equation and using the kinematic boundary condition on the free surface og þ uf rg ¼ wf ot
8t > 0 and 8ðx; yÞ 2 X; z ¼ gðx; y; tÞ;
ð3Þ
as well as the following condition on the bottom surface z ¼ hðx; yÞ e Vf n f ¼ Q
8t > 0 and ðx; yÞ 2 X; z ¼ hðx; yÞ:
ð4Þ
76
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
e is taken equal to zero, i.e., a no-flux condition Usually the bottom surface is impermeable, thus Q through the bottom is imposed, yielding a null right-hand side in the fourth equation of (2). e to be different from zero accounting for the presence However, in view of the coupling we allow Q of sinks and/or sources. In this paper we neglect the effect of the wind on the free surface i.e., ouf ¼ 0 8t > 0 and 8ðx; yÞ 2 X; z ¼ gðx; y; tÞ oz and moreover on a fictitious water–water boundary (Copen that we suppose to be vertical) we assign the normal component of the velocity. In the following nf ¼ ðnxyf ; nzf ÞT and np ¼ ðnxyp ; nzp ÞT will denote the outward unit normal vectors b p , respectively; nxy and nz are the horizontal and vertical components of these normal b f and o X to o X vectors. On the interface C we have nf ¼ np . 3. Interface conditions In 1967 Beavers and Joseph [3], observed that the interfacial velocities of the exterior fluid and of the fluid in the porous medium could be related by an ad hoc matching condition which admits a discontinuity in the tangential component. Using this idea to model the friction effect on the bottom we propose the following interface conditions on C 8 Vf nf ¼ Vp np ; > > > pffiffiffi < ouf aBJ 3 ð5Þ ¼ pffiffiffiffiffiffiffiffi ðuf up Þ; > oz > tr K > : pf ¼ qf gU ¼ qf gH þ pp ; where pf is the pressure exerted by the fluid on the bottom and H ¼ g h is the total height of the b f. fluid in X The parameter aBJ is an empiric dimensionless coefficient whose value depends on the properties of the porous medium. Its range of values is 0.01–5 [10]. The first condition in (5) enforces the continuity of the normal velocities on C. The set of boundary conditions we have advocated is alternative to the following one proposed in [12] for the coupling of the Stokes equations with the Darcy’s law (still to be held on C) 8 Vf nf ¼ Vp np ; > > > pffiffiffi > < ouf aBJ 3 ¼ pffiffiffiffiffiffiffiffi ðuf up Þ; ð6Þ oz tr K > > > > : p ow ¼ p : f p oz It must be noticed that, contrarily to what done by Ene and Sanchez-Palencia [5], Jones [9] and Ochoa-Tapia and Whitaker [11], in both sets of interface conditions the pressure is allowed to be discontinuous across the interface as advocated also by Bayada and Chambat [1] and J€ager and Mikelic [7,8].
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
77
4. Subdomain iterations between surface and subsurface flows To decouple the computation of the physical variables in the free surface fluid domain from that in the porous medium we adopt an iterative scheme which is inspired to the Dirichlet– Neumann method in heterogeneous domain decomposition methods [13]. To this aim the first and the second condition of (5) are assigned to the free surface fluid: in e appearing on the particular the first condition of (5) is enforced through the sink/source term Q e right-hand side of the fourth equation of (2) i.e., Q ¼ Vp np , while the second condition of (5) substitutes the Chezy or Manning formulas accounting for the bottom friction. Finally the third condition of (5) is used as a Dirichlet condition for the porous media. It is worthwhile to notice that these choices correspond to the specification of the tangential component of the normal stress and the normal component of the velocity field for the free surface fluid and a Dirichlet condition on the piezometric head for the porous media. More precisely, the iterative method reads as follows: 8t > 0, given k0 , find for j P 0 the sojþ1 ðtÞ, qjþ1 ðtÞ, Ujþ1 ðtÞg of lution fðVjþ1 f ðtÞ, g 8 oUjþ1 > > e qjþ1 S ¼ r > 0 > < ot e jþ1 qjþ1 ¼ KrU > jþ1 > qgU ¼ kj > > : þ b:c:
b p; in X b p; in X on C; b p =C on o X
ð7Þ
and 8 ! jþ1 jþ1 > Du o ou > f > > mv f ¼ g rgjþ1 rqjþ1 þ > > oz Dt oz > > ! > > jþ1 jþ1 > > Dwf oqjþ1 o ow > > mv f ¼
þ > > oz Dt oz oz > > > > jþ1 < ow þ f ¼0 r ujþ1 f > > Zozgjþ1 > h i jþ1 > og > jþ1 jþ1 > > þ r u dz ¼ V n p > p f > C ot > h pffiffiffi > > jþ1 > ou mv aBJ 3 > > >
ujþ1 mv f ¼ pffiffiffiffiffiffiffiffi ðujþ1 > f p Þ > oz tr K > > : þ b:c:
b f; in X b f; in X b f; in X
ð8Þ
in X; on C; b f =C; on o X
finally setting kjþ1 ¼ hpfjþ1 þ ð1 hÞkj on C where h 2 ½1=2; 1 is a relaxation parameter. As for the time discretization the semi implicit Euler scheme proposed in [4] can be adopted: in this context when solving the problem at time tnþ1 we have g ¼ gðtn Þ and similarly k0 ¼ pf jC ðtn Þ. To summarize, the iterative procedure to solve the coupled problem (7) and (8) reads as follows: for every t, let k0 be the an initial guess for the interface condition on C; then iterate on j till kkjþ1 kj kL2 ðCÞ =kkjþ1 kL2 ðCÞ 6 :
78
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
(i) solve problem (7), (ii) solve problem (8) using Uj computed at step (i), (iii) compute the new value kjþ1 and go back to step (i).
4.1. Weak formulation and finite element approximation In this section we derive the weak formulation of the two problems. The porous domain is decomposed using a mesh (Tp ) of prisms and each prism of this mesh will be denoted by P. We adopt the lowest order Raviart–Thomas finite elements (RT0 ) [14] for the space discretization while the time advancing is achieved using the implicit Euler scheme [15]. Neumann condition. For simplicity we assume CpD ¼ ; and on CpN n n o we take a homogeneous 2 b b p Þjs 2 RT0 ðP Þ; Let Qph and Wph ¼ s h 2 HC~ p ;div ð X hjP o ¼ uh 2 L ð X p ÞjuhjP 2 P0 ðP Þ; 8P 2 Tp 8P 2 Tp
where
b pÞ ¼ HC~ p ;div ð X
3 o b b p Þ and s n ¼ 0 on C ep e s 2 L2 ð X s 2 L ð X p Þ r
2
and P0 is the space of constant functions. In the following j is the index of the iterations by subdomains. The weak form of the problem for the porous media reads: for all t > 0, find ðqh ; Uh Þ 2 Wph Qph such that 8Z Z Z > jþ1 e j
1 jþ1 b b > > < ^ K qh s h d X p ^ Uh r sh d X p þ k sh np dC ¼ 0 8sh 2 Wph ; C Xp X Z p Z ð9Þ jþ1 oUh > jþ1 b > b e > d X þ u d X ¼ 0 8u 2 Q : r q u S p p ph h h h : 0 ^ h ot ^p Xp X The velocity in the porous media can be calculated a-posteriori when the discharge vector is jþ1 known using the relation Vjþ1 ph ¼ qh =n. As for the 3D-NH-SW model the domain is discretized using a multilayer approach [6] which induces a grid Tf of prism; each prism P can be seen as the cartesian product of a triangle T 2 Txy (where Txy is a subdivision into triangles of X) with a vertical interval dz. The time advancing is achieved using a fractional step scheme with a Lagrangian treatment of the convective terms [4]. To derive the weak form of the 3D-NH-SW model we assume Copen ¼ ; and we introduce the following spaces n o b f ; @z Þ; r ^ b f Þj^ bf ; u 2 H 1ðX u 2 L2 ð X u nxy ¼ 0 on o X Vf ¼ ^ ) ^ o w 2 b f Þ; b f Þjw bf ; ^ 2 L ðX ^ ¼ 0 on o X w 2 L ðX oz
( Wf ¼
2
b f ; @z Þ is the following anisotropic Sobolev space [16] where H 1 ð X 2 ow 2 1 b 2 b 2 b ^ ðx; y; zÞ 2 ðL ð X f ÞÞ H ð X f ; @z Þ ¼ w 2 ðL ð X f ÞÞ : oz
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
79
The discrete n spaces are o uh 2 Vf j^ uhjP 2 RT0 ðT Þ P1 ðdzÞ 8P 2 Tf ; Vfh ¼ ^ n o ^ h 2 Wf jw ^ hjP 2 P0 ðT Þ P1 ðdzÞ 8P 2 Tf ; Wfh ¼ w n o ghjT 2 P0 ðT Þ 8T 2 Txy ; Pfh ¼ g^h 2 L2 ðXÞj^ n o b f Þj^ qhjP 2 P0 ðP Þ 8P 2 Tf ; Qfh ¼ q^h 2 L2 ð X where P1 is the space of linear functions in the z variable. The weak form reads: for all t > 0, find ðufh ; wfh ; gh ; qh Þ 2 Vfh Wfh Pfh Qfh such that 8Z Z > Dujþ1 oujþ1 o^ uh b fh > b > ^ uh d X f þ dXf mv fh > > Dt oz oz ^f ^f > X X > p ffiffi ffi > Z > > mv aBJ 3 jþ1 > >
pffiffiffiffiffiffiffiffi ðufh ujþ1 uh dX > ph Þnz ^ > > tr K X > Z Z > > > jþ1 bf
bf > uh d X qjþ1 uh d X > h r^ > ^ ggh r ^ > ^f Xf X > Z > > > > b f 8^ > ¼ f^ uh d X uh 2 Vf ; > > ^ > Xf > Z
qjþ1 d X f ¼ 0 8w > h > oz > ^f X > > Z Z > > owjþ1 > jþ1 fh b b f ¼ 0 8^ > r ufh q^h d X f þ qh 2 Qf ; q^ d X > > > oz h ^f ^f X X > > Z Z Z > gh > ogjþ1 > h > ^ ^h dX > g dX þ r ujþ1 h fh dz g > > ot > X X h Z > > > > : ¼ ðVjþ1 gh dX 8^ gh 2 Pf ; ph np Þ^ X
jþ1 þ ð1 hÞkj . and finally k ¼ hpnh The weak formulation of the two problems shows that the interface conditions we have proposed are easily incorporated in the models we have used. jþ1
5. Numerical results In this section we present some preliminary numerical results concerning the coupling of free surface and groundwater flows. In this tests we have considered the situation presented in Fig. 2 in which the interface and the bottom of the porous media are flat. In the first simulation we consider a water channel 1 m wide and 10 m long with a water depth of 1 m covering a porous media whose height is 9 m. At the inlet of the channel a discharge of 10 4 m3 /s is prescribed while at outlet a non-reflective boundary condition is enforced. The porous
80
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
b in the case of a flat interface. Fig. 2. Schematic representation of the domain X
media is homogeneous (K ¼ dI, where d ¼ 10 5 and I is the identity tensor) and aBJ ¼ 1. The time step is set to Dt ¼ 600 s. The mesh for the porous media is composed by 15 320 prisms while the mesh for the free surface flows is composed by 13 788 prisms. In Fig. 3 the velocity profile in the center of the channel is shown: the set of boundary conditions (5) using the Beavers and Joseph model for the bottom friction gives rise to a parabolic velocity profile as expected. In the second test we consider the situation shown in Fig. 4 in which a water basin of 20 m side and 10 m depth covers a porous media whose height is 9 m. In the basin the free surface is initially flat and the velocity is null while in the porous media the initial potential head is 0. The mesh for the porous media is composed by 23 706 prisms while the mesh for the free surface flows is
Fig. 3. Vertical profile of the velocity field in the center of the basin (longitudinal section).
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
81
Fig. 4. Setting for the second test case.
Fig. 5. Vertical section of the velocity profile at the beginning of the simulation in the free surface fluid.
composed by 26 340 prisms. Due to the weight of the water a filtration in the porous media takes place: in Figs. 5 and 6 the vertical section of the velocity profile is shown at the beginning of the simulation and after two days, while Fig. 7 shows the piezometric head and the velocity field after t ¼ 48 h (the scales of the velocity vectors in the figure for the free surface fluid and for the porous media are different). In Fig. 8 the solid and the dashed lines show the water volume in the basin as
Fig. 6. Vertical section of the velocity profile after two days in the free surface fluid.
82
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
Fig. 7. Piezometric head and velocity field after two days.
Fig. 8. Volume balance.
a function of time and the net flux of water through the bottom, respectively: it can be appreciated that the total mass of water is conserved. In the numerical results presented in this paper the iteration by subdomain is performed at each time step. Nevertheless the algorithm we have presented can be used also in those situations where the time scales in the porous media and in the free surface fluid are very different (for instance in presence of porous media with very low permeability). In this case the iteration by subdomain is performed only every N time steps of the free surface fluid and N is such that Dtp ¼ N Dtf , where Dtp and Dtf are the time steps for the porous media and for the free surface fluid, respectively. Remark 1. The numerical results show that the iterative scheme we have proposed is convergent: in particular, we have taken a tolerance ¼ 10 3 and the convergence is obtained in about five iterations. From the theoretical point of view the convergence analysis of this procedure will rely on the formulation of a Steklov–Poincare problem for the unknown U at the interface C and will be the subject of future work.
E. Miglio et al. / Computers & Fluids 32 (2003) 73–83
83
6. Conclusion In this paper we have addressed the problem of coupling free surface and groundwater flows. A set of boundary conditions based on physical considerations is proposed; to model the bottom stress in the free surface fluid the Beavers and Joseph formula is employed. A suitable iterative scheme to decouple the solution of the full problem is proposed. The numerical tests show that the iterative scheme is convergent and the effectiveness of the interface conditions we have used.
Acknowledgement We thank Eng. T. Ponta for having performed some simulations.
References [1] Bayada G, Chambat M. On the interface conditions for a thin film flow past a porous medium. SIAM J Math Anal 1995;26:1113–29. [2] Bear J. Hydraulics of groundwater. New York: McGraw-Hill; 1979. [3] Beavers GS, Joseph DD. Boundary conditions at a naturally permeable wall. J Fluid Mech 1967;30:197–207. [4] Causin P, Miglio E, Saleri F. Algebraic factorizations for 3D non-hydrostatic free surface flows. Accepted for publication in Comput Visual Sci. quations et phenomenes de surface pour l’ecoulement dans un modele de milieu [5] Ene HI, Sanchez-Palencia E. E poreaux. J Mecanique 1975;14:73–108. [6] Fontana L, Miglio E, Quarteroni A, Saleri F. A finite element method for 3D hydrostatic water flows. Comput Visual Sci 1999;2(2/3):85–93. [7] J€ager W, Mikelic A. On the boundary conditions at the contact interface between a porous medium and a free fluid. Ann Scuola Norm Sup Pisa Cl Sci 1996;23:403–65. [8] J€ager W, Mikelic A. On the interface boundary condition of Beavers, Joseph and Saffman. SIAM J Appl Math 2000;60:1111–27. [9] Jones IP. Low Reynolds number flow past a porous spherical shell. Proc Cambridge Philos Soc 1973;73:231–8. [10] Nield DA, Bejan A. Convection in porous media. New York: Springer; 1999. [11] Ochoa-Tapia JA, Whitaker S. Momentum transfer at the boundary between a porous medium and a homogeneous fluid. I. Theoretical development. Int J Heat Mass Transfer 1995;38:2635–46. [12] Payne LE, Straughan B. Analysis of the boundary condition at the interface between a viscous fluid and a porous medium and related modelling questions. J Math Pures Appl 1998;77:317–54. [13] Quarteroni A, Valli A. Domain decomposition method for partial differential equations. Oxford: Oxford Science Publications; 1999. [14] Raviart PA, Thomas JM. A mixed finite element method for 2nd order elliptic problems. In: Dold A, Eckmann B, editors. Mathematical aspects of finite elements methods. Lecture Notes 606. Berlin, Heidelberg, New York: Springer; 1977. [15] Thomee V. Galerkin finite elements methods for parabolic problems. Berlin, Heidelberg: Springer; 1997. [16] Temam R. Sur la stabilite et la convergence de la methode des pas fractionnaires. Ann Mat Pura et Appl IV Ser 1968;79:191–379. [17] Wood WL. Introduction to numerical methods for water resources. Oxford: Oxford Science Publications; 1993.