# B Free Finite Element Approach for Saturated Porous Media: Consolidation

Authors: M. M. Stickle; A. Yague; M. Pastor
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

TheB 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.


---

## Body

## 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, theB 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 withinB 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 operatorB, 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 theB 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.TheB 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 u-p 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 fourth-order tensors will never appear explicitly, irrespective of the material type considered.An important point is whether theB 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 Section2, 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.


## 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 particlesPαα=f,s. The superscript α denotes pore fluid for α=f and solid skeleton for α=s.The kinematics in TPM is based on the assumption [4] that the spatial position x in the current configuration at time t is simultaneously occupied by particles Pα of both constituents (α=f,s). Each particle Pαα=f,s proceeds from its own reference position Xα in the reference configuration at time t0; therefore, each phase has its own motion χα and at the common spatial position x the phases interact (Figure 1).Figure 1
Kinematical basic assumption within the theory of porous media.The following relation holds at timet:(1)x=χαXα,t.An important concept within the modern TPM, responsible for the notion of superimpose phases over the control space, is the volume fractionnα(x,t) of the phase α, defined by(2)nαx,t=dvαdv,where dvα is the element volume of the phase α in the current configuration within a representative element volume dv of the overall porous media. The volume fractions nα(x,t) in (2) satisfy the following relation:(3)ns+nf=1.Partial densitiesρα of each phase (α=f,s) are defined as the averaged density with the volume fraction nα; that is,(4)ρα=nαραR,where ραR 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,(5)ραR=constant.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α: (6)DαραDt+ραdivvα=0.Balance of momentum of phaseα: (7)divσα+ραb-aα+rα=0.In expression (6), vα is the velocity of the particles Pα defined by(8)vα=∂χαXα,t∂t,while Dα□/Dt is the material derivative with respect to the phase α defined by(9)Dα□Dt=∂□∂t+grad⁡□vα.In expression (7), σα is the partial Cauchy stress second-order tensor of the phase α, b is the gravity acceleration vector, and aα is the acceleration vector of the particles Pα 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:(10)rs+rf=0.The symbolgrad⁡(□) 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(11)DαnαDt+nαdiv⁡vα=0.Now, by (3), (9), and (11), the volume balance of the mixture is obtained:(12)div⁡nfvf+nsvs=0.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 vfs of the fluid with respect to the velocity of the solid is introduced:(13)vfs=vf-vs.Replacing (13) into (12), keeping (3) in mind, yields(14)divvs+nfvfs=0.As a consequence of the incompressibility constraint (5) for α=s, the partial stress tensor for the solid σs can be additively decomposed into two terms:(15)σs=σ′-nspfI,where σ′ is the Cauchy effective stress second-order tensor, I is the second-order identity tensor, and pf is the pore fluid pressure. Due to assumption (4), the partial stress tensor σf associated with this phase can be taken as(16)σf=-nfpfI.By adding the partial stresses (15) and (16), the total Cauchy stress second-order tensor is recovered via Terzaghi’s principle:(17)σ=σ′-pfI.Replacing (15), (16) in (7) with α=s and α=f, we arrive to the expressions(18)divσ′-grad⁡nspf+ρsb-as+rs=0,(19)-grad⁡nfpf+ρfb-af+rf=0.The volume fractionns in (18), and by (3) also nf in (19), can be computed by solving analytically the ordinary differential equation (11) with α=s. Assuming a small strain regime, this yields(20)ns=n0s1-div⁡us.However, sincediv(us)≪1 in a small strain setting, ns may be approximated by n0s, 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 α=f and α=s and for the Cauchy effective stress tensor σ′.Concerning the interaction vectorrf, the relative velocity defined in (13) allows its definition as explained in [4, 11, 18]. In the case of isotropic permeability, the vector rf can be expressed as(21)rf=pfgrad⁡nf-nf2gρfRkfvfs,where g is the gravity acceleration and kf is the Darcy permeability coefficient of the porous medium. Vector rs is obtained from (21) through (10). Replacing these expressions for the interaction force vectors in (18) and (19) yields(22)div⁡σ′-nsgrad⁡pf+ρsb-as+nf2gρfRkfvfs=0,-nfgrad⁡pf+ρfb-af-nf2gρfRkfvfs=0.Regarding the constitutive relation for the Cauchy effective stress tensorσ′, it is convenient to introduce the solid skeleton displacement vector us defined by(23)us=x-Xs.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:(24)σ′=Ce:ε,where(25)Ce=λI⊗I+2μIsis the fourth-order elastic tensor, while ε=1/2(grad⁡us+gradT(us)) is the small strain tensor. In (25), I is the second-order identity tensor and Is 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(26)div⁡σ′-nsgrad⁡pf-ρsas+nf2gρfRkfvfs=0,-nfgrad⁡pf-ρfaf-nf2gρfRkfvfs=0.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(27)div⁡σ′-nsgrad⁡pf+nf2gρfRkfvfs=0,(28)-nfgrad⁡pf-nf2gρfRkfvfs=0.Adding (27) and (28), solving for vfs in (28), and taking into account relation (17) yield(29)div⁡σ=0,(30)nfvfs=-kfρfRggrad⁡pf,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 us and pf 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:(31)div⁡σ=0.Mixture volume balance:(32)div⁡vs+nfvfs=0together with the following equations to close the system:Darcy’s law:(33)nfvfs=-kfρfRggrad⁡pf.Constant fluid volume fraction:(34)nf=1-n0s.Terzaghi’s principle:(35)σ=σ′-pfI.Solid constitutive relation:(36)σ′=Ce:ε.Solid strain tensor:(37)ε=12grad⁡us+gradT⁡us.Solid displacement vector:(38)us=x-Xs.Solid velocity vector:(39)vs=u˙s,where the overdot in (39) indicates partial derivative with respect time t. The primary variables for this system are us and pf. 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]: (40)∂Ω=∂Ωus∪∂Ωt¯=∂Ωpf∪∂Ωq¯,where the overline denotes closure while ∂Ωus and ∂Ωt are nonintersecting parts of the total boundary ∂Ω on which solid displacement and total tractions are prescribed, respectively. On the other hand, ∂Ωpf and ∂Ωq 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:(41)us=u^son∂Ωus.Total traction specified:(42)σn=t^on∂Ωt.Pore fluid pressure specified:(43)pf=p^fon∂Ωpf.Fluid flux specified:(44)nfvfs·n=q^on∂Ωq,where u^ 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 nfvfs is the Darcy velocity.The trial solution spacesSp(t) and Su(t) and the variation spaces Vp and Vu that are needed to define the weak or variational formulation of the problem can be described as follows:(45)Sut=us·,t:Ω⟶R3∣us·,t∈H1,us·,t=u^s·,ton∂Ωus,(46)Spt=pf·,t:Ω⟶R∣pf·,t∈H1,pf·,t=p^f·,ton∂Ωpf,(47)Vu=φ:Ω⟶R3∣φ∈H1,φ=0on∂Ωus,(48)Vp=φ:Ω⟶R∣φ∈H1,φ=0on∂Ωpf,where H1 and H1 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:(49)∫∂Ωtφ·t^d∂Ω-∫Ωσ:grad⁡φdΩ=0,(50)∫Ωφdiv⁡vsdΩ+∫∂Ωqφq^d∂Ω-∫Ωgrad⁡φ·nfvfsdΩ=0.In standard finite element procedure, thegrad⁡(φ) in (49) is usually changed by grads(φ), 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 us and pf. With this objective, it is assumed that uns and pnf are known at time tn and determination of us and pf at time t=tn+Δt is wanted (subscript “n+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 t=tn+Δt. 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 vs at time t=tn+Δt is considered:(51)vs=α1us-α1uns+α2vns,where α1=1/θΔt and α2=1-θ/θ, while θ∈[0,1]. If θ⩾1/2, an unconditionally stable method is obtained. Also in (51), uns and vns are known values of us and vs at time tn, respectively. Substituting (51) into (50), the variational problem reads as follows:(52)Findus∈Sut,pf∈Sptat  timets.t.L=M=0,∀φ∈Vp,∀φ∈Vu,where Su(t), Sp(t), with t=tn+Δt, are defined in (45) and (46), respectively. Spaces Vu and Vp are defined in (47) and (48), respectively, while L and M are defined by(53)L=∫∂Ωtφ·t^d∂Ω-∫Ωσ:grad⁡φdΩ,M=∫Ωφdiv⁡vsdΩ+∫∂Ωqφq^d∂Ω-∫Ωgrad⁡φ·nfvfsdΩ.For the solution of the variational problem (52), at a fixed time t=tn+Δt, the standard Bubnov-Galerkin approach is applied. In the following section, the B 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 Ne elements Ωe, connecting Nnd nodes. Mixed finite element formulation with equal order of interpolation for displacements us and fluid pressure pf is employed. In this way, each node has nsd⁡+1 degrees of freedom, where nsd⁡=spatial dimension. The independent variables are the vector of nodal displacements uαs and nodal fluid pressure pαf, where α∈1,…,Nnd.Isoparametric shape functionsNα:Ωr→R are defined in standard fashion where Ωr is the reference element associated with the local orthogonal coordinate system ξ: that is, for every x∈Ωe, there is one and only one ξ∈Ωr such that x=x(ξ). These shape functions span the finite dimensional subspaces Sh and Vh of S and V (subscripts u and p are omitted for brevity), respectively; that is, (54)φ≈∑α=1NndNαξφα,φ≈∑α=1NndNαξφα,(55)us≈∑α=1NndNαξuαs,pf≈∑α=1NndNαξpαf,where uαs and pαf are nodal values of the trial functions, already evaluated at a fixed time t=tn+Δt, 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(56)∫∂ΩtNαφα·t^d∂Ω-∫Ωσ:grad⁡NαφαdΩ=0,taking into account the following relations from tensor calculus:(57)grad⁡Nαφα=φα⊗grad⁡Nα,(58)σ:φα⊗grad⁡Nα=φα·σgrad⁡Nα,where ⊗ stands for the standard dyadic product operation, and, by bilinearity of the scalar product, (56) reads as(59)∫∂ΩtNαt^d∂Ω-∫Ωσgrad⁡NαdΩ·φα=0by noting that the nodal values φα of the weighting functions are arbitrary, it follows(60)∫∂ΩtNαt^d∂Ω-∫Ωσgrad⁡NαdΩ=0.Making use of Terzaghi’s principle (35) yields the first set of nodal equilibrium equations:(61)RαL=0s.t.RαL=∫∂ΩtNαt^d∂Ω-∫Ωσ′grad⁡NαdΩ+∫Ωpfgrad⁡NαdΩ.The following remarks might be stated regarding (61). Firstly, as isoparametric elements are considered, grad⁡(Nα) is computed by the following expression:(62)grad⁡Nα=F^-Tgrad^Nα,where dx=F^dξ, while grad^ 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⁡(Nα); that is, the strain operator B is unnecessary for such a multiplication.Now, replacing the approximation expression (54)2 into (50) yields(63)∫ΩNαφαdiv⁡vsdΩ+∫∂ΩqNαφαq^d∂Ω-∫Ωgrad⁡Nαφα·nfvfsdΩ=0again, recalling from tensor calculus the following relation:(64)grad⁡Nαφα=grad⁡Nαφα.The second set of nodal equilibrium equations is obtained as(65)RαM=0s.t.RαM=∫ΩNαdiv⁡vsdΩ+∫∂ΩqNαq^d∂Ω-∫Ωgrad⁡Nα·nfvfsdΩ.Again, some aspects need to be clarified. Firstly, keeping in mind expressions (55)1 and (39), the term div⁡(vs) in (65) can be computed as(66)div⁡vs=vαs·grad⁡Nα,where the nodal values vαs are obtained from uαs via (51). Secondly, Darcy’s law (33) and relation (64) yield(67)nfvfs=-kfρfRggrad⁡Nαpαf.Simultaneous solution of algebraic system of (61) and (65) for the unknown nodal displacements uαs and nodal fluid pressure pαf 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αs,pαfi with i=1,2,…, of displacement and pressure values at node α by the following iterative scheme:(68)∂RαL∂uβs∂RαL∂pβf∂RαM∂uβs∂RαM∂pβfiδuβsδpβf=-RαLRαMi,(69)uβspβfi+1=uβspβfi+δuβsδpβf.The matrix at the left-hand side in (68) is the tangent matrix of the coupled problem. The subscript i in (68) and (69) means evaluation at uβs,pβfi.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αMi 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αs, ∂RαL/∂pαf, ∂RαM/∂uαs, and ∂RαM/∂pαf in (68) should be specified.Starting by∂RαL/∂uβs, consider (61), with integration performed at element level; then,(70)∂RαL∂uβs=-∫Ωe∂σ′grad⁡Nα∂uβsdΩ.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 tensorCe, replacing approximation (55)1, and applying relation (57), (36) yields(71)σ′=Ce:ε=Ce:grad⁡us=Ce:uβs⊗grad⁡Nβ.The internal virtual work density in expression (61) can thus be written as(72)σ′grad⁡Nα=Ce:uβs⊗grad⁡Nβgrad⁡Nα,where the right-hand side in (72) can be expressed as(73)Ce:uβs⊗grad⁡Nβgrad⁡Nα=Cegrad⁡Nα,grad⁡Nβuβs.In expression (73), the second-order tensor Cegrad⁡Nα,grad⁡Nβ is defined by the relation [12]:(74)a·Cv,wb=a⊗v:Cb⊗wfor any vectors a, b, v, w, and fourth-order tensor C. Replacing now the right-hand side of (73) in (70) yields(75)
∂
R
α
L
∂
u
β
s
=
-
∫
Ω
e
C
e
grad
⁡
N
α
,
grad
⁡
N
β
d
Ω
.Taking into account (74) and keeping in mind the definition of the fourth-order elastic tensor (25), the following expression is obtained:(76)Cegrad⁡Nα,grad⁡Nβ=λgrad⁡Nα⊗grad⁡Nβ+μgrad⁡Nα·grad⁡NβI+μgrad⁡Nβ⊗grad⁡Nα.Some remarks follow. Firstly, due to expression (76), no fourt-order tensors appear within the tangent stiffness block ∂RαL/∂uβs 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, Δσ′=Cep:Δε, the tangent stiffness block ∂RαL/∂uβs will have the same structure as (75). However, in this case the expression for Cepgrad⁡Nα,grad⁡Nβ 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 Cu,v for any pair of vectors u and v.The block∂RαL/∂pβf of the tangent matrix (68) will be described next. Again, by (61), the following expression, at element level, is obtained:(77)∂RαL∂pβf=∫Ωe∂pfgrad⁡Nα∂pβfdΩ.Taking into account the approximation (55)2, expression (77) becomes(78)
∂
R
α
L
∂
p
β
f
=
∫
Ω
e
grad
⁡
N
α
N
β
d
Ω
.For the following block∂RαM/∂uβs, let us consider (65). Due to the Darcy’s law (33), the only term appearing in (65) depending on us is div(vs); therefore, (79)∂RαM∂uβs=∫Ωe∂Nαdiv⁡vs∂uβsdΩ.By (66) and (51), (79) yields(80)
∂
R
α
M
∂
u
β
s
=
∫
Ω
e
N
α
α
1
grad
⁡
N
β
d
Ω
.Finally, the block∂RαM/∂pβf is described. By (65), this last block can be rewritten as(81)∂RαM∂pβf=-∫Ωe∂grad⁡Nα·nfvfs∂pβfdΩ.By (67), (81) becomes(82)
∂
R
α
M
∂
p
β
f
=
∫
Ω
e
k
f
ρ
f
R
g
grad
⁡
N
α
·
grad
⁡
N
β
d
Ω
.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αs and nodal fluid pressure pαf might be obtained at time t=tn+Δt.It is worth mentioning some remarks regarding equal order of interpolation for displacementsus and fluid pressure pf employed in the present work. Firstly, it is well known that, in the incompressible and nearly incompressible regimes, fluid phase incompressible ((5) with α=f) 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 us and fluid pressure pf fields.


## 5. Numerical Application

In the previous sections, theB 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 L, composed by a fluid saturated poroelastic material and subjected to a time dependent loading ft. 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 theB 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 thez axis, which is considered normal to the top boundary and pointing downward (Figure 2).Figure 2
Geometry and boundary conditions of the investigated problem.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(83)λ+2μuzz-pz=0,utz-kρgpzz=0.In (83), the vertical value of us has been substituted by u, the fluid pressure pf by p, the Darcy’s permeability kf by k, and the intrinsic density of the fluid phase ρfR by ρ. Subscript z, like in pz or in uzz, means partial differentiation with respect to the spatial coordinate z, while subscript t, for example in utz, 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 forp are (84)p0,t=0,pzL,t=0.Boundary conditions foru are (85)σ0,t=ft,uL,t=0.Initial condition foru is (86)uz,0=0.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:(87)uz,t=-ftLλ+2μ1-zL+∑n=0∞Tntcos⁡1+2nπz2L,pz,t=ρgf′tz2-2Lz2kλ+2μ+∑n=0∞2LρgTn′t1+2nπksin⁡1+2nπz2L,where(88)Tnt=e-ℵntTn0+∫0te-ℵnt-τhnτdτ,(89)Tn0=8Lf0π22n+12λ+2μ,(90)hnτ=8Lf′τπ22n+12λ+2μ,and ℵn=kλ+2μ/ρg1+2nπ/2L2 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 gravityg=10 m/s2, length of the column L=10 m, a harmonic time dependent loading ft=1-cos⁡2πt kN/m2, and the following values for the material parameters:(91)λ=8.375MN/m2,μ=5.583MN/m2,ρ=1000kg/m3.Regarding the permeability, three different values have been considered from a high permeabilityk=10-2 m/s to a moderate value k=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 Figure3.Figure 3
Geometry, boundary conditions, mesh, and loading path for the numerical simulation.
(a) (b)As it can be observed in Figure3, eight node mixed quadrilateral finite elements are considered. This type of elements is employed to interpolate both displacements us and fluid pressure pf 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 Δt = 5·10−2 s is considered.In Figures4, 5, and 6, the vertical displacement at the top boundary and the fluid pressure distribution along the column can be observed for permeabilities k=10-2 m/s, k=10-3 m/s, and k=10-4 m/s.Figure 4
Numerical versus semianalytical response fork=10-2 m/s. (a) Vertical displacement at the top boundary and (b) pore pressure distribution with depth at t=0.5 s.
(a) (b)Figure 5
Numerical versus semianalytical response fork=10-3 m/s. (a) Vertical displacement at the top boundary and (b) pore pressure distribution with depth at t=0.5 s.
(a) (b)Figure 6
Numerical versus semianalytical response fork=10-4 m/s. (a) Vertical displacement at the top boundary and (b) pore pressure distribution with depth at t=0.5 s.
(a) (b)It is clear that theB 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 withk=10-4 m/s, where the steeper pressure gradient is observed, the numerical solution closely follows the semianalytical one.


## 6. Conclusions

In this work, theB 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-calledB 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.



---

*Source: 4256079-2016-06-05.xml*