Mixed Method for Isogeometric Analysis of Coupled Flow and Deformation in Poroelastic Media

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, for problems where plastic deformations are negligible. 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, [1]. 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 one-dimensional case by Terzaghi [2,3]. The formulations were later generalized for a three-dimensional case and extended by Biot [4,5,6]. The mathematical formulations have since been extensively studied by several researchers. 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 [7]. Hwang et al. [8] 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 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 [9] 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 [10] 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. [11]. 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 [12] 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 [13] 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. [14]. 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 problems in poroelasticity and the mathematical properties of the governing equations and the numerical solution were studied in a further great detail. Murad and Loula [15] 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 [16], 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 [17,18]. Automatic time increments were selected such that the temporal discretization error in the displacements is close to a specified tolerance. Ferronato et al. [19] 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. [20] studied the numerical performance of projection methods in finite element consolidation models. Dureisseix et al. [21] 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. [22]. Korsawe et al. [23] 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. [24], 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 [25,26]. They also studied a coupling of mixed and discontinuous Galerkin finite-element methods, [27]. Haga et al. [28] studied the causes of pressure oscillations in low-permeable and low-compressible media by presenting two, three and four field 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, as a numerical stabilization technique for ill-conditioning or solution oscillations, has been studied by very few researchers. Larsson and Runesson [29] presented a novel approach for space-time adaptive finite element analysis for the coupled consolidation problem in geomechanics. El-Hamalawi and Bolton [30] 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 [31]. Adaptive isogeometric finite element analysis using LR B-Splines [32] and recovery based error estimators [33], was presented for steady-state groundwater flow problems by Bekele et al. [34].

Isogeometric finite element analysis of poroelasticity was first presented by Irzal et al. [35]. 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. However, 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. Mixed finite element and NURBS-based isogeometric analysis for Darcy and Darcy–Brinkman flow through deformable porous media was presented in [36]. Other related studies include isogeometric analysis of fluid-saturated porous media including flow in the cracks [37], mixed isogeometric analysis of strongly coupled diffusion in porous materials [38] and isogeometric analysis of fracture propagation in saturated porous media due to a pressurised non-Newtonian fluid [39]. Regarding mixed methods for isogeometric finite elements we want to mention the use of so-called Mixed Integration Point (MIP) recently proposed in [40] to achieve the advantages of mixed formulation for avoiding locking in isogeometric analysis of thin shell problems [41,42] without explicitly introducing the mixed set of variables.

In this paper, we present mixed formulation for isogeometric analysis of poroelasticity. We demonstrate through numerical examples that mixed isogeometric analysis can overcome pressure oscillations with the right choice of continuity and meshing, without introducing additional field variables in the mathematical formulation. The results are discussed in comparison with related numerical studies 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,5] 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:

$$\nabla ·\underset{=\mathbf{\sigma }}{\underbrace{\left({\mathbf{\sigma }}^{\prime }+\alpha p^{\mathrm{f}}\mathit{I}\right)}}+\rho \mathit{b}=\mathbf{0}$$

(1)

where $\mathbf{\sigma }$ is the total stress, ${\mathbf{\sigma }}^{\prime }$ is the effective stress, $\alpha$ is Biot’s coefficient, $p^{\mathrm{f}}$ is the fluid pressure, $\mathit{I}$ is an identity matrix, $\rho$ is the overall density of the porous medium and $\mathit{b}$ represents body forces. The Biot coefficient $\alpha$ can be calculated from:

$$\alpha =1-\frac{K_{\mathrm{t}}}{K_{\mathrm{s}}}$$

(2)

where $K_{\mathrm{t}}$ and $K_{\mathrm{s}}$ are the bulk moduli of the porous medium and solid particles, respectively.

The constitutive equation for poroelasticity relates stress and strain linearly as:

$${\mathbf{\sigma }}^{\prime }=\mathit{D}:\mathit{\epsilon }$$

(3)

where $\mathit{D}$ is a fourth-order stiffness tensor. Small deformations are also assumed, so the strain $\mathbf{\epsilon }$ satisfies a linear first-order equation with respect to the displacement $\mathit{u}$,

$$\mathit{\epsilon }=\frac{1}{2}\left(\nabla \mathit{u}+{\nabla }^{⊺}\mathit{u}\right)$$

(4)

