Two finite element approaches for Darcy and Darcy–Brinkman flow through deformable porous media—Mixed method vs. NURBS based (isogeometric) continuity

Abstract

Porous media problems with coupled finite deformations and flow are still quite challenging. If large deformations occur, a non homogeneous porosity field is almost inevitable, as the porosity can vary due to deformation or due to flow and pore pressure. From a numerical point of view this significantly changes the mathematical description of the problem. In this contribution, we will analyze the consequences of the occurrence of porosity gradients on the finite element formulations of saturated, two phase porous media for Darcy and Darcy–Brinkman flow. We will present two approaches which are capable of fulfilling the numerical requirements for convergence. The first will be an approach based on NURBS shape functions. The second is a mixed formulation including the porosity as primary solution field. In addition, we will propose a method to include no-penetration or no-slip constraints into the monolithic system using a dual Lagrange multiplier method. The potential of the proposed methods for solving porous flow structure interaction undergoing large deformations will be illustrated by numerical examples.

Introduction

Porous media models are of great interest due to their potential application in various fields of engineering, like geophysics, civil engineering, physical chemistry, acoustics and biomechanics. In geophysics such problems occur for example in melting of ice or snow layers (see  [1]), or flow within magma chambers (see  [2, chap. 11]). Additionally, flow and transport in fractured rock (see  [3]) and the analysis of seismic attenuation (see  [4]) are of interest in this field. In civil engineering, applications include flow in porous media, e.g. while analyzing seepage through soil, walls of water reservoirs or dams (see  [5], [6]). In chemical engineering, such approaches are used to model mass and heat transfer through membranes in fuel cells (see, e.g.  [7]) or in packed-bed reactor columns (see  [8]). Also, a large number of models for acoustic and poroelastic wave propagation have been developed (see, e.g.  [9], [10]). In biomechanics, porous media models have been utilized to model bones or soft tissues (see, e.g.  [11], [12]). On a smaller scale, porous media models have successfully been applied to tumor growth models (see, e.g.  [13], [14]) and cell mechanics (see, e.g.  [15], [16]).

Porous media consist of one structure phase (the skeleton) and one or several fluid phases, flowing through connected pores. A fully resolved model would lead to a surface coupled fluid–structure-interaction system. However, such a separate analysis of both structure and fluid phases would be computationally highly demanding due to the complex and in almost every application unknown geometry of the pore structure. Yet, in many cases such a direct numerical simulation is also not needed to answer the relevant questions. Therefore, special methods for describing porous media on a macroscopic level have been developed. Based on the mathematical homogenization theory, so called “Reference Volume Averaging” has been hugely applied in engineering science (see e.g.  [17], [18], [19], [20]). This method leads to a continuous description of the porous medium, where fluid and solid are modeled as overlapping continua and, hence, the actual interface is not resolved explicitly. Finally, a volume coupled fluid–structure-interaction problem is derived, which enables modeling of a porous medium on a macroscopic scale without presuming detailed knowledge of the pore geometry. Also, apart from Reference Volume Averaging, other approaches to obtain a change of scale of the considered problem exist. The so called asymptotic expansion (see e.g.  [21], [22], [23]) uses a multi-scale perturbation method to derive the macroscopic equations. In the thermodynamically constrained averaging theory microscale thermodynamical expressions are included into the averaging process  [24].

Depending on the assumptions made, different flow equations are obtained. The most common flow equation is the Darcy equation. The flow is driven by a pressure gradient and a reactive term accounts for the resistance opposing the flow when passing the solid phase. The Darcy–Brinkman equation, on the other hand, further includes viscous effects. In both cases, but most significantly for Darcy–Brinkman flow, special care needs to be taken if the porosity is to vary in space, due to deformation or pore pressure, for instance. This leads to stricter continuity requirements for the ansatz spaces in the finite element formulation. In this work, we will analyze the resulting equations and present two approaches which are suitable to solve such problems.

Motivated by the need for smooth ansatz functions, we will introduce a first formulation based on special shape functions having high continuity, namely Non-Uniform Rational Basis Splines (NURBS). NURBS have been the industrial standard for CAD-systems (Computer Aided Design). Recently they gained a lot of popularity also in the Computational Mechanics community through the work of Hughes et al.  [25]. There, NURBS were used as shape functions in the context of the finite element method and they even gave a name to this special version of the FEM, namely isogeometric analysis. Since then, isogeometric analysis has been applied to a variety of problems. One application that is of special interest to our work is given in  [26] where isogeometric analysis has been applied to poroelasticity but still was yet restricted to linear deformations. Hence, to our knowledge, this specific part of our contribution also presents the first application of an isogeometric approach to nonlinear, instationary, large deformation poroelastic problems. As NURBS basis functions of order $p$ can be up to $C^{p-1}$-continuous by construction, they perfectly meet the requirements posed by our specific porous medium problem.

The second formulation will be based on a mixed formulation. Thereby, we do not only choose the solid displacements, the fluid velocity and fluid pressure as primary fields (which is, of course, a mixed formulation itself), but also incorporate the porosity as primary variable, by weakly enforcing the constitutive relation between porosity, pressure and displacements. By doing so, the continuity requirements are reduced and occurring porosity gradients can be represented by first order shape functions for Darcy–Brinkman flow also. Mixed finite element approaches for porous media problems are not new. Most commonly they are used in extended porous medium problems, such as thermo-poroelasticity, transport of chemical substances or multiphase flow within porous media (see, e.g.  [13], [27]). In this work we analyze our mixed approach specifically regarding the continuity requirements of Darcy–Brinkman flow and compare it to the NURBS based approach.

Another difficulty is the application of complex boundary conditions. Here, especially moving, impermeable boundaries are of interest. In case of Darcy flow this implies a no-penetration constraint (no relative flow in normal direction of the boundary) and in case of Darcy–Brinkman flow a no-slip constraint (no relative flow at the boundary) is often applied. We will demonstrate how such constraints can be included into the monolithic solution scheme using a Lagrange-multiplier method. We will further condense the additional degrees of freedom of the Lagrange multiplier out of the linear system and therefore maintain or even reduce the overall system size.

The outline of this contribution is as follows: In the next section we will give a brief summary of the definition of NURBS. The third section will present the governing equation of incompressible flow through saturated porous media undergoing large deformations. Subsequently, in Section  4 the two finite element formulations will be derived, including the incorporation of impermeability constraints on the fluid within the monolithic framework. The proposed methodology will be tested by numerical examples in Section  5 followed by concluding remarks in Section  6.