B Free Finite Element Approach for Saturated Porous Media : Consolidation

The B 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 B 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.


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][2][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][6][7][8][9][10][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 B free finite element approach recently proposed in [12] is of paramount significance.The main ideas behind the B free finite element approach are the following ones: (1) Use of Voigt algebra is avoided within B free approach.
Under Voigt algebra, second-order symmetric tensors in 3 dimensions are represented as a 6-dimensional vectors and a fourth-order symmetric tensor is rewritten as a 6 × 6 matrix.Instead of this vectorization approach, which might be error-prone [12], the so-called B 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 B, which relates the strain vector (increment) to the nodal displacements (increment), is not required.In the same way, the constitutive D matrix, appearing in the tangent stiffness matrix in the standard finite element implementation, is also not required.
The main contribution of the B 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, B 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 B free finite element approach is related to [13].See also [14,15].In [13], the authors proposed an implementation where the use of B and D 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 B 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 B free approach, the Cartesian components of a fourthorder tensors will never appear explicitly, irrespective of the material type considered.
An important point is whether the B 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 B 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 B 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 B free approach, a direct comparison with an analytical solution is presented.

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.
(3) Neither mass nor heat exchanges between the solid and the fluid phase are considered.
(4) Viscous shear stress within the fluid is neglected.
(6) Porous solid is considered isotropic and linear elastic under small strain theory.
(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 x in the current configuration at time  is simultaneously occupied by particles   of both constituents ( = , ).Each particle   ( = , ) proceeds from its own reference position X  in the reference configuration at time  0 ; therefore, each phase has its own motion   and at the common spatial position x the phases interact (Figure 1).

The following relation holds at time 𝑡:
x =   (X  , ) . (1) An important concept within the modern TPM, responsible for the notion of superimpose phases over the control space, is the volume fraction   (x, ) of the phase , defined by where V  is the element volume of the phase  in the current configuration within a representative element volume V of the overall porous media.The volume fractions   (x, ) 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), k  is the velocity of the particles   defined by while   ◻/ is the material derivative with respect to the phase  defined by In expression (7),   is the partial Cauchy stress secondorder tensor of the phase , b is the gravity acceleration vector, and a  is the acceleration vector of the particles   which is defined in accordance with (9).The vector r  represents the interaction force between phase  and the other phase and satisfies the summing rule: The symbol grad(◻) in ( 9) means gradient operator with respect to the spatial position x, while the operator div(◻) in (7) is the divergence associated with grad(◻).
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 k  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, I is the second-order identity tensor, and   is the pore fluid pressure.Due to assumption (4), 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 div 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 div(u  ) ≪ 1 in a small strain setting,   may be approximated by   0 , 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 r  for  =  and  =  and for the Cauchy effective stress tensor   .
Concerning the interaction vector r  , the relative velocity defined in (13) allows its definition as explained in [4,11,18].In the case of isotropic permeability, the vector r  can be expressed as where  is the gravity acceleration and   is the Darcy permeability coefficient of the porous medium.Vector r  is obtained from ( 21) through (10).Replacing these expressions for the interaction force vectors in ( 18) and ( 19) Regarding the constitutive relation for the Cauchy effective stress tensor   , it is convenient to introduce the solid skeleton displacement vector u  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: where is the fourth-order elastic tensor, while  = (1/2)(grad(u  ) + grad T (u  )) is the small strain tensor.In (25), I is the secondorder identity tensor and I  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 div 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 div Adding ( 27) and (28), solving for k  in (28), and taking into account relation (17) where  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 u  and   in a consolidation process where loading is exclusively due to external forces.

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 u  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 Ω u  and Ω t 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 t are given space and time vector functions, while p and q are given space and time scalar functions.The vector n is the outward unit normal vector to the boundary Ω and   k  is the Darcy velocity.
The trial solution spaces S  () and S u () and the variation spaces V  and V u that are needed to define the weak or variational formulation of the problem can be described as follows: where  1 and H 1 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 grad() in ( 49) is usually changed by grad  (), which is the symmetric part of the tensor grad().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 B. As no such operator exists in the B 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 u  and   .With this objective, it is assumed that u   and    are known at time   and determination of u  and   at time  =   + Δ is wanted (subscript " + 1" 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 k  at time  =   + Δ is considered: where  1 = 1/Δ and  2 = (1 − )/, while  ∈ [0, 1].
If  ⩾ 1/2, an unconditionally stable method is obtained.Also in (51), u   and k   are known values of u  and k  at time   , respectively.Substituting (51) into (50), the variational problem reads as follows: where S u (), S  (), with  =   + Δ, are defined in (45) and (46), respectively.Spaces V u and V  are defined in (47) and (48), respectively, while L and M 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 B free finite element approach [12] will be described in detail.

B Free Finite Element Approach
To obtain the finite element discretization of the problem (52), the domain Ω is subdivided into   elements Ω  , connecting  nd nodes.Mixed finite element formulation with equal order of interpolation for displacements u  and fluid pressure   is employed.In this way, each node has  sd + 1 degrees of freedom, where  sd = spatial dimension.The independent variables are the vector of nodal displacements u   and nodal fluid pressure    , where  ∈ {1, . . .,  nd }.Isoparametric shape functions   : Ω  → R are defined in standard fashion where Ω  is the reference element associated with the local orthogonal coordinate system : that is, for every x ∈ Ω  , there is one and only one  ∈ Ω  such that x = x().These shape functions span the finite dimensional subspaces S ℎ and V ℎ of S and V (subscripts u and  are omitted for brevity), respectively; that is, where u   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 yield taking 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 as The following remarks might be stated regarding (61).Firstly, as isoparametric elements are considered, grad(  ) is computed by the following expression: where x = F , while ĝrad 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 3 × 3 matrix, representing the matrix of components of   , with a 3 × 1 vector grad(  ); that is, the strain operator B is unnecessary for such a multiplication.Now, replacing the approximation expression (54) 2 into (50) yields again, 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 div(k  ) in (65) can be computed as where the nodal values k   are obtained from u   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 u   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 [u   ,    ]  with  = 1, 2, . .., 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 [u   , ]  .Given a small value TOL (usual ∼10 −10 ), the Newton-Raphson process finished whenever the value of an adequate norm of the vector [R L  , R M  ]  is less than TOL.As the present work is restricted to an isotropic linear elastic setting, convergence is achieved in the first iteration.
B 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 R L  /u (70) 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 C  , replacing approximation (55) 1 , and applying relation ( 57), (36) yields (71) The internal virtual work density in expression (61) can thus be written as where the right-hand side in (72) can be expressed as In expression (73), the second-order tensor C  {grad(  ), grad(  )} is defined by the relation [12]: for any vectors a, b, k, w, and fourth-order tensor C.
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: (76) Some remarks follow.Firstly, due to expression (76), no fourt-order tensors appear within the tangent stiffness block R L  /u   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, Δ  = C  : Δ, the tangent stiffness block R L  /u   will have the same structure as (75).However, in this case the expression for C  {grad(  ), grad(  )} must be modified in accordance with the constitutive relation [12].Thirdly, for a given material defined by Δ  = C : Δ, the implementation of the tangent block (75) requires only writing a subroutine that computes C{u, k} for any pair of vectors u and k.
The block R L  /   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 R M  /u   , let us consider (65).Due to the Darcy's law (33), the only term appearing in (65) depending on u  is div(k  ); therefore, By ( 66) and ( 51), (79) yields Finally, the block R M  /   is described.By (65), this last block can be rewritten as By ( 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 u   and nodal fluid pressure    might be obtained at time  =   + Δ.It is worth mentioning some remarks regarding equal order of interpolation for displacements u  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 B 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 B free approach, and an analytical solution is presented, eight node mixed quadrilateral finite elements will be considered to interpolate both displacements u  and fluid pressure   fields.

Numerical Application
In the previous sections, the B 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][23][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 B 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 u  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  (0, ) = 0,   (, ) = 0.  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: where and ℵ  = (( + 2)/)((1 + 2)/2) 2 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.(91) Regarding the permeability, three different values have been considered from a high permeability  = 10 −2 m/s to a moderate value  = 10 −4 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.
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 u  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  = 0.51; 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  = 10 −2 m/s,  = 10 −3 m/s, and  = 10 −4 m/s.
It is clear that the B 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  = 10 −4 m/s, where the steeper pressure gradient is observed, the numerical solution closely follows the semianalytical one.

Conclusions
In this work, the B 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 B or D matrices.Even though this aspect is crucial in the method, the main advantage of the B 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.

Figure 1 :
Figure 1: Kinematical basic assumption within the theory of porous media.

Figure 2 :
Figure 2: Geometry and boundary conditions of the investigated problem.

Figure 3 :
Figure 3: Geometry, boundary conditions, mesh, and loading path for the numerical simulation.

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Numerical versus semianalytical response for  = 10 −2 m/s.(a) Vertical displacement at the top boundary and (b) pore pressure distribution with depth at  = 0.5 s.