where $1/2(\nabla +{\nabla }^{⊺})$ is the symmetrized gradient operator i.e.,

$${\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

(5)

In the following, it will be convenient to lower tensors and higher differential operators to Voigt notation, which represents the symmetric $d\times d$ tensor ${\mathbf{\sigma }}^{\prime }$ as a $d(d+1)/2$-vector, which we will denote with a tilde:

$$\underset{{\mathbf{\sigma }}^{\prime }}{\underbrace{\left[\begin{array}{ccc}{\sigma }_{xx}^{\prime }& {\sigma }_{xy}^{\prime }& {\sigma }_{xz}^{\prime }\\ {\sigma }_{xy}^{\prime }& {\sigma }_{yy}^{\prime }& {\sigma }_{yz}^{\prime }\\ {\sigma }_{xz}^{\prime }& {\sigma }_{yz}^{\prime }& {\sigma }_{zz}^{\prime }\end{array}\right]}}⟺\underset{{\tilde{\mathbf{\sigma }}}^{\prime }}{\underbrace{{\left\{\begin{array}{cccccc}{\sigma }_{xx}^{\prime }& {\sigma }_{yy}^{\prime }& {\sigma }_{zz}^{\prime }& {\sigma }_{yz}^{\prime }& {\sigma }_{xz}^{\prime }& {\sigma }_{xy}^{\prime }\end{array}\right\}}^{⊺}}}$$

A similar conversion takes place for the strains, where the shear strains are replaced by the engineering shear strains:

$$\underset{\mathbf{\epsilon }}{\underbrace{\left[\begin{array}{ccc}{\epsilon }_{xx}& {\epsilon }_{xy}& {\epsilon }_{xz}\\ {\epsilon }_{xy}& {\epsilon }_{yy}& {\epsilon }_{yz}\\ {\epsilon }_{xz}& {\epsilon }_{yz}& {\epsilon }_{zz}\end{array}\right]}}⟺\underset{\tilde{\mathbf{\epsilon }}}{\underbrace{{\left\{\begin{array}{cccccc}{\epsilon }_{xx}& {\epsilon }_{yy}& {\epsilon }_{zz}& 2{\epsilon }_{yz}& 2{\epsilon }_{xz}& 2{\epsilon }_{xy}\end{array}\right\}}^{⊺}}}$$

Voigt notation allows us to express the equilibrium equation and the stress-strain equation using the same differential operator $\mathit{L}$,

$${\mathit{L}}^{⊺}=\left[\begin{array}{cccccc}\frac{\partial }{\partial x}& 0& 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial z}\\ 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial x}& \frac{\partial }{\partial z}& 0\\ 0& 0& \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]$$

(6)

Using $\mathit{L}$ yields the following equilibrium equation in terms of the two primary unknowns $\mathit{u}$ and $p^{\mathrm{f}}$,

$${\mathit{L}}^{⊺}\tilde{\mathit{D}}\mathit{L}\mathit{u}-\alpha \nabla p^{\mathrm{f}}+\rho \mathit{b}=\mathbf{0},$$

(7)

where $\tilde{\mathit{D}}$ is the Voigt notation equivalent of $\mathit{D}$, taking into account the aforementioned engineering shear strains. We will generally assume isotropic materials, where $\tilde{\mathit{D}}$ takes the block form (in terms of Young’s modulus E and Poisson’s ratio $\nu$)

$$\tilde{\mathit{D}}=\frac{E}{(1+\nu )(1-2\nu )}\left[\begin{array}{cc}{\tilde{\mathit{D}}}_{11}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathit{D}}}_{22}\end{array}\right]$$

(8)

where the two blocks are given as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tilde{\mathit{D}}}_{11}& =(1-2\nu )\mathit{I}+\nu \mathit{1}\hfill \\ \hfill {\tilde{\mathit{D}}}_{22}& =\frac{1-2\nu }{2}\mathit{I}\hfill \end{array}\end{array}$$

(9)

