In this paper, we consider coupled hydro-mechanical processes in variably saturated porous media governed by the coupled nonlinear partial differential equations:
- Fluid mass conservation equation
(1)
- Darcy's law
(2)
- Equilibrium equation
(3)
- Generalized Hooke's law
(4)
- Van Genuchten (1980) water content - pressure head relationship
(5)
- Mualem (1976) conductivity-saturation relationship
(6)
with the saturation [-], the current water content [L3L−3], the saturated water content [L3L−3], [-] the effective saturation, the residual water content [L3L−3], the specific mass storativity related to head changes [L-1], the hydraulic head [L], the pressure head, the pressure [Pa], the fluid density [ML-3], the gravity acceleration [LT-2], the upward vertical coordinate [L], the specific moisture capacity [L−1], the Darcy velocity [LT-1], the sink term [T-1], the hydraulic conductivity [LT-1], the permeability tensor [L2], the fluid dynamic viscosity [ML-1T-1], the relative conductivity [-], [L-1] and [-] the van Genuchten parameters, , the Biot coefficient, the bulk density, the density of the solid [ML-3], the effective stress tensor, the Bishop's function typically set equal to , the displacement field [L], and the Lamé coefficients.
Several challenges emerge from the nonlinear character and the high coupling between the equations when solving the system (1)-(6). In this work, the equations are solved simultaneously using the method of lines (MOL) which avoids operator-splitting errors. The MOL is an efficient procedure for solving highly nonlinear systems of equations (Miller et al., 2006). With MOL, all the spatial derivatives are discretized while the time derivatives are kept in their continuous form.
For practical reasons, related to the computational burden of the system (1)-(6), the spatial discretization is often limited to low-order approximation methods, such as standard Galerkin Finite Elements (GFE). However, the standard GFE method can produce unstable and oscillatory pressure results, which is known as locking.
In this work, we use the low-order Crouzeix-Raviart (CR) finite elements for both the hydraulic head of the fluid phase and the displacement of the solid phase. The CR method uses P1 linear test functions with the degrees of freedom allocated to center of the edges, rather than to the vertices. Contrarily to the standard GFEs, the CR method is locally conservative. Further, using the connection between the Discontinuous Galerkin method and the CR method, Hansbo and Larson (2003) developed a locking-free CR formulation for elasticity. This formulation is employed here for poroelasticity in unsaturated porous media.
For the fluid flow, the nonlinear Richard's equation is also discretized in space using the CR method which is equivalent to the Lumped Raviart-Thomas Mixed Finite Element Method developed in Younes et al. (2006). The time discretization of the obtained system is performed using high-order integration methods and an efficient adaptive time stepping scheme with the DASPK time solver.
Numerical results in saturated and unsaturated conditions are presented to validate the new model and to show the effectiveness of the model to overcome nonphysical pressure oscillations.