On Mixed Isogeometric Analysis of Poroelasticity

On Mixed Isogeometric Analysis of Poroelasticity Yared W. Bekele∗ 1 , Eivind Fonn2 , Trond Kvamsdal2,3 , Arne M. Kvarving2 and Steinar Nordal4 1 Rock and Soil Mechanics Group, SINTEF Building and Infrastructure, Trondheim, Norway 2 Department of Applied Mathematics, SINTEF ICT, Trondheim, Norway 3 arXiv:1706.01275v1 [math.NA] 5 Jun 2017 4 Department of Mathematical Sciences, NTNU, Trondheim, Norway Department of Civil and Environmental Engineering, NTNU, Trondheim, Norway Abstract Pressure oscillations at small time steps have been known to be an issue in poroelasticity simulations. A review of proposed approaches to overcome this problem is presented. Critical time steps are specified to alleviate this in finite element analyses. We present a mixed isogeometric formulation here with a view to assessing the results at very small time steps. Numerical studies are performed on Terzaghi’s problem and consolidation of a layered porous medium with a very low permeability layer for varying polynomial degrees, continuities across knot spans and spatial discretizations. Comparisons are made with equal order simulations. 1 Introduction The study of porous materials, where the flow of fluid and solid deformation are coupled, is essential in several areas of science and engineering. The theory of poroelasticity is a mathematical formulation developed to describe these coupled processes and predict the response of fluid saturated/unsaturated porous media to external loading. There are different types of porous materials that are studied under this theory such as soil, rock, concrete and other man-made materials. Poroelasticity has a wide range of applications in different disciplines of engineering mechanics and natural sciences. Some of the application areas include geomechanics, biomechanics, reservoir engineering and earthquake engineering. In addition to these diverse areas of application, it is gaining popularity in the study of modern man-made porous media in material science. The mathematical formulations describing the fluid-solid coupled processes are developed based on porous media theory where the multiphase medium is approximated as a continuum, [7]. The volume fraction concept is used for averaging the properties of the multiphase medium in a continuum formulation. The governing partial differential equations of poroelasticity were first developed for a onedimensional case by Terzaghi [33, 34]. The formulations were later generalized for a threedimensional case and extended by Biot [4, 2, 3]. The mathematical formulations have been studied extensively by several researchers since then. Various analytical and numerical studies have been proposed in the literature. Analytical solutions were obtained for problems with simplified material domains and boundary conditions. Application to boundary value problems with complex material domains and boundary conditions required the use of numerical methods. The emergence of the finite element method opened the door for a detailed numerical study of poroelasticity and for application to arbitrary geometries and boundary conditions. The finite element method was first applied to the governing equations of poroelasticity to solve the initial boundary value problem of flow in a saturated porous elastic medium by Sandhu and Wilson [29]. Hwang et al. [16] also used the finite element method for plane strain consolidation problems and verified the results against closed form solutions. The application of the finite element method started gaining momentum afterwards and several researchers engaged themselves ∗ Corresponding author. Email addresses: [email protected] (Yared W. Bekele), [email protected] (Eivind Fonn), [email protected] (Trond Kvamsdal), [email protected] (Arne M. Kvarving), [email protected] (Steinar Nordal) 1 not only on application problems but also in the investigation of the numerical properties of the method within the context of poroelasticity. Ghaboussi and Wilson [13] applied the finite element method to partially saturated elastic porous media and first noticed the ill-conditioning of the matrix equations that may result when an incompressible fluid is assumed to occupy the pore spaces. Booker and Small [5] investigated the stability of the numerical solution when the finite element method is applied to Biot’s consolidation equations. The stability was studied for different numerical integration schemes and time-step sizes. The numerical performance of some finite element schemes for analysis of seepage in porous elastic media was studied by Sandhu et al. [28]. They studied various spatial and temporal discretization schemes and evaluated the numerical performances against the analytical solution of Terzaghi’s one-dimensional consolidation problem. Triangular and quadrilateral elements with equal and mixed orders of interpolation for the displacement and pressure were considered. It was shown that the elements with equal orders of interpolation showed oscillatory behavior in the solution. Vermeer and Verruijt [35] derived a lower bound for the time-step size in the analysis of consolidation by finite elements in terms of the mesh size and the coefficient of consolidation. They showed that there is an accuracy condition in the finite element analysis of consolidation by using a critical time-step, below which oscillatory solutions are observed. The derived critical time-step is strictly valid for a one-dimensional case and a uniform finite element mesh. Reed [26] analyzed the numerical errors in the analysis of consolidation by finite elements. It was shown that the use of a mixed formulation for the field variables helps in reducing the pore pressure oscillations but may not remove them entirely. They instead used Gauss point smoothing to eliminate the pore pressure oscillations. Special finite elements for the analysis of consolidation were proposed by Sandhu et al. [27]. They presented “singularity” elements to model pore pressures in the vicinity of free-draining loaded surfaces immediately after application of loads. The elements were special in that they use special interpolation schemes which reflect the actual variation of the field variables. The finite element method became a well-established method for the analysis of poroelasticity problems and the mathematical properties of the governing equations and the numerical solution were studied in a further great detail. Murad and Loula [21] presented numerical analysis and error estimates of finite element approximations of Biot’s consolidation problem. They used a mixed formulation and improved the rates of convergence by using a sequential Galerkin Petrov-Galerkin post-processing technique. In a further study, [22], they investigated the stability and convergence of finite elements approximations of poroelasticity. They derived decay functions showing that the pore pressure oscillations, arising from an unstable approximation of the incompressibility constraint on the initial conditions, decay in time. Finite element analysis of consolidation with automatic time-stepping and error control was presented by Sloan and Abbo [30, 31]. Automatic time increments were selected such that the temporal discretization error in the displacements is close to a specified tolerance. Ferronato et al. [11] studied the ill-conditioning of finite element poroelasticity equations with a focus on the instabilities that may affect the pore pressure solution. They claim that the origin of most instabilities is due to the assumption that, for initial conditions, the porous medium behaves as an incompressible medium if the pore fluid is incompressible. They also argue that oscillatory pore pressure solutions may not always be observed for very stiff and low permeable materials depending on the critical time step. Gambolati et al. [12] studied the numerical performance of projection methods in finite element consolidation models. Dureisseix et al. [8] proposed a large time increment (LATIN) computational strategy for problems of poroelasticity to improve the efficiency of the finite element analysis. A finite element formulation to overcome spatial pore pressure oscillations caused by small time increments was proposed by Zhu et al. [36]. Korsawe et al. [18] compared standard and mixed finite element methods for poroelasticity. In particular, Galerkin and least-squares mixed finite element methods were compared. They claim that Galerkin’s method is able to preserve steep pressure gradients but overestimates the effective stresses. On the other hand, a least-squares mixed method was noticed to have the advantage of direct approximation of the primary variables and explicit approximation of Neumann type boundary conditions but to be computationally more expensive. A mixed least-squares finite element method for poroelasticity was also proposed by Tchonkova et al. [32], claiming that pore pressure oscillations are eliminated for different temporal discretizations. A coupling of mixed and continuous Galerkin finite element methods for poroelasticity was investigated for continuous and discrete in time cases by Phillips and Wheeler [23, 24]. They also studied a coupling of mixed and discontinuous Galerkin finite-element methods, [25]. Haga et al. [14] studied the causes of pressure oscillations in low-permeable and low-compressible media by presenting two, three and four field 2 mixed formulations in terms of the field variables displacement, pore fluid pressure, fluid velocity and solid skeleton stress. A posteriori error estimation and adaptive refinement in poroelasticity has been studied by very few researchers. Larsson and Runesson [19] presented a novel approach for space-time adaptive finite element analysis for the coupled consolidation problem in geomechanics. El-Hamalawi and Bolton [9] proposed an a posteriori error estimator for plane-strain geotechnical analyses based on superconvergent patch recovery with application to Biot’s consolidation problem. They later extended the application of the a posteriori estimator for axisymmetric geotechnical analyses in [10]. Adaptive isogeometric finite element analysis, with LR B-Splines, for steady-state groundwater flow problems was presented by Bekele et al. [1] but application to poroelasticity still remains as a task to study. Isogeometric finite element analysis of poroelasticity was first presented by Irzal et al. [17]. The advantages of the smoothness of the basis functions in isogeometric analysis were highlighted in their application. One of the advantages of higher continuity is that the numerical implementation results in a locally mass conserving flow between knotspans, analogous to elements in finite element analysis. But the formulation presented relied on equal orders of interpolation for the field variables in poroelasticity, namely displacement and pore fluid pressure. Such a formulation, while still useful for several applications without significant numerical challenges, has limitations when it comes to problems where the material properties or boundary conditions are problematic. In this paper, we present a mixed isogeometric formulation for poroelasticity. To our best knowledge, this is the first time that a mixed isogeometric formulation for poroelasticity is presented. The paper is structured as follows. In Section 2, the governing equations of poroelasticity are presented. The fundamentals of isogeometric analysis and its particular features of interest within the current context are discussed in Section 3. Numerical examples are given in Section 4 and the observations are summarized with concluding remarks in Section 5. 2 Governing Equations Biot’s poroelasticity theory [4, 2] couples elastic solid deformation with fluid flow in the porous medium where the fluid flow is assumed to be governed by Darcy’s law. The governing equations of the theory, the necessary boundary conditions, weak formulation and Galerkin finite element discretization are presented in the following sections. 2.1 Linear Momentum Balance Equation The linear momentum balance equation for a fluid-saturated porous medium is given by:
∇ · σ ′ + αpf I +ρb = 0 | {z } (1) =σ where σ is the total stress, σ is the effective stress, α is Biot’s coefficient, pf is the fluid pressure, I is an identity matrix, ρ is the overall density of the porous medium and b represents body forces. The Biot coefficient α can be calculated from: ′ α=1− Kt Ks (2) where Kt and Ks are the bulk moduli of the porous medium and solid particles, respectively. The constitutive equation for poroelasticity relates stress and strain linearly as: σ′ = D : ε (3) where D is a fourth-order stiffness tensor. Small deformations are also assumed, so the strain ε satisfies a linear first-order equation with respect to the displacement u, ε= 1 (∇u + ∇⊺ u) 2 where 1/2(∇ + ∇⊺ ) is the symmetrized gradient operator i.e.   1 ∂ui ∂uj εij = . + 2 ∂xj ∂xi 3 (4) (5) In the following, it will be convenient to lower tensors and higher differential operators to Voigt notation, which represents the symmetric d × d tensor σ ′ as a d(d+1)/2-vector, which we will denote with a tilde:   ′ ′ ′ σxx σxy σxz  ′ ⊺ ′ ′ ′  ′ ′ ′ ′ ′ σxy σyz σyy σyy σzz σyz σxz σxy ⇐⇒ σxx {z } | ′ ′ ′ σxz σyz σzz σ̃ ′ | {z } σ′ A similar conversion takes place for the strains, where the shear strains are replaced by the engineering shear strains:   εxx εxy εxz  εxy εyy εyz  ⇐⇒ εxx εyy εzz 2εyz 2εxz 2εxy ⊺ | {z } εxz εyz εzz ε̃ | {z } ε Voigt notation allows us to express the equilibrium equation and the stress-strain equation using the same differential operator L,   ∂ ∂ ∂ 0 0 ∂y 0 ∂z ∂x  ∂ ∂ ∂ 0 ∂x 0 (6) L⊺ =  0 ∂y  ∂z ∂ ∂ ∂ 0 0 ∂z 0 ∂y ∂x Using L yields the following equilibrium equation in terms of the two primary unknowns u and pf , L⊺ D̃Lu − α∇pf + ρb = 0, (7) where D̃ is the Voigt notation equivalent of D, taking into account the aforementioned engineering shear strains. We will generally assume isotropic materials, where D̃ takes the block form (in terms of Young’s modulus E and Poisson’s ratio ν)   E D̃11 0 (8) D̃ = D̃22 (1 + ν)(1 − 2ν) 0 where the two blocks are given as D̃11 = (1 − 2ν)I + ν1 1 − 2ν D̃22 = I 2 (9) and 1 is a matrix of ones. 2.2 Mass Balance Equation A mass conservation equation together with the equilibrium equation in (7) completes the governing equations of poroelasticity. The fluid content ζ is given by ζ = α∇ · u + cpf where c is the storativity or specific storage coefficient at constant strain. It is given by n α−n + c= Ks Kf (10) (11) where Kf is the bulk modulus of the fluid and n is the porosity of the material. The change in the fluid content ζ satisfies the equation ∂ζ +∇·w =0 (12) ∂t where w is the fluid flux, which is given by Darcy’s law as:
1 w = − k · ∇pf − ρf b (13) γf where γf is the unit weight of the fluid, ρf its density and k is the hydraulic conductivity matrix. The final equation of mass balance is then  
∂pf 1 α∇ · u̇ + c (14) + ∇ · − k · ∇pf − ρf b = 0. ∂t γf 4 2.3 Boundary Conditions The governing linear momentum and mass balance equations in (7) and (14), respectively, are accompanied by the usual boundary conditions in the formulation of bounary value problems. Let (ΓuD , ΓpD ) and (ΓuN , ΓpN ) be two partitions of the boundary ∂Ω of domain Ω, for representing Dirichlet and Neumann boundary conditions, respectively. The Dirichlet boundary conditions for the equilibrium (7) and mass balance (14) equations are ( u = u on ΓuD , (15) pf = pf on ΓpD , where u and pf are the prescribed displacement and pressure, respectively. The Neumann boundary conditions are ( σ · n = t on ΓuN , w · n = q on ΓpN , (16) where n is the outward pointing normal vector, t is the surface traction and q is the fluid flux on the boundary. 2.4 Variational Formulation To derive the variational formulations of equations (7) and (14), we introduce a vector-valued test function δu, which vanishes on ΓuD , and a scalar test function δp, which vanishes on ΓpD . We start with the total stress formulation of the linear momentum balance equation, which from equation (7) is given by ∇ · σ + ρb = 0. (17) Multiplying by the test function δu and integrating over the domain Ω gives Z Z δu⊺ ∇ · σdΩ + δu⊺ ρbdΩ = 0. Ω (18) Ω The first term in the above equation contains a double derivative of the unknown displacement, and is relaxed using a form of Green’s theorem, Z XZ δu⊺ ∇ · σdΩ = δui ∇ · σi dΩ Ω Ω i = XZ = δui σi · ndΓ − ∂Ω i Z XZ i δu⊺ tdΓ − Γu N Z ∇δui · σi dΩ Ω (19) ∇δu : σdΩ. Ω Due to the symmetry of the stress tensor, the last term is expressible in Voigt notation, ⊺ ∇δu : σ = (Lδu) σ, (20) yielding the weak form of (7) as Z Z Z Z ⊺ ⊺˜ f (Lδu) D̃(Lu)dΩ − α (Lδu) Ip dΩ = δu⊺ ρbdΩ + Ω Ω Ω δu⊺ tdΓ (21) Γu N where we have used I˜ as the Voigt notation identity operator, which for a general three-dimensional case is given by: ⊺ I˜ = {1, 1, 1, 0, 0, 0} (22) For the mass balance equation, multiplying (14) by the scalar test function δp and integrating over the domain Ω, we get   Z Z Z
1 ∂pf (23) dΩ + δp∇ · − k · ∇pf − ρf b dΩ = 0. α δp∇ · u̇dΩ + c δp ∂t γf Ω Ω Ω 5 Again, by applying Green’s theorem to the last term, we obtain     Z Z Z

1 1 f f ∇δp · − k · ∇p − ρf b dΩ. δp∇ · − k · ∇p − ρf b dΩ = δpqdΓ − γf γf Ω Ω Γp N Thus, the weak form of the mass balance equation, (14), is Z Z Z Z Z 1 1 ∂pf δpqdΓ. ∇δp⊺ kρf bdΩ − ∇δp⊺ k∇pf dΩ = dΩ + δp δp∇ · u̇dΩ + c α ∂t γf γf Ω Γp Ω Ω Ω N 2.5 (24) (25) Galerkin Finite Element Formulation With a suitable number N of basis functions defined, let Np : Ω → R1×N and Nu : Ω → Rd×dN be the basis interpolation matrices for the pressure and displacement respectively. The unknowns and the test functions can then be represented using coefficient vectors: u = Nu uc , δu = Nu δuc , p f = Np p c , δp = Np δpc (26) where uc and pc are the control point values of the displacement and pressure field variables. Application of (26) to the weak form of the linear momentum balance equation in (21) results in the matrix the discrete system of equations (after canceling δuc and δpc , as equations (7) and (14) are supposed to be valid for any choice of these) Kuc − Qpc = fu (27) where the stiffness matrix K, the coupling matrix Q and the vector of body forces and surface tractions fu are given by Z K= B ⊺ D̃BdΩ Ω Z ˜ p dΩ B ⊺ αIN Q= (28) Ω Z Z Nu⊺ ρbdΩ + fu = Nu⊺ tdΓ Γu N Ω Here B = LNu is the strain-displacement matrix. Similarly, using (26) in the weak form of the mass balance equation in (25) results in the discrete system of equations ∂uc ∂pc +S + P pc = fp (29) ∂t ∂t where the storage matrix S, the permeability matrix P and the vector of fluid body forces and fluxes fp are given by Z Np⊺ cNp dΩ S= Ω Z 1 ∇Np⊺ k∇Np dΩ P = (30) γf Ω Z Z 1 fp = ∇Np⊺ kρf bdΩ − Np⊺ qdΓ. γf Ω Γp N Q⊺ Combining equations (27) and (29) results in the coupled system of equations for poroelasticity    c      0 0 u̇ K −Q uc fu . (31) + = Q⊺ S ṗc 0 P fp pc A symmetric system of equations can be obtained by time-differentiating the first equation and multiplying one of the equations by −1, [20]:        c  −K Q u̇c 0 0 u −f˙u + (32) = Q⊺ S ṗc 0 P pc fp In this formulation, it is important that time-dependent quantities involved in fu , such as traction and body forces, are “ramped up” from an initial equilibrium instead of being applied immediately. This can be done in the first time step. 6 2.6 Temporal Discretization The generalized trapezoidal rule (GTR) is applied for the temporal discretization of the coupled system of matrix equations in (32). Representing the vector of unknowns by X = {uc , pc }⊺ , we have the GTR approximation ∂X ∂t = n+θ Xn+1 − Xn ∆t (33) Xn+θ = (1 − θ)Xn + θXn+1 where θ is a time integration parameter which has limits 0 ≤ θ ≤ 1 and n is a time step identifier. Adopting backward Euler time stepping (θ = 1) with time step ∆t and applying (33) to (32) we obtain the system of equations    c      −K Q u −K Q uc −f˙u = + ∆t (34) Q⊺ S + ∆tP pc n+1 Q⊺ S pc n fp n+1 which is a linear system in this case, for poroelasticity, as the coefficient matrices are independent of the unknowns. 3 3.1 Isogeometric Analysis Introduction Since its first introduction by Hughes et al. [15], isogeometric analysis (IGA) has been successfully applied to several areas of engineering mechanics problems. The fundamental aim for the introduction of IGA was the idea of bridging the gap between finite element analysis (FEA) and computer-aided design (CAD). The main concept behind the method is the application of the same basis functions used in CAD for performing finite element analysis. In the process of its application to various engineering problems, IGA has shown advantages over the conventional finite element method, for instance the ease of performing finite element analysis using higher order polynomials. We briefly present the fundamentals behind B-Splines and Non-Uniform Rational B-Splines (NURBS) in the next section and highlight the features of IGA that are important in our context. 3.2 Fundamentals on B-Splines and NURBS We start the discussion on B-Splines and NURBS by first defining a knot vector. A knot vector in one dimension is a non-decreasing set of coordinates in the parameter space, written Ξ = {ξ1 , ξ2 , ..., ξn+p+1 }, where ξi ∈ R is the ith knot, i is the knot index, i = 1, 2, ..., n + p + 1, p is the polynomial order, and n is the number of basis functions. Knot vectors may be uniform or non-uniform depending on whether the knots are equally spaced in the parameter space or not. A univariate B-Spline curve is parametrized by a linear combination of n B-Spline basis functions, {Ni,p }ni=1 . The coefficients corresponding to these functions, {Xi }ni=1 , are referred to as control points. The B-Spline basis functions are recursively defined starting with piecewise constants (p = 0): ( 1 if ξi ≤ ξ < ξi+1 (35) Ni,0 (ξ) = 0 otherwise For higher-order polynomial degrees (p ≥ 1), the basis functions are defined by the Cox-de Boor recursion formula: Ni,p (ξ) = ξi+p+1 − ξ ξ − ξi Ni,p−1 (ξ) + Ni+1,p−1 (ξ) ξi+p − ξi ξi+p+1 − ξi+1 (36) B-Spline geometries, curves, surfaces and solids, are constructed from a linear combination of B-Spline basis functions. Given n basis functions Ni,p and corresponding control points Pi ∈ Rd , i = 1, 2, ..., n, a piecewise polynomial B-Spline curve is given by: C(ξ) = n X i=1 7 Ni,p (ξ)Pi (37) C −1 C3 C2 C1 C0 C −1 1 22 333 4444 55555 1 00000 Figure 1: Different continuities across knotspans, after [6]. Similarly, for a given control net Pi,j , i = 1, 2, ..., n, j = 1, 2, ..., m, polynomial orders p and q, and knot vectors Ξ = {ξ1 , ξ2 , ..., ξn+p+1 }, and H = {η1 , η2 , ..., ηm+q+1 }, a tensor product B-Spline surface is defined by: S(ξ, η) = m n X X Ni,p (ξ)Mj,q (η)Pi,j (38) i=1 j=1 B-Spline solids are defined in a similar way as B-Spline surfaces from tensor products over a control lattice. NURBS are built from B-Splines to represent a wide array of objects that cannot be exactly represented by polynomials. A NURBS entity in Rd is obtained by projective transformation of a B-Spline entity in Rd+1 . The control points for the NURBS geometry are found by performing exactly the same projective transformation to the control points of the B-Spline curve. A detailed treatment of B-Splines and NURBS can be referred from Cottrell et al. [6]. 3.3 Important Features in Current Context IGA has a number of advantages over FEA such as the ability to represent exact geometries of structures or domains, non-negative basis functions and isoparametric mapping at patch level. In the context of the current work, we focus on the features of IGA that are especially important. These features are improved continuity because of the smoothness of the basis functions and the ability to perform simulations with high continuity and high regularity meshes. We look closely into each here. 3.3.1 Continuity One of the most distinctive and powerful features of IGA is that the basis functions will be C p−m continuous across knotspans (analogous to elements in FEA), where p is the polynomial degree and m is the multiplicity of the knot. This means that the continuity across knotspans can be controlled by the proper choice of p and m. The continuity can be decreased by repeating a knot - important to model non-smooth geometry features or to facilitate the application of boundary conditions. For instance, quadratic (p = 2) splines are C 1 continuous over nonrepeated knots while quadratic Lagrange finite element bases are only C 0 continuous. If we consider the quartic (p = 4) basis functions constructed from the open, non-uniform knot vector Ξ = {0, 0, 0, 0, 0, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5}, we get different continuities across knotspans as shown in Figure 1. 3.3.2 k-refinement IGA and FEA both allow h- and p-refinements i.e. increasing the number of knotspans by knot insertion (increasing the number of elements in FEA) and raising the polynomial order. The noncommutativity of knot insertion and polynomial order elevation results in a type of refinement that is unique to IGA, called k-refinement. This is achieved by performing polynomial order elevation followed by knot insertion. This results in a high continuity mesh with the least number of degrees of freedom i.e. high regularity. 8 p=1 p=2 p=3 Figure 2: Number of control points for a given element on a simple B-Spline surface with different polynomial degrees. The element is highlighted and the blue squares represent control points. 3.4 Mixed Isogeometric Formulation A mixed formulation is constructed by first defining the knot vectors and basis functions defining the geometry of the domain. The polynomial order defining the geometry is used as the polynomial degree for one of the field variables and is raised by the desired degree for the other field variable. In our context, the polynomial order for the pressure, pp , is defined by the geometry construction and the polynomial order for the displacement, pu , is raised by one. Both pp and pu can then be raised to the desired degree starting from the initial definition. For example, a simple two-dimensional geometry defined by the knot vectors Ξ = {0, 0, 1, 1} and H = {0, 0, 1, 1} implies pp = 1 and pu = 2 with 4 and 9 control points, respectively. The number of control points, location of degrees of freedom in IGA, on a B-Spline surface for different polynomial degrees is shown in Figure 2. 4 Numerical Examples In this section, the performance of a mixed isogeometric formulation is investigated for some numerical examples. We first consider Terzaghi’s classical one-dimensional consolidation problem for verification and mesh convergence studies. Consolidation of a layered medium with a low permeability layer sandwiched between two high permeability layers is studied. The mixed formulation results are compared with equal order simulation. 4.1 Terzaghi’s Problem Terzaghi’s problem is a classical one-dimensional consolidation problem with an analytical solution, which makes it suitable for code validation. A saturated porous medium subjected an external loading under plane-strain condition is considered where the fluid is allowed to dissipate only at the top boundary, hence resulting in a one-dimensional consolidation. A no flux boundary condition is assumed for the lateral and bottom boundaries. The displacement boundary conditions are such that the lateral sides are constrained from horizontal deformation and the bottom boundary is fixed in both the horizontal and vertical directions. The external load is applied as a Neumann traction p0 at the top boundary. The domain and boundary conditions considered are shown in Figure 3. The analytical solution for the pressure field as a function of time and space is given by:   ∞ h 2 4 X (−1)i−1 πy i pf (t, y) 2 π ts = exp −(2i − 1) cos (2i − 1) p0 π i=1 2i − 1 4 2h where the dimensionless time ts is given as a function of the consolidation coefficient cv and drainage path h (total height for one-way drainage) by: ts = cv t. h2 The consolidation coefficient cv is given by: 9 (39) ty = −p0 , ux = 0, pf = 0 ux = 0 ux = 0 h = 8 mm q=0 q=0 u = 0, q = 0 Figure 3: Terzaghi’s problem: Domain and boundary conditions. Table 1: Terzaghi’s problem: Load and material parameters. Parameter Value External load, p0 Hydraulic conductivity, k Biot’s coefficient, α Young’s modulus, E Poisson’s ratio, ν Storativity, c Body forces, b cv = Unit 6 1.0 × 10 1.962 × 10−14 1.0 6.0 × 106 0.4 0 0 (1 − ν)Eκ (1 + ν)(1 − 2ν) Pa m2 − Pa − Pa−1 N (40) The material parameters used for this problem are given in Table 1, as used in [17]. The choice of the storativity value c = 0 effectively corresponds to assuming incompressible solid grains and an incompressible fluid. The Terzaghi verification problem is simulated in a mixed and equal order formulation for comparison. The polynomial degrees considered for the pressure are pp = 1, 2, 3. The corresponding values for the displacement in a mixed formulation are pu = 2, 3, 4. The number of elements used in the simulation is Ne = 72. Critical and sub-critical time step sizes are considered to study the sensitivity of the simulations to temporal discretization and to evaluate accuracy of the solution for small time step sizes. The critical time step is calculated according to the relation derived in [35]. The results from a simulation using the critical time step are shown in Figure 4. A linear solution space is used for the pressure and a quadratic space for the displacement. The results from simulations with a sub-critical time step are shown in Figure 5 for mixed and equal order cases. The results with the time step size equal to the critical time step show no oscillations in the pressure values. On the other hand, slight oscillations are visible for the sub-critical time step case. These oscillations at very small time steps appear worse for the equal order simulations compared to the mixed simulation. In both cases, the results are observed to improve with increasing polynomial degrees. 10 1 t t t t t pf/p0 0 0 = = = = = 2∆tc 100∆tc 500∆tc 2000∆tc 5000∆tc 1 y/h Figure 4: Numerical solution to the Terzhagi problem with pp = 1, pu = 2 and and Ne = 72 using critical time step. pp = 1 pf/p0 pp = 2 1 0 pf/p0 0 1 pp = 3 1 0 0 y/h 1 y/h (a) Mixed (b) Equal order Figure 5: Numerical solution to the Terzaghi problem with Ne = 72 using a sub-critical time step of ∆t = 0.1∆tc for different polynomial degrees. All plots are shown for the first time step. 4.2 Terzaghi’s Problem: Convergence Study Next, a simplified version of the Terzaghi problem is used as a convergence study. We consider a domain with dimensions of w × h = 1 × 1 with the same boundary conditions as in the previous case. For simplicity we choose the following material parameters: α = 1, c = 0, E = 2/3, ν = 0.25 and κ = 1. The external load applied is p0 = 1 and we assume no body forces i.e. b = 0. This case was run with an increasing number of degrees of freedom using polynomial degrees pp = 1, 2, 3 for the pressure and correspondingly pu = 2, 3, 4 for the displacement. In all cases, the time step was kept sufficiently small for the spatial discretization error to dominate and we look at the results at the end of the first time step. The convergence study is performed by calculating the relative L2 error of the pressure field. The relative error based on the computed pressure values, ρh , is calculated from ρh = kpfh − pf kL2 kpf kL2 (41) where pfh and pf are the computed and analytical solution pressures, respectively. The results from the mesh convergence study are shown in Figure 6 in terms of plots of the relative error versus the total number of degrees of freedom. The expected convergence rate based on the analytical solution is also shown. We observe from the results that optimal convergence rates are obtained for all polynomial degrees considered. 11 Relative Error, ρh 10−2 pp = 1 O(N −1 ) pp = 2 10−3 O(N −3/2 ) pp = 3 O(N −2 ) 10−4 10−5 10−6 103 104 Number of DOFs 105 Figure 6: Convergence rates in the relative L2 norm of pressure, for three different polynomial degrees. 4.3 Low Permeability Layer The next example we consider is the consolidation of a very low permeability layer sandwiched between two high permeability layers, as presented in [14]. A one-dimensional consolidation is assumed by applying the appropriate boundary conditions. The fluid is allowed to dissipate at the top boundary and a no flux condition is defined at the lateral and bottom boundaries. The bottom boundary is fixed from vertical and horizontal displacement and the domain is allowed to deform only in the vertical direction. An external load p0 is applied at the top boundary. The problem setup with the boundary conditions is shown in Figure 7. The material parameters for this problem are given in Table 2. Simplified material properties are assumed to focus on the permeability differences of the middle and the bounding layers. The low permeability layer problem is studied using mixed and equal order simulations. The polynomial degrees for the pressure are increased continuously from linear to quartic i.e. pp = 1, 2, 3, 4. The corresponding polynomial degrees for the displacement in a mixed formulation are pu = 2, 3, 4, 5. The continuities at the boundaries between the layers are also varied. We consider C 0 and C pp −1 continuities at these interfaces. In addition, simulations are performed for uniform and graded meshes. The results are presented for these different combinations. The results from simulations with a uniformly refined mesh are shown in Figure 8 for the mixed and equal order cases. Severe pressure oscillations are observed within the low permeability layer for the equal order simulations. Due to its high permeability, the fluid in the top layer dissipates very quickly for the time step size considered here i.e. ∆t = 1s. The pressure oscillations start ty = −p0 , ux = 0, pf = 0 κ1 ux = 0 q=0 κ2 ux = 0 q=0 κ1 u = 0, q = 0 Figure 7: The Haga problem: Domain and boundary conditions. 12 Table 2: The Haga problem: Load and material parameters. Parameter Value External load, p0 Darcy coefficient, k1 /γf Darcy coefficient, k2 /γf Biot’s coefficient, α Young’s modulus, E Poisson’s ratio, ν Storativity, c Body forces, b 1.0 1.0 1.0 × 10−8 1.0 0.67 0.25 0 0 2 Unit Pa m2 /Pa s m2 /Pa s − Pa − Pa−1 N 2 pf/p0 pf/p0 1 0 1 0.25 0.5 0 0.75 0.25 0.5 y/h 0.75 y/h 2 2 pf/p0 pf/p0 1 0 1 0.25 0.5 0 0.75 0.25 0.5 y/h pp = 1 0.75 y/h pp = 2 pp = 3 pp = 4 Figure 8: Numerical solution to the Haga problem using Ne = 60 uniform elements and ∆t = 1 s. All figures are shown after two time steps. On the left the mixed order method, and on the right the equal order method. The continuity in the boundary layer is C pp −1 in the top row, and C 0 in the bottom row. as soon as the fluid in the low permeability layer starts dissipating. The results improve with increasing polynomial degrees but some oscillations are still seen for a quartic solution space for the pressure, pp = 4. The results with C 0 continuities at the material interfaces improve slightly better than with C pp −1 continuity since a C 0 continuity is a more accurate representation of material interfaces. The pressure oscillations in the mixed simulations are less severe and are localized at the boundary between the low permeability and bottom layers. These again decrease with increasing polynomial degrees and a C 0 continuity at the material interfaces. Simulations with a graded mesh are also performed for the different combination of polynomial degrees and interface continuities. The graded mesh is generated such that more elements are concentrated at the material interfaces. The results from this case are shown in Figure 9. The pressure oscillations in the equal order case improve significantly in this case compared to the results from uniform mesh refinement. However, the oscillations still occur throughout the low permeability layer. The equal order results for linear basis functions show a slightly strange behavior in that the oscillations are lesser within the low permeability layer than for higher order elements, but show slightly higher oscillations at the top material interface. The results are again better with a C 0 continuity at the material interfaces. The mixed simulation results also improve with a graded mesh. Almost no oscillations are noticed for combinations of higher polynomial degrees and C 0 continuities at the material interfaces. 13 2 2 pf/p0 pf/p0 1 0 1 0.25 0.5 0 0.75 0.25 0.5 y/h 0.75 y/h 2 2 pf/p0 pf/p0 1 0 1 0.25 0.5 0 0.75 0.25 0.5 y/h pp = 1 0.75 y/h pp = 2 pp = 3 pp = 4 Figure 9: Numerical solution to the Haga problem using a graded mesh with small elements near the boundary layer and ∆t = 1 s. All figures are shown after two time steps. On the left the mixed order method, and on the right the equal order method. The continuity in the boundary layer is C pp −1 in the top row, and C 0 in the bottom row. 5 Conclusions Mixed isogeometric analysis of poroelasticity is presented where different order of polynomials are used for the displacement and pore pressure field variables. Numerical studies on Terzaghi’s classical one-dimensional consolidation problem and consolidation of a layered soil with a middle low permeability layer are presented. The results from mixed polynomial order simulations are compared with equal order analyses. For Terzaghi’s one-dimensional consolidation problem, the pore pressure oscillations are investigated when a time step size less than the critical value is used. The oscillations were observed to be higher in the equal order simulations compared to the mixed order results. The oscillations are not completely removed in the mixed isogeometric simulations but it is observed that the they tend to decrease with increasing polynomial orders for the pore pressure. This is illustrated by the convergence of the relative L2 norm of the pore pressure error for varying polynomial orders. The low permeability layer problem showed similar trends in the pore pressure oscillations i.e. the equal order simulations resulted in worse pore pressure oscillations compared to the mixed results. Again, in both cases, the oscillations decreased with increasing polynomial orders. The use of a graded mesh, where the knot spans are concentrated at the interfaces between the low permeability and other layers, resulted in much lower oscillations both in the equal order and mixed cases. This indicates the potential of adaptive refinement for such class of problems. Acknowledgement This work is financially supported by the Research Council of Norway and industrial partners through the research project SAMCoT, Sustainable Arctic Marine and Coastal Technology. The authors gratefully acknowledge the support. References [1] Yared W Bekele, Trond Kvamsdal, Arne M Kvarving, and Steinar Nordal. Adaptive isogeometric finite element analysis of steady-state groundwater flow. International Journal for Numerical and Analytical Methods in Geomechanics, 2015. 14 [2] M Av Biot. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 26(2):182–185, 1955. [3] MA Biot. General solutions of the equations of elasticity and consolidation for a porous material. J. appl. Mech, 23(1):91–96, 1956. [4] Maurice A Biot. General theory of three-dimensional consolidation. Journal of applied physics, 12(2):155–164, 1941. [5] John R Booker and JC Small. An investigation of the stability of numerical solutions of Biot’s equations of consolidation. International Journal of Solids and Structures, 11(7):907–917, 1975. [6] J. Austin Cottrell, Thomas J. R. Hughes, and Yuri Bazilevs. Isogeometric Analysis : Toward Integration of CAD and FEA. Wiley, Chichester, West Sussex, U.K., Hoboken, NJ, 2009. [7] R de Boer and W Ehlers. A historical review of the formulation of porous media theories. Acta Mechanica, 74(1-4):1–8, 1988. [8] David Dureisseix, Pierre Ladevèze, and Bernard A Schrefler. A LATIN computational strategy for multiphysics problems: application to poroelasticity. International Journal for Numerical Methods in Engineering, 56(10):1489–1510, 2003. [9] Ashraf El-Hamalawi and MD Bolton. An a posteriori error estimator for plane-strain geotechnical analyses. Finite elements in analysis and design, 33(4):335–354, 1999. [10] Ashraf El-Hamalawi and MD Bolton. A-posteriori error estimation in axisymmetric geotechnical analyses. Computers and Geotechnics, 29(8):587–607, 2002. [11] Massimiliano Ferronato, Giuseppe Gambolati, and Pietro Teatini. Ill-conditioning of finite element poroelasticity equations. International Journal of Solids and Structures, 38(34):5995– 6014, 2001. [12] Giuseppe Gambolati, Giorgio Pini, and Massimiliano Ferronato. Numerical performance of projection methods in finite element consolidation models. International Journal for Numerical and Analytical Methods in Geomechanics, 25(14):1429–1447, 2001. [13] Jamshid Ghaboussi and Edward L Wilson. Flow of compressible fluid in porous elastic media. International Journal for Numerical Methods in Engineering, 5(3):419–442, 1973. [14] Joachim Berdal Haga, Harald Osnes, and Hans Petter Langtangen. On the causes of pressure oscillations in low-permeable and low-compressible porous media. International Journal for Numerical and Analytical Methods in Geomechanics, 36(12):1507–1522, 2012. [15] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(3941):4135–4195, 2005. [16] CoT Hwang, NR Morgenstern, and DW Murray. On solutions of plane strain consolidation problems by finite element methods. Canadian Geotechnical Journal, 8(1):109–118, 1971. [17] Faisal Irzal, Joris JC Remmers, Clemens V Verhoosel, and René de Borst. Isogeometric finite element analysis of poroelasticity. International Journal for Numerical and Analytical Methods in Geomechanics, 37(12):1891–1907, 2013. [18] Johannes Korsawe, Gerhard Starke, Wenqing Wang, and Olaf Kolditz. Finite element analysis of poro-elastic consolidation in porous media: Standard and mixed approaches. Computer Methods in Applied Mechanics and Engineering, 195(9):1096–1115, 2006. [19] Fredrik Larsson and Kenneth Runesson. A sequential-adaptive strategy in space-time with application to consolidation of porous media. Computer Methods in Applied Mechanics and Engineering, 288:146–171, 2015. [20] R.W. Lewis and B.A. Schrefler. The finite element method in the static and dynamic deformation and consolidation of porous media. Numerical methods in engineering. John Wiley, 1998. 15 [21] Márcio A Murad and Abimael FD Loula. Improved accuracy in finite element analysis of Biot’s consolidation problem. Computer Methods in Applied Mechanics and Engineering, 95(3):359–382, 1992. [22] Márcio A Murad and Abimael FD Loula. On stability and convergence of finite element approximations of Biot’s consolidation problem. International Journal for Numerical Methods in Engineering, 37(4):645–667, 1994. [23] Phillip Joseph Phillips and Mary F Wheeler. A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous-in-time case. Computational Geosciences, 11(2):131–144, 2007. [24] Phillip Joseph Phillips and Mary F Wheeler. A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discrete-in-time case. Computational Geosciences, 11(2):145–158, 2007. [25] Phillip Joseph Phillips and Mary F Wheeler. A coupling of mixed and discontinuous Galerkin finite-element methods for poroelasticity. Computational Geosciences, 12(4):417–435, 2008. [26] MB Reed. An investigation of numerical errors in the analysis of consolidation by finite elements. International journal for numerical and analytical methods in geomechanics, 8(3):243– 257, 1984. [27] Ranbir S Sandhu, Shyan Chyun Lee, and Hwie-Ing The. Special finite elements for analysis of soil consolidation. International journal for numerical and analytical methods in geomechanics, 9(2):125–147, 1985. [28] Ranbir S Sandhu, Honho Liu, and Kamar J Singh. Numerical performance of some finite element schemes for analysis of seepage in porous elastic media. International Journal for Numerical and Analytical Methods in Geomechanics, 1(2):177–194, 1977. [29] Ranbir S Sandhu and Edward L Wilson. Finite-element analysis of seepage in elastic media. Journal of the Engineering Mechanics Division, 95(3):641–652, 1969. [30] Scott W Sloan and Andrew J Abbo. Biot consolidation analysis with automatic time stepping and error control part 1: theory and implementation. International Journal for Numerical and Analytical Methods in Geomechanics, 23(6):467–492, 1999. [31] Scott W Sloan and Andrew J Abbo. Biot consolidation analysis with automatic time stepping and error control: Part 2: Applications. International Journal for Numerical and Analytical Methods in Geomechanics, 23(6):493–529, 1999. [32] Maria Tchonkova, John Peters, and Stein Sture. A new mixed finite element method for poro-elasticity. International journal for numerical and analytical methods in geomechanics, 32(6):579–606, 2008. [33] K von Terzaghi. Die berechnung der durchlassigkeitsziffer des tones aus dem verlauf der hydrodynamischen spannungserscheinungen. Sitzungsberichte der Akademie der Wissenschaften in Wien, Mathematisch-Naturwissenschaftliche Klasse, Abteilung IIa, 132:125–138, 1923. [34] Karl Terzaghi et al. Erdbaumechanik auf bodenphysikalischer grundlage. 1925. [35] PA Vermeer and A Verruijt. An accuracy condition for consolidation by finite elements. International Journal for numerical and analytical methods in geomechanics, 5(1):1–14, 1981. [36] Guofu Zhu, Jian-Hua Yin, and Shun-tim Luk. Numerical characteristics of a simple finite element formulation for consolidation analysis. Communications in Numerical Methods in Engineering, 20(10):767–775, 2004. 16