and $\mathit{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 $\zeta$ is given by

$$\zeta =\alpha \nabla ·\mathit{u}+c{p}^{\mathrm{f}}$$

(10)

where c is the storativity or specific storage coefficient at constant strain. It is given by

$$c=\frac{\alpha -n}{K_{\mathrm{s}}}+\frac{n}{K_{\mathrm{f}}}$$

(11)

where $K_{\mathrm{f}}$ is the bulk modulus of the fluid and n is the porosity of the material. The change in the fluid content $\zeta$ satisfies the equation

$$\frac{\partial \zeta }{\partial t}+\nabla ·\mathit{w}=0$$

(12)

where $\mathit{w}$ is the fluid flux, which is given by Darcy’s law as:

$$\mathit{w}=-\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}·\left(\nabla p^{\mathrm{f}}-{\rho }_{\mathrm{f}}\mathit{b}\right)$$

(13)

where ${\gamma }_{\mathrm{f}}$ is the unit weight of the fluid, ${\rho }_{\mathrm{f}}$ its density and $\mathit{k}$ is the hydraulic conductivity matrix.

The final equation of mass balance is then

$$\alpha \nabla ·\dot{\mathit{u}}+c\frac{\partial p^{\mathrm{f}}}{\partial t}+\nabla ·\left[-\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}·\left(\nabla p^{\mathrm{f}}-{\rho }_{\mathrm{f}}\mathit{b}\right)\right]=0.$$

(14)

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 boundary value problems. Let $({\mathrm{\Gamma }}_D^{\mathrm{u}},{\mathrm{\Gamma }}_D^{\mathrm{p}})$ and $({\mathrm{\Gamma }}_N^{\mathrm{u}},{\mathrm{\Gamma }}_N^{\mathrm{p}})$ be two partitions of the boundary $\partial \Omega$ of domain $\Omega$, for representing Dirichlet and Neumann boundary conditions, respectively.

The Dirichlet boundary conditions for the equilibrium (7) and mass balance (14) equations are

