Abstract
The free finite element approach is applied to the governing equations describing the consolidation process in saturated poroelastic medium with intrinsically incompressible solid and fluid phases. Under this approach, where Voigt notation is avoided, the finite element equilibrium equations and the linearization of the coupled governing equations are fully derived using tensor algebra. In order to assess the free approach for the consolidation equations, direct comparison with analytical solution of the response of a homogeneous and isotropic water-saturated poroelastic finite column under harmonic load is presented. The results illustrate the capability of this finite element approach of reproducing accurately the response of quasistatic phenomena in a saturated porous medium.
1. Introduction
The range of interesting problems involving porous media has increased dramatically in recent years. These problems might arise from different environmental fields like agriculture, petroleum engineering, or natural hazards, to mention a few [1–3]. A key aspect in order to develop accurate models, capable of reproducing the principal features of this wide class of phenomena, is the adequate description of multiphase porous media. Such multiphase materials can be properly described within the thermodynamically consistent theory of porous media (TPM), where all constituents are assumed to occupy simultaneously the space [4].
In order to obtain accurate, robust, and efficient solutions of the coupled partial differential equations that come up from TPM, numerical techniques are usually required. In this context, the finite element method is one that has achieved significant results [5–11].
Due to the widespread application of the finite element technique in engineering modeling, it is of great relevance how this numerical procedure is taught in graduate and undergraduate courses. In this context, the free finite element approach recently proposed in [12] is of paramount significance. The main ideas behind the free finite element approach are the following ones:(1)Use of Voigt algebra is avoided within free approach. Under Voigt algebra, second-order symmetric tensors in 3 dimensions are represented as a -dimensional vectors and a fourth-order symmetric tensor is rewritten as a matrix. Instead of this vectorization approach, which might be error-prone [12], the so-called free approach is tensor based. Only vectors and second-order tensors and their natural operations are considered, at least in solid mechanics. Fourth-order tensors are employed through development of the finite element equations; however, the Cartesian components of such type of tensors will never appear within the implementation.(2)From the previous point, it can be inferred that the strain operator , which relates the strain vector (increment) to the nodal displacements (increment), is not required. In the same way, the constitutive matrix, appearing in the tangent stiffness matrix in the standard finite element implementation, is also not required.
The main contribution of the free finite element approach is that it allows a direct translation of the continuum mechanics operations into the finite element implementation, simplifying the standard approach.
To the best of the authors knowledge, free finite element approach has not been considered before within the context of the TPM, even though some tensorial based finite element approaches have been developed.
The free finite element approach is related to [13]. See also [14, 15]. In [13], the authors proposed an implementation where the use of and matrices is avoided. However, this work is restricted to isotropic and linear elastic nonporous material. Also, the fourth-order tensor of classical elasticity appears in the expression of the tangent stiffness matrix, making use of indicial notation to operate with it. However, the free approach can be employed irrespective of the material type, which applies to the computation of internal forces, and no fourth-order tensors ever appear within the implementation.
Regarding the numerical resolution of saturated porous media equations, it is worth mentioning the tensorial form of the finite element discretization of the Generalized Biot field equations [5] presented by Zienkiewicz and coworkers [7]. The tensorial approach considered by these authors follows [13]. However, as stated by the authors, Voigt notation is preferred for the finite element computation. In other words, Zienkiewicz and coworkers presented the tensorial form for the sake of completeness of the referred text, but it seems that this approach was avoided in the finite element implementation.
It is also worth mentioning [16]. In this work, the author proposed a tensorial form of the finite element method for the formulation within TPM framework. Again the tensorial approach presented in this work presents similarities with [13]. However, in [16], the tensorial form is extended to deal with unsaturated porous media. Also the material model considered for the solid skeleton was extended to deal with a Drucker-Prager elastic-perfectly plastic model. On the other hand, a fourth-order tensor (consistent operator) appears in the expression of the tangent stiffness matrix in [16], making use of indicial notation to operate with it. However, in the free approach, the Cartesian components of a fourth-order tensors will never appear explicitly, irrespective of the material type considered.
An important point is whether the free approach is more efficient than the standard one. In [12], the authors stated that this is true in certain problems. However, this claim is not easy to sustain as there are specific aspects of the implementation, sparsity of the matrices and parallelization of process are only two of them, that might complicate the task of efficiency assessment. In the present work, this type of assessment is not accomplished.
The paper is structured as follows. For the sake of completeness, in Section 2, the consolidation equations of a fluid saturated porous medium within the theory of porous media (TPM) [4, 17] will be derived. Assumptions considered to obtain the full set of equations will be stated first. Then, the balance equations of mass and momentum of each phase will be presented under the scope of the assumptions made. In Section 3, the variational formulation of the consolidation equations is presented. In Section 4, the tensor based free finite element approach is developed in detail for the present problem. Evaluation of the accuracy of numerical procedures is of paramount importance in practical engineering. Therefore, in Section 5, the free numerical technique is applied to investigate the response of a one-dimensional column of finite length, composed by a fluid saturated poroelastic material and subjected to a time dependent loading. In order to assess the accuracy of the free approach, a direct comparison with an analytical solution is presented.
2. Consolidation Equations for Saturated Soils
The assumptions considered to obtain the full set of governing equations are the following ones:(1)Porous medium is fully saturated.(2)Solid and fluid phases are intrinsically incompressible.(3)Neither mass nor heat exchanges between the solid and the fluid phase are considered.(4)Viscous shear stress within the fluid is neglected.(5)Permeability is considered isotropic.(6)Porous solid is considered isotropic and linear elastic under small strain theory.(7)All dynamic terms, including convective terms, are neglected.(8)Only loading by external forces is considered. Body forces are excluded.
Within the framework of TPM, a fluid saturated porous medium is described by two superimposed but immiscible species or phases with particles . The superscript denotes pore fluid for and solid skeleton for .
The kinematics in TPM is based on the assumption [4] that the spatial position in the current configuration at time is simultaneously occupied by particles of both constituents . Each particle proceeds from its own reference position in the reference configuration at time ; therefore, each phase has its own motion and at the common spatial position the phases interact (Figure 1).
The following relation holds at time :
An important concept within the modern TPM, responsible for the notion of superimpose phases over the control space, is the volume fraction of the phase , defined bywhere is the element volume of the phase in the current configuration within a representative element volume of the overall porous media. The volume fractions in (2) satisfy the following relation:
Partial densities of each phase are defined as the averaged density with the volume fraction ; that is,where is the real or intrinsic density of the phase . As the phases are assumed incompressible in this work, the intrinsic density of each phase is kept constant during the deformation process; that is,
Excluding mass and heat exchanges between the solid and the fluid phases, the balance equations for each phase are given as follows [17]: Balance of mass of phase : Balance of momentum of phase :
In expression (6), is the velocity of the particles defined bywhile is the material derivative with respect to the phase defined by
In expression (7), is the partial Cauchy stress second-order tensor of the phase , is the gravity acceleration vector, and is the acceleration vector of the particles which is defined in accordance with (9). The vector represents the interaction force between phase and the other phase and satisfies the summing rule:
The symbol in (9) means gradient operator with respect to the spatial position , while the operator in (7) is the divergence associated with .
Combining (4), (5), and (6), it follows
Now, by (3), (9), and (11), the volume balance of the mixture is obtained:
It is worth recalling that the motion of fluid within a porous medium is usually expressed in terms of an Eulerian coordinate system with respect to the current configuration of the solid skeleton [6]. Under this point of view, the relative velocity of the fluid with respect to the velocity of the solid is introduced:
Replacing (13) into (12), keeping (3) in mind, yields
As a consequence of the incompressibility constraint (5) for , the partial stress tensor for the solid can be additively decomposed into two terms:where is the Cauchy effective stress second-order tensor, is the second-order identity tensor, and is the pore fluid pressure. Due to assumption , the partial stress tensor associated with this phase can be taken as
By adding the partial stresses (15) and (16), the total Cauchy stress second-order tensor is recovered via Terzaghi’s principle:
Replacing (15), (16) in (7) with and , we arrive to the expressions
The volume fraction in (18), and by (3) also in (19), can be computed by solving analytically the ordinary differential equation (11) with . Assuming a small strain regime, this yields
However, since in a small strain setting, may be approximated by , which is the solid volume fraction in the reference configuration.
In order to fulfill the set of (18) and (19), constitutive relations must be formulated for the interaction force vector for and and for the Cauchy effective stress tensor .
Concerning the interaction vector , the relative velocity defined in (13) allows its definition as explained in [4, 11, 18]. In the case of isotropic permeability, the vector can be expressed aswhere is the gravity acceleration and is the Darcy permeability coefficient of the porous medium. Vector is obtained from (21) through (10). Replacing these expressions for the interaction force vectors in (18) and (19) yields
Regarding the constitutive relation for the Cauchy effective stress tensor , it is convenient to introduce the solid skeleton displacement vector defined by
This work is restricted to an isotropic linear elastic porous solid in small strain regime. Therefore, the effective stress tensor is determined by the Hookean elasticity law:whereis the fourth-order elastic tensor, while is the small strain tensor. In (25), is the second-order identity tensor and is the symmetric fourth-order identity tensor, while and are the macroscopic Lamé constants. The double dot in (24) means standard double contraction operation.
As the present work is developed within the framework of small strain theory, the superposition principal holds; that is, the loading by body forces and the external forces can be treated separately. Taking into account only loading by external forces, the equations of motion (22) can be written as
As a quasistatic slow motion phenomena is considered in the present work, it is reasonable to neglect all dynamic terms in (26) [6, 7, 11, 15]. These yield
Adding (27) and (28), solving for in (28), and taking into account relation (17) yieldwhere is the total Cauchy stress second-order tensor and (30) is no other than Darcy’s law. The set of coupled partial differential equations (29), (30), and (14) allows the determination of and in a consolidation process where loading is exclusively due to external forces.
3. Variational Formulation of the Consolidation Equations
The full set of coupled partial differential governing equations is given by the following: Mixture momentum balance: Mixture volume balance:together with the following equations to close the system: Darcy’s law: Constant fluid volume fraction: Terzaghi’s principle: Solid constitutive relation: Solid strain tensor: Solid displacement vector: Solid velocity vector:where the overdot in (39) indicates partial derivative with respect time . The primary variables for this system are and . Appropriate initial and boundary conditions should be specified for this set of equations in order to have a well pose initial boundary value problem.
Concerning the boundary conditions a fully saturated mixture is considered occupying the domain , where its boundary is decomposed in the following two different ways [7]: where the overline denotes closure while and are nonintersecting parts of the total boundary on which solid displacement and total tractions are prescribed, respectively. On the other hand, and are nonintersecting parts of the total boundary on which pore fluid pressure and fluid flux are prescribed, respectively. The following boundary conditions are prescribed in the above-mentioned parts of the boundary: Solid displacement specified: Total traction specified: Pore fluid pressure specified: Fluid flux specified:where and are given space and time vector functions, while and are given space and time scalar functions. The vector is the outward unit normal vector to the boundary and is the Darcy velocity.
The trial solution spaces and and the variation spaces and that are needed to define the weak or variational formulation of the problem can be described as follows:where and are, respectively, the scalar and vectorial Sobolev spaces of degree 1 [14].
Making use of Gauss divergence theorem and the definitions (47) and (48) of the variational spaces, the variational formulation of (31) and (32) reads as follows:
In standard finite element procedure, the in (49) is usually changed by , which is the symmetric part of the tensor . This can be done because Cauchy total stress tensor is symmetric. This small modification of the equations points towards the construction of the strain operator . As no such operator exists in the free approach, this change is avoided in the present work.
Integration in time of (49) and (50) is performed to obtain the discrete evolutions of and . With this objective, it is assumed that and are known at time and determination of and at time is wanted (subscript “” is omitted for brevity).
Equation (49) is an elliptic equation for which there is no proper time integration but evaluation of the variables at time . However, (50) is first order in time for which a generalized trapezoidal time integration method is considered [14]. In this context, the following approximation for at time is considered:where and , while . If , an unconditionally stable method is obtained. Also in (51), and are known values of and at time , respectively. Substituting (51) into (50), the variational problem reads as follows:where , , with , are defined in (45) and (46), respectively. Spaces and are defined in (47) and (48), respectively, while and are defined by
For the solution of the variational problem (52), at a fixed time , the standard Bubnov-Galerkin approach is applied. In the following section, the free finite element approach [12] will be described in detail.
4. B Free Finite Element Approach
To obtain the finite element discretization of the problem (52), the domain is subdivided into elements , connecting nodes. Mixed finite element formulation with equal order of interpolation for displacements and fluid pressure is employed. In this way, each node has degrees of freedom, where . The independent variables are the vector of nodal displacements and nodal fluid pressure , where .
Isoparametric shape functions are defined in standard fashion where is the reference element associated with the local orthogonal coordinate system : that is, for every , there is one and only one such that . These shape functions span the finite dimensional subspaces and of and (subscripts and are omitted for brevity), respectively; that is, where and are nodal values of the trial functions, already evaluated at a fixed time , while and are nodal values of the weighting functions.
Replacing the approximation expression (54)1 into (49) and adopting Einstein convention for repeated indices yieldtaking into account the following relations from tensor calculus:where stands for the standard dyadic product operation, and, by bilinearity of the scalar product, (56) reads asby noting that the nodal values of the weighting functions are arbitrary, it follows
Making use of Terzaghi’s principle (35) yields the first set of nodal equilibrium equations:
The following remarks might be stated regarding (61). Firstly, as isoparametric elements are considered, is computed by the following expression:where , while is the gradient operator with respect to the local orthogonal coordinate system . Secondly, the calculation of the internal force vector in (61) can be implemented as the multiplication of a matrix, representing the matrix of components of , with a vector ; that is, the strain operator is unnecessary for such a multiplication.
Now, replacing the approximation expression (54)2 into (50) yieldsagain, recalling from tensor calculus the following relation:The second set of nodal equilibrium equations is obtained as
Again, some aspects need to be clarified. Firstly, keeping in mind expressions (55)1 and (39), the term in (65) can be computed aswhere the nodal values are obtained from via (51). Secondly, Darcy’s law (33) and relation (64) yield
Simultaneous solution of algebraic system of (61) and (65) for the unknown nodal displacements and nodal fluid pressure will usually demand an iterative process like Newton Raphson.
Newton-Raphson iteration process for the system of (61) and (65) provides a vector sequence with , of displacement and pressure values at node by the following iterative scheme:
The matrix at the left-hand side in (68) is the tangent matrix of the coupled problem. The subscript in (68) and (69) means evaluation at .
Given a small value TOL (usual ~10−10), the Newton-Raphson process finished whenever the value of an adequate norm of the vector is less than TOL. As the present work is restricted to an isotropic linear elastic setting, convergence is achieved in the first iteration.
free finite element approach will be completed once the tangent matrix in (68) is described at element level. To this end, each of the four blocks , , , and in (68) should be specified.
Starting by , consider (61), with integration performed at element level; then,
Instead of computing the partial derivative of the integrand in (70) directly, by means of indicial notation [7, 13, 15, 16], the tensorial based procedure suggested in [12] will be applied.
Making use of the minor symmetries of the fourth-order elastic tensor , replacing approximation (55)1, and applying relation (57), (36) yields
The internal virtual work density in expression (61) can thus be written aswhere the right-hand side in (72) can be expressed as
In expression (73), the second-order tensor is defined by the relation [12]:for any vectors , , , , and fourth-order tensor . Replacing now the right-hand side of (73) in (70) yields
Taking into account (74) and keeping in mind the definition of the fourth-order elastic tensor (25), the following expression is obtained:
Some remarks follow. Firstly, due to expression (76), no fourt-order tensors appear within the tangent stiffness block making its computation simpler. The consideration of material symmetries and the use of intrinsic tensor algebra (relation (74)) prevent from using indicial notation with fourth-order tensors. Secondly, if an elastoplastic material is considered in a small strain setting; that is, , the tangent stiffness block will have the same structure as (75). However, in this case the expression for must be modified in accordance with the constitutive relation [12]. Thirdly, for a given material defined by , the implementation of the tangent block (75) requires only writing a subroutine that computes for any pair of vectors and .
The block of the tangent matrix (68) will be described next. Again, by (61), the following expression, at element level, is obtained:
Taking into account the approximation (55)2, expression (77) becomes
For the following block , let us consider (65). Due to the Darcy’s law (33), the only term appearing in (65) depending on is ; therefore, By (66) and (51), (79) yields
Finally, the block is described. By (65), this last block can be rewritten asBy (67), (81) becomes
After considering numerical integration, such as Gauss quadrature, and performing the standard assembling procedure in (75), (78), (80), and (82), and also for (61) and (65), new values of nodal displacements and nodal fluid pressure might be obtained at time .
It is worth mentioning some remarks regarding equal order of interpolation for displacements and fluid pressure employed in the present work. Firstly, it is well known that, in the incompressible and nearly incompressible regimes, fluid phase incompressible ((5) with ) with very low permeability equal-order interpolation for displacements and pressure cause spurious oscillation in the pressure field, unless some form of stabilization is utilized [7, 19, 20]. In the present work, we do not employ any stabilization technique for the sake of clarity of the free approach. In order to avoid the above-mentioned spurious oscillation in the simulations performed within the present work, very low values of permeability are avoided.
Secondly, as explained in [21], due to the strong coupling presented in the governing equations (31)–(39), linear interpolation of fluid pressure might yield wrong pressure values when high pressure gradient is observed, thus falsifying the overall displacement solution. In Section 5, where a comparative study between a numerical solution, obtained by the free approach, and an analytical solution is presented, eight node mixed quadrilateral finite elements will be considered to interpolate both displacements and fluid pressure fields.
5. Numerical Application
In the previous sections, the free finite element approach has been described within the context of the theory of porous media without dynamic terms. In the present section, this numerical technique is applied to investigate the response of a one-dimensional column of length , composed by a fluid saturated poroelastic material and subjected to a time dependent loading . The load is applied at the upper boundary of the column which is perfectly drained. This problem has been investigated previously by other researches [22–24].
The fluid pressure nearby the perfectly drained boundary, where the load is applied, might undergo a steep gradient. Therefore, to obtain an accurate fluid pressure distribution along with skeleton displacement at positions close to this boundary might become a challenging task for any numerical approach. In order to assess the accuracy of the free approach, direct comparison with an analytical solution is performed.
The motion of the solid skeleton within the column is constrained to take place in the vertical direction, defined by the axis, which is considered normal to the top boundary and pointing downward (Figure 2).
The bottom of the column is considered impermeable and fixed. Zero vertical displacement is assumed along the column at initial time.
For a one-dimensional problem, the governing equations (31)–(39) can be expressed as
In (83), the vertical value of has been substituted by , the fluid pressure by , the Darcy’s permeability by , and the intrinsic density of the fluid phase by . Subscript , like in or in , means partial differentiation with respect to the spatial coordinate , while subscript , for example in , means partial differentiation with respect to time. Appropriate initial and boundary conditions should be specified along with (83). Due to the above description of the problem, the following boundary and initial conditions are considered: Boundary conditions for are Boundary conditions for are Initial condition for is
Applying the eigenfunction expansion method, the following analytical solution for the initial boundary value problem defined by the system (83) subjected to conditions (84), (85), and (86) might be obtained:whereand in (88).
Finite number of terms in the series appearing in (87) should be considered in order to obtain a semianalytical solution to be compared with. This number of terms is established in such a way that further inclusion of more terms produces no changes in the solution.
The comparison has been performed with an acceleration of gravity m/s2, length of the column m, a harmonic time dependent loading kN/m2, and the following values for the material parameters:
Regarding the permeability, three different values have been considered from a high permeability m/s to a moderate value m/s. These sets of geometrical, loading, and material parameters yield a case of study that can be well classified as quasistatic slow motion phenomena [22].
In order to reproduce numerically the above case of study, a plane strain confined compression simulation is performed where the mixture domain is surrounded by impermeable, frictionless but rigid boundaries, except for the loaded top side which is perfectly drained. Geometry, boundary conditions, mesh, and the loading path are illustrated in Figure 3.
(a)
(b)
As it can be observed in Figure 3, eight node mixed quadrilateral finite elements are considered. This type of elements is employed to interpolate both displacements and fluid pressure fields. No pressure stabilization technique is considered, while full integration is employed for all the elements. For the generalized trapezoidal scheme, we have considered ; thus, the time integration method is unconditionally stable. A time step = 5·10−2 s is considered.
In Figures 4, 5, and 6, the vertical displacement at the top boundary and the fluid pressure distribution along the column can be observed for permeabilities m/s, m/s, and m/s.
(a)
(b)
(a)
(b)
(a)
(b)
It is clear that the free finite element approach is capable of reproducing accurately the response of a fluid saturated mixture under a harmonic loading. The comparison displayed in Figures 4, 5, and 6 shows good agreement both in vertical displacement and pore pressure distribution.
Keeping in mind that only two finite elements are considered per meter, the close profile between the numerical and semianalytical solutions is remarkable. Even for the case with m/s, where the steeper pressure gradient is observed, the numerical solution closely follows the semianalytical one.
6. Conclusions
In this work, the free finite element approach is applied to the governing equations describing the consolidation process in saturated poroelastic medium. A detailed description of the finite element discretization under this approach has been presented.
Only vectors and second-order tensors that are directly obtained from the continuum mechanics derivation appear within the finite element equations. In this way, the residual equations and their linearization have no extra ingredients, such as the so-called or matrices. Even though this aspect is crucial in the method, the main advantage of the free approach appears, while computing the tangent stiffness matrix. The consideration of material symmetries and the use of intrinsic tensor algebra prevent from using indicial notation with fourth-order tensors.
For a given material, the implementation of the tangent matrix, under this numerical approach, requires only writing easy subroutines based in standard tensor algebra operations.
It has been shown that this approach can be considered to reproduce accurately the response of quasistatic phenomena in a saturated poroelastic medium in small strains. The procedure presented in this work can be easily extended to deal with dynamic terms as well as with more advanced poroelastoplastic media.
Competing Interests
The authors declare no competing interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the support from the Spanish Ministry of Science and Innovation under Grant BIA2012-37020-C02-01 (GEOFLOW). In addition, M. M. Stickle would like to thank Jose Entrecanales Ibarra Foundation for the mobility grant funding his visiting scholar stay at Stanford University.