$$\left\{\begin{array}{cc}\mathit{u}=\overline{\mathit{u}}\hfill & \mathrm{on}{\mathrm{\Gamma }}_D^{\mathrm{u}},\hfill \\ p^{\mathrm{f}}={\overline{p}}^{\mathrm{f}}\hfill & \mathrm{on}{\mathrm{\Gamma }}_D^{\mathrm{p}},\hfill \end{array}\right.$$

(15)

where $\overline{\mathit{u}}$ and ${\overline{p}}_{\mathrm{f}}$ are the prescribed displacement and pressure, respectively.

The Neumann boundary conditions are

$$\left\{\begin{array}{cc}\mathbf{\sigma }·\mathit{n}=\overline{\mathit{t}}\hfill & \mathrm{on}{\mathrm{\Gamma }}_N^{\mathrm{u}},\hfill \\ \mathit{w}·\mathit{n}=\overline{q}\hfill & \mathrm{on}{\mathrm{\Gamma }}_N^{\mathrm{p}},\hfill \end{array}\right.$$

(16)

where $\mathit{n}$ is the outward pointing normal vector, $\overline{\mathit{t}}$ is the surface traction and $\overline{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 $\delta \mathit{u}$, which vanishes on ${\mathrm{\Gamma }}_D^{\mathrm{u}}$, and a scalar test function $\delta p$, which vanishes on ${\mathrm{\Gamma }}_D^{\mathrm{p}}$.

We start with the total stress formulation of the linear momentum balance equation, which from Equation (7) is given by

$$\nabla ·\mathbf{\sigma }+\rho \mathit{b}=0.$$

(17)

Multiplying by the test function $\delta \mathit{u}$ and integrating over the domain $\Omega$ gives

$${\int }_{\Omega }\delta {\mathit{u}}^{⊺}\nabla ·\mathbf{\sigma }\mathrm{d}\Omega +{\int }_{\Omega }\delta {\mathit{u}}^{⊺}\rho \mathit{b}\mathrm{d}\Omega =\mathbf{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,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\int }_{\Omega }\delta {\mathit{u}}^{⊺}\nabla ·\mathbf{\sigma }\mathrm{d}\Omega & =\sum _i{\int }_{\Omega }\delta u_i\nabla ·{\mathbf{\sigma }}_i\mathrm{d}\Omega \hfill \\ \hfill & =\sum _i{\int }_{\partial \Omega }\delta u_i{\mathbf{\sigma }}_i·\mathit{n}\mathrm{d}\mathrm{\Gamma }-\sum _i{\int }_{\Omega }\nabla \delta u_i·{\sigma }_i\mathrm{d}\Omega \hfill \\ \hfill & ={\int }_{{\mathrm{\Gamma }}_N^{\mathrm{u}}}\delta {\mathit{u}}^{⊺}\overline{\mathit{t}}\mathrm{d}\mathrm{\Gamma }-{\int }_{\Omega }\nabla \delta \mathit{u}:\mathbf{\sigma }\mathrm{d}\Omega .\hfill \end{array}\end{array}$$

(19)

Due to the symmetry of the stress tensor, the last term is expressible in Voigt notation,

$$\nabla \delta \mathit{u}:\mathbf{\sigma }={\left(\mathit{L}\delta \mathit{u}\right)}^{⊺}\mathbf{\sigma },$$

(20)

yielding the weak form of (7) as

$${\int }_{\Omega }{\left(\mathit{L}\delta \mathit{u}\right)}^{⊺}\tilde{\mathit{D}}\left(\mathit{L}\mathit{u}\right)\mathrm{d}\Omega -\alpha {\int }_{\Omega }{\left(\mathit{L}\delta \mathit{u}\right)}^{⊺}\tilde{\mathit{I}}{p}^{\mathrm{f}}\mathrm{d}\Omega ={\int }_{\Omega }\delta {\mathit{u}}^{⊺}\rho \mathit{b}\mathrm{d}\Omega +{\int }_{{\mathrm{\Gamma }}_N^{\mathrm{u}}}\delta {\mathit{u}}^{⊺}\overline{\mathit{t}}\mathrm{d}\mathrm{\Gamma }$$

(21)

where we have used $\tilde{\mathit{I}}$ as the Voigt notation identity operator, which for a general three-dimensional case is given by:

$$\tilde{\mathit{I}}={\left\{1,1,1,0,0,0\right\}}^{⊺}$$

(22)

For the mass balance equation, multiplying (14) by the scalar test function $\delta p$ and integrating over the domain $\Omega$, we get

$$\alpha {\int }_{\Omega }\delta p\nabla ·\dot{\mathit{u}}\mathrm{d}\Omega +c{\int }_{\Omega }\delta p\frac{\partial p^{\mathrm{f}}}{\partial t}\mathrm{d}\Omega +{\int }_{\Omega }\delta p\nabla ·\left[-\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}·\left(\nabla p^{\mathrm{f}}-{\rho }_{\mathrm{f}}\mathit{b}\right)\right]\mathrm{d}\Omega =0.$$

(23)

Again, by applying Green’s theorem to the last term, we obtain

$$\begin{array}{c}\hfill {\int }_{\Omega }\delta p\nabla ·\left[-\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}·\left(\nabla p^{\mathrm{f}}-{\rho }_{\mathrm{f}}\mathit{b}\right)\right]\mathrm{d}\Omega ={\int }_{{\mathrm{\Gamma }}_N^{\mathrm{p}}}\delta p\overline{q}\mathrm{d}\mathrm{\Gamma }\\ \hfill -{\int }_{\Omega }\nabla \delta p·\left[-\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}·\left(\nabla p^{\mathrm{f}}-{\rho }_{\mathrm{f}}\mathit{b}\right)\right]\mathrm{d}\Omega .\end{array}$$

(24)

Thus, the weak form of the mass balance equation, (14), is

$$\begin{array}{c}\hfill \alpha {\int }_{\Omega }\delta p\nabla ·\dot{\mathit{u}}\mathrm{d}\Omega +c{\int }_{\Omega }\delta p\frac{\partial p^{\mathrm{f}}}{\partial t}\mathrm{d}\Omega +{\int }_{\Omega }\nabla \delta p^{⊺}\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}\nabla p^{\mathrm{f}}\mathrm{d}\Omega \\ \hfill ={\int }_{\Omega }\nabla \delta p^{⊺}\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}{\rho }_{\mathrm{f}}\mathit{b}\mathrm{d}\Omega -{\int }_{{\mathrm{\Gamma }}_N^{\mathrm{p}}}\delta p\overline{q}\mathrm{d}\mathrm{\Gamma }.\end{array}$$

(25)

2.5. Galerkin Finite Element Formulation

With a suitable number N of basis functions defined, let ${\mathit{N}}_{\mathrm{p}}:\Omega \to {\mathbb{R}}^{1\times N}$ and ${\mathit{N}}_{\mathrm{u}}:\Omega \to {\mathbb{R}}^{d\times 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:

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill \mathit{u}& ={\mathit{N}}_{\mathrm{u}}{\mathit{N}}^{\mathrm{c}},\hfill & \hfill \delta \mathit{u}& ={\mathit{N}}_{\mathrm{u}}\delta {\mathit{u}}^{\mathrm{c}},\hfill \\ \hfill p^{\mathrm{f}}& ={\mathit{N}}_{\mathrm{p}}{\mathit{p}}^{\mathrm{c}},\hfill & \hfill \delta p& ={\mathit{N}}_{\mathrm{p}}\delta {\mathit{p}}^{\mathrm{c}}\hfill \end{array}\end{array}$$

(26)

where ${\mathit{u}}^{\mathrm{c}}$ and ${\mathit{p}}^{\mathrm{c}}$ 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 $\delta {\mathit{u}}^{\mathrm{c}}$ and $\mathbf{\delta }{p}^{\mathrm{c}}$, as Equations (7) and (14) are supposed to be valid for any choice of these)

$$\mathit{K}{\mathit{u}}^{\mathrm{c}}-\mathit{Q}{\mathit{p}}^{\mathrm{c}}={\mathit{f}}_{\mathrm{u}}$$

(27)

where the stiffness matrix $\mathit{K}$, the coupling matrix $\mathit{Q}$ and the vector of body forces and surface tractions ${\mathit{f}}_u$ are given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathit{K}& ={\int }_{\Omega }{\mathit{B}}^{⊺}\tilde{\mathit{D}}\mathit{B}\mathrm{d}\Omega ,\hfill \\ \hfill \mathit{Q}& ={\int }_{\Omega }{\mathit{B}}^{⊺}\alpha \tilde{\mathit{I}}{\mathit{N}}_{\mathrm{p}}\mathrm{d}\Omega \hfill \\ \hfill {\mathit{f}}_{\mathrm{u}}& ={\int }_{\Omega }{\mathit{N}}_{\mathrm{u}}^{⊺}\rho \mathit{b}\mathrm{d}\Omega ,+{\int }_{{\mathrm{\Gamma }}_N^{\mathrm{u}}}{\mathit{N}}_{\mathrm{u}}^{⊺}\overline{\mathit{t}}\mathrm{d}\mathrm{\Gamma }.\hfill \end{array}\end{array}$$

(28)

Here $\mathit{B}=\mathit{L}{\mathit{N}}_{\mathrm{u}}$ 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

$${\mathit{Q}}^{⊺}\frac{\partial {\mathit{u}}^{\mathrm{c}}}{\partial t}+\mathit{S}\frac{\partial {\mathit{p}}^{\mathrm{c}}}{\partial t}+\mathit{P}{\mathit{p}}^{\mathrm{c}}={\mathit{f}}_{\mathrm{p}}$$

(29)

where the storage matrix $\mathit{S}$, the permeability matrix $\mathit{P}$ and the vector of fluid body forces and fluxes ${\mathit{f}}_{\mathrm{p}}$ are given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathit{S}& ={\int }_{\Omega }{\mathit{N}}_{\mathrm{p}}^{⊺}c{\mathit{N}}_{\mathrm{p}}\mathrm{d}\Omega \hfill \\ \hfill \mathit{P}& ={\int }_{\Omega }\nabla {\mathit{N}}_{\mathrm{p}}^{⊺}\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}\nabla {\mathit{N}}_{\mathrm{p}}\mathrm{d}\Omega \hfill \\ \hfill {\mathit{f}}_{\mathrm{p}}& ={\int }_{\Omega }\nabla {\mathit{N}}_{\mathrm{p}}^{⊺}\frac{1}{{\gamma }_{\mathrm{f}}}\mathit{k}{\rho }_{\mathrm{f}}\mathit{b}\mathrm{d}\Omega -{\int }_{{\mathrm{\Gamma }}_N^{\mathrm{p}}}{\mathit{N}}_{\mathrm{p}}^{⊺}\overline{q}\mathrm{d}\mathrm{\Gamma }.\hfill \end{array}\end{array}$$

(30)

Combining Equations (27) and (29) results in the coupled system of equations for poroelasticity

$$\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ {\mathit{Q}}^{⊺}& \mathit{S}\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{u}}}^{\mathrm{c}}\\ {\dot{\mathit{p}}}^{\mathrm{c}}\end{array}\right\}+\left[\begin{array}{cc}\mathit{K}& -\mathit{Q}\\ \mathbf{0}& \mathit{P}\end{array}\right]\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{c}}\\ {\mathit{p}}^{\mathrm{c}}\end{array}\right\}=\left\{\begin{array}{c}{\mathit{f}}_{\mathrm{u}}\\ {\mathit{f}}_{\mathrm{p}}\end{array}\right\}.$$

(31)

A symmetric system of equations can be obtained by time-differentiating the first equation and multiplying one of the equations by $-1$, [43]:

$$\left[\begin{array}{cc}-\mathit{K}& \mathit{Q}\\ {\mathit{Q}}^{⊺}& \mathit{S}\end{array}\right]\left\{\begin{array}{c}{\dot{\mathit{u}}}^{\mathrm{c}}\\ {\dot{\mathit{p}}}^{\mathrm{c}}\end{array}\right\}+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathit{P}\end{array}\right]\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{c}}\\ {\mathit{p}}^{\mathrm{c}}\end{array}\right\}=\left\{\begin{array}{c}-{\dot{\mathit{f}}}_{\mathrm{u}}\\ {\mathit{f}}_{\mathrm{p}}\end{array}\right\}$$

(32)

In this formulation, it is important that time-dependent quantities involved in ${\mathit{f}}_u$, 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.

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 $\mathit{X}={\left\{{\mathit{u}}^c,{\mathit{p}}^c\right\}}^{⊺}$, we have the GTR approximation

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left.\frac{\partial \mathit{X}}{\partial t}\right|}_{n+\theta }& =\frac{{\mathit{X}}_{n+1}-{\mathit{X}}_n}{\Delta t}\hfill \\ \hfill {\mathit{X}}_{n+\theta }& =(1-\theta ){\mathit{X}}_n+\theta {\mathit{X}}_{n+1}\hfill \end{array}\end{array}$$

(33)

where $\theta$ is a time integration parameter which has limits $0\le \theta \le 1$ and n is a time step identifier. Adopting backward Euler time stepping ($\theta =1$) with time step $\Delta t$ and applying (32) and (33) we obtain the system of equations

$$\left[\begin{array}{cc}-\mathit{K}& \mathit{Q}\\ {\mathit{Q}}^{⊺}& \mathit{S}+\Delta t\mathit{P}\end{array}\right]{\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{c}}\\ {\mathit{p}}^{\mathrm{c}}\end{array}\right\}}_{n+1}=\left[\begin{array}{cc}-\mathit{K}& \mathit{Q}\\ {\mathit{Q}}^{⊺}& \mathit{S}\end{array}\right]{\left\{\begin{array}{c}{\mathit{u}}^{\mathrm{c}}\\ {\mathit{p}}^{\mathrm{c}}\end{array}\right\}}_n+\Delta t{\left\{\begin{array}{c}-{\dot{\mathit{f}}}_u\\ {\mathit{f}}_p\end{array}\right\}}_{n+1}$$

(34)

which is a linear system in this case, for poroelasticity, as the coefficient matrices are independent of the unknowns.

3. Isogeometric Analysis

3.1. Introduction

Since its first introduction by Hughes et al. [44], 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 $\Xi =\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{n+p+1}\}$, where ${\xi }_i\in \mathbb{R}$ is the ith knot, i is the knot index, $i=1,2,\dots ,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, ${\left\{N_{i,p}\right\}}_{i=1}^n$. The coefficients corresponding to these functions, ${\left\{{\mathit{X}}_i\right\}}_{i=1}^n$, are referred to as control points. The B-Spline basis functions are recursively defined starting with piecewise constants ($p=0$):

$$N_{i,0}\left(\xi \right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{\xi }_i\le \xi <{\xi }_{i+1}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(35)

For higher-order polynomial degrees ($p\ge 1$), the basis functions are defined by the Cox-de Boor recursion formula:

$$N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}\left(\xi \right)+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}\left(\xi \right)$$

(36)

B-Spline geometries, curves, surfaces and solids, are constructed from a linear combination of B-Spline basis functions. Given n basis functions $N_{i,p}$ and corresponding control points ${\mathit{P}}_i\in {\mathbb{R}}^d,i=1,2,\dots ,n$, a piecewise polynomial B-Spline curve is given by:

$$\mathit{C}\left(\xi \right)=\sum _{i=1}^n{N}_{i,p}\left(\xi \right){\mathit{P}}_i$$

(37)

Similarly, for a given control net ${\mathit{P}}_{i,j},i=1,2,\dots ,n,j=1,2,\dots ,m,$ polynomial orders p and q, and knot vectors $\Xi =\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{n+p+1}\}$, and $\mathcal{H}=\{{\eta }_1,{\eta }_2,\dots ,{\eta }_{m+q+1}\}$, a tensor product B-Spline surface is defined by:

$$\mathit{S}(\xi ,\eta )=\sum _{i=1}^n\sum _{j=1}^m{N}_{i,p}\left(\xi \right)M_{j,q}\left(\eta \right){\mathit{P}}_{i,j}$$

(38)

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 ${\mathbb{R}}^d$ is obtained by projective transformation of a B-Spline entity in ${\mathbb{R}}^{d+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. [45].

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 non-repeated 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 $\Xi =\left\{0,0,0,0,0,1,2,2,3,3,3,4,4,4,4,5,5,5,5,5\right\}$, 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 non-commutativity 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.

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, $p_{\mathrm{p}}$, is defined by the geometry construction and the polynomial order for the displacement, $p_{\mathrm{u}}$, is raised by one. Both $p_{\mathrm{p}}$ and $p_{\mathrm{u}}$ 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 $\Xi =\left\{0,0,1,1\right\}$ and $\mathcal{H}=\left\{0,0,1,1\right\}$ implies $p_{\mathrm{p}}=1$ and $p_{\mathrm{u}}=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 to 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 $p_0$ at the top boundary. Note that we use x and y to denote the horizontal and vertical directions, respectively, while it is also common to use z for the vertical direction in one-dimensional poroelastic problems as in [28]. 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:

$$\begin{array}{cc}\hfill \frac{p^{\mathrm{f}}(t,y)}{p_0}& =\frac{4}{\pi }\sum _{i=1}^{\infty }\frac{{(-1)}^{i-1}}{2i-1}exp\left[-{(2i-1)}^2\frac{{\pi }^2{t}_{\mathrm{s}}}{4}\right]cos\left[(2i-1)\frac{\pi y}{2h}\right]\hfill \end{array}$$

where the dimensionless time $t_{\mathrm{s}}$ is given as a function of the consolidation coefficient $c_{\mathrm{v}}$ and drainage path h (total height for one-way drainage) by:

$$t_{\mathrm{s}}=\frac{c_{\mathrm{v}}}{h^2}t.$$

(39)

The consolidation coefficient $c_{\mathrm{v}}$ is given by:

$$c_{\mathrm{v}}=\frac{(1-\nu )E\kappa }{(1+\nu )(1-2\nu )}$$

(40)

The material parameters used for this problem are given in

Table 1. Terzaghi’s problem: Load and material parameters.

Parameter	Value	Unit
External load, $p_0$	$1.0\times {10}^6$	$\mathrm{Pa}$
Hydraulic conductivity, k	$1.962\times {10}^{-14}$	${\mathrm{m}}^2$
Biot’s coefficient, $\alpha$	1.0	−
Young’s modulus, E	$6.0\times {10}^6$	$\mathrm{Pa}$
Poisson’s ratio, $\nu$	0.4	−
Storativity, c	0	${\mathrm{Pa}}^{-1}$
Body forces, $\mathit{b}$	$\mathbf{0}$	$\mathrm{N}$

, as used in [35]. 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 $p_{\mathrm{p}}=1,2,3$. The corresponding values for the displacement in a mixed formulation are $p_{\mathrm{u}}=2,3,4$. The number of elements used in the simulation is $N_{\mathrm{e}}=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 [12].

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.

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\times h=1\times 1$ with the same boundary conditions as in the previous case. For simplicity we choose the following material parameters: $\alpha =1$, $c=0$, $E=2/3$, $\nu =0.25$ and $\kappa =1$. The external load applied is $p_0=1$ and we assume no body forces i.e., $\mathit{b}=\mathbf{0}$.

This case was run with an increasing number of degrees of freedom using polynomial degrees $p_{\mathrm{p}}=1,2,3$ for the pressure and correspondingly $p_{\mathrm{u}}=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 $L^2$ error of the pressure field. The relative error based on the computed pressure values, ${\rho }_h$, is calculated from

$${\rho }_h=\frac{\parallel p_h^{\mathrm{f}}-p^{\mathrm{f}}{\parallel }_{L^2}}{\parallel p^{\mathrm{f}}{\parallel }_{L^2}}$$

(41)

where $p_h^{\mathrm{f}}$ and $p^{\mathrm{f}}$ 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.

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 [28]. 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 $p_0$ 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. The Haga problem: Load and material parameters.

Parameter	Value	Unit
External load, $p_0$	$1.0$	$\mathrm{Pa}$
Darcy coefficient, $k_1/{\gamma }_{\mathrm{f}}$	$1.0$	${\mathrm{m}}^2/\mathrm{Pa}\mathrm{s}$
Darcy coefficient, $k_2/{\gamma }_{\mathrm{f}}$	$1.0\times {10}^{-8}$	${\mathrm{m}}^2/\mathrm{Pa}\mathrm{s}$
Biot’s coefficient, $\alpha$	1.0	−
Young’s modulus, E	$0.67$	$\mathrm{Pa}$
Poisson’s ratio, $\nu$	0.25	−
Storativity, c	0	${\mathrm{Pa}}^{-1}$
Body forces, $\mathit{b}$	$\mathbf{0}$	$\mathrm{N}$

. 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., $p_{\mathrm{p}}=1,2,3,4$. The corresponding polynomial degrees for the displacement in a mixed formulation are $p_{\mathrm{u}}=2,3,4,5$. The continuities at the boundaries between the layers are also varied. We consider $C^0$ and $C^{p_{\mathrm{p}}-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. For short and early time steps coupled with a very low permeability, the solution to Equation (14) requires divergence-free displacements close to the dissipation boundary and in the layer with the low permeability. Requiring to satisfy this with equal order elements for displacement and pressure locks out most of the displacement degrees of freedom leading to oscillations in the pressure solution. Due to its high permeability, the fluid in the top layer dissipates very quickly for the time step size considered here i.e., $\Delta t=1s$. The pressure oscillations start 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, $p_{\mathrm{p}}=4$. The results with $C^0$ continuities at the material interfaces improve slightly better than with $C^{p_{\mathrm{p}}-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.

To compare our IGA analysis results with standard finite elements, we refer to the extensive finite element study on pressure oscillations presented in [28] where the low permeability layer problem is one of the numerical examples investigated. Haga et al. presented a systematic case study where they considered three different formulations for the mathematical model of poroelasticity and several combinations of finite elements. The mathematical formulations considered are (1) two-field formulation (displacement and fluid pressure), similar to our paper, (2) three-field formulation (displacement, velocity and fluid pressure) and (3) four-field formulation (displacement, soild pressure, velocity and fluid pressure). These different formulations were numerically investigated using several combinations of standard triangular and quadrilateral finite elements with different combinations of polynomial orders or continuities for the field variables. Lagrangian, Raviart–Thomas, Crouzeix–Raviart (triangular) and Rannacher–Turek (quadrilateral) elements were used for the numerical studies. For the low permeability problem, and as expected, equal order triangular and quadrilateral Lagrange finite elements (both first- and second-order) were found to fail in removing pressure oscillations. Mixed triangular and quadrilateral Lagrange elements (second-order for displacement and first-order for fluid pressure) were found to improve the pressure solution, but local pressure spikes were still exhibited. Satisfactory solutions were obtained only when using three- and four-field formulations with different element types for the field variables for example when representing the velocity field with linear Raviart–Thomas elements where the displacement and fluid pressure use second- and lowest-order triangular and quadrilateral Lagrangian elements. Our IGA analysis results produced satisfactory results based on the two-field formulation, without the need of introducing additional field variables and combinations of finite element types for the field variables.

5. Conclusions

Mixed isogeometric analysis of poroelasticity is presented where different orders 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 they tend to decrease with increasing polynomial orders for the pore pressure. This is illustrated by the convergence of the relative $L^2$ 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.

Our isogeometric mixed finite elements with maximum continuity show superior results compared to earlier studies of these problems using $C^0$-continuous mixed Lagrange finite elements. In particular, for the problem of consolidation of a very low permeability layer, the improved performance of the isogeometric finite element formulation (due to the higher regularity) is significant.