

HybridTrefftz stress and displacement elements for transient analysis of incompressible saturated porous mediaAbstractThe stress and displacement models of the hybridTrefftz finite element formulation are developed for the analysis of incompressible saturated porous media. They are used to model the frequency and time domain response of soft tissue specimens under twodimensional and axisymmetric testing conditions. The alternative stress and displacement models are complementary in terms of the supporting approximation criteria. The stress (displacement) model is derived from the direct approximation of the stress and pressure fields (the displacements in the solid and fluid phases) in the domain of the element. Two sets of numerical tests are performed; namely the tests implemented in frequency domain and those in time domain, for both twodimensional and axisymmetric problems. The tests on the frequency domain are designed to assess the basic aspects of the performance of the element. The tests on the time domain are obtained using a nonperiodic spectral decomposition time integration procedure.1. IntroductionThe biphasic theory of Mow et al. reported in 1980 [1], and known as the KLM (Kuei, Lai, Mow) biphasic theory, is based on the coupling of interstitial fluid flow and solid matrix deformation. The KLM biphasic theory models soft hydrated tissues, such as articular cartilage and the intervertebral disk, as two immiscible, incompressible phases. The governing equations of linear biphasic theory are mathematically equivalent to those for Biot’s theory of linear quasistatic poroelasticity in the case of incompressible constituents [2,3].Besides the hybridTrefftz method reported in [410] and summarized in this paper, there have been four methods previously reported in the literature on the formulation of biphasic finite elements: the penalty formulation of Suh et al. [11,12], the mixedpenalty formulation of Spilker and Maxian [13], the hybrid formulation of Vermilyea and Spilker [14,15] and the displacementpressure formulation of Wayne et al. [16] or the similar velocitypressure formulation of Oomens et al. [17]. Penalty and mixedpenalty formulations are based on imposing the continuity condition using a penalty parameter. The hybrid formulations are designed to satisfy exactly the momentum equation of the mixture (but not of each phase) by using adequate approximations for the stress and pressure fields. The displacement and velocitypressure formulations of the finite element method, typically developed from the virtual work principle, are based on the direct approximation of the displacement or velocity fields in the solid matrix and of the pressure field in the fluid phase. Their values at the nodes of the element are the finite element degreesoffreedom. In hybridTrefftz formulations, the use of a Trefftz approximation basis ensures that all domain conditions of the problem are locally satisfied, namely the equilibrium (in each phase), the constitutive and the straindisplacement relations and the mixture incompressibility condition. The stress model develops from the direct (and naturally hierarchical) approximation of the stress and pressure fields in the domain of the element, while the displacements of the solid phase and the normal displacement of the fluid phase are approximated independently on the boundary of the element and used to enforce in weak form the interelement and boundary equilibrium conditions on the forces in the solid phase and on the fluid pressure. In the alternative displacement model, the fields approximated directly in the domain of the element are the displacements in each phase of the mixture and the boundary approximation is on the solid surface forces and the fluid pressure, which is used to enforce on average the interelement and boundary displacement continuity conditions. The domain approximation bases are extracted from the formal solution of the system of differential equations that govern the behaviour of incompressible saturated porous media. 2. Basic equationsLet V and Γ denote the domain and the boundary of the element, which may not be convex, simply connected or bounded. Five distinct regions are identified on the boundary of the element, Γ, namely the interelement boundary, Γ_{i}, the edges (or surfaces, in 3D applications) whereon forces or displacement are prescribed on the solid, Γ_{t} and Γ_{u}, respectively, and the edges whereon the pressure or the outward displacement is prescribed on the fluid, Γ_{p} and Γ_{w}, respectively. After implementation of the time integration method, see [18], the system of equations governing the response of a saturated porous element subject to a typical forcing frequency ω is stated as follows, where i is the imaginary unit and subscripts s and f identify quantities associated with the solid and fluid phases of the mixture, respectively:
In the Neumann boundary conditions (6) and (7), is the pressure prescribed on the fluid and matrix N collects the components of the unit outward normal vector, n, that extract the contribution of the stress vector components in equilibrium with the prescribed forces . In the Dirichlet conditions (8) and (9), vector defines the displacements prescribed on the solid matrix and is the outward normal component of the displacement in the fluid. The equations above can be written to account for mixed boundary conditions, which may be written in terms of velocity components. They hold for interelement boundary conditions, in which case the prescribed term defines the displacements and the (reaction) forces on the connecting element. In addition, the elasticity condition can be extended to include creep and stress relaxation terms. It is noted that the definition of the Neumann boundary is extended to contain the interelement boundary of the stress element defined below. Therefore, the prescribed terms in conditions (6) and (7) are extended to include the boundary force on connecting elements. Conversely, the Dirichlet boundary of the alternative displacement element contains its interelement boundary, meaning that the identification of the prescribed terms in conditions (8) and (9) is extended to include the boundary displacements of connecting elements. 3. Trefftz basesCombination of equations (1) to (5) leads to the mixed (NavierBeltrami) system of differential equations,
The homogeneous system (10) has three sets of solutions, namely constant pressure modes,
The Trefftz method, as applied in the present finite element context, develops from the condition of using approximation bases extracted from the solution set of the homogeneous system (10), corrected by particular solution terms associated with the prescribed body forces and the initial conditions. In order to exploit fully the hierarchical nature of the potential functions used to set up the approximation bases, it is convenient to express the estimates for the different fields as direct linear combinations of the independent solutions of the homogeneous differential equation, which can be extended to include terms associated with the (prescribed) particular solution. Thus, in the following description of the Trefftz approximation on the displacement, stress and strain fields,
The potential functions present in the general expressions (11) to (13), and the associated pressure, P, stress, S, strain, E, and solid and fluid displacement modes, U and W, respectively, are defined in [10]. Consequent upon the Trefftz constraint, the approximation modes satisfy locally the domain equations (1)(5), to yield the following relations on the complementary solution terms:
4. Alternative finite element modelsTwo basic models can be established in the formulation of the finite element method. The displacement model is designed to produce (weak) kinematically admissible solutions, that is, solutions that satisfy on average the domain and boundary conditions (3), (4), (8) and (9). The alternative stress model is designed to produce (weak) statically admissible solutions, which satisfy on average the domain and boundary conditions (1), (2), (6) and (7).The Trefftz approximations defined above satisfy the domain equilibrium and compatibility conditions in strong form, but they are not constrained to satisfy locally either the Neumann or the Dirichlet boundary conditions. When the hierarchical nature of the Trefftz basis is to be exploited through adaptive refinement, it is in general not feasible, or computationally effective, to combine the Trefftz solutions so to enforce locally the boundary and interelement continuity conditions. This leads to the independent approximation of boundary fields, namely the surface force on the solid and the fluid pressure on the Dirichlet boundary of the displacement element, and the boundary displacements on the solid and the outward normal fluid displacement on the Neumann boundary of the alternative stress element: thus the hybrid nature of the formulation used here. 4.1 Displacement modelThe displacement model of the hybridTrefftz finite element formulation is derived approximating directly the solid and fluid displacement fields in the domain of the element, as stated by equations (14) and (15), respectively, and the force and pressure fields on the Dirichlet boundary, that is, the boundaries of the element whereon the surface forces are unknown. This independent boundary approximation is stated as follows,
The primary approximations (14) and (15) satisfy locally the mixture incompressibility condition (4) through constraint (23) and its particular solution term (28). In order to ensure the local enforcement of the compatibility condition (3), as it is typical of displacement elements, a compatible strain field (18) is assumed as a dependent approximation. 4.1.1 Dual variables As it is shown below, the generalized body forces defined by the dual transformation of the displacement approximation,
Similarly, the generalized stress defined by the dual transformation of the dependent strain approximation (18),
4.1.2 Equilibrium condition The displacement element equilibrium equation represents the weak form of the static admissibility conditions (1), (2), (6) and (7):
Implementation of these operations yields the following expression for the (Hermitian) damping matrix, , and for vector , which collects the equivalent velocity contribution from the particular solution:
4.1.3 Compatibility and elasticity conditions As conditions (3) and (4) are locally satisfied by the dependent approximations (18) and (19), the finite element kinematic admissibility condition reduces to the average enforcement of the Dirichlet conditions (8) and (9) for the assumed displacements (14) and (15). The following expression, the dual of the equilibrium condition (36),
The finite element elasticity equation is obtained inserting the local condition (5) in definition (33) for the generalised effective stress vector and enforcing next the dependent strain approximations (18) and (19),
4.1.4 Solving sytem The solving system for the hybridTrefftz displacement model is obtained combining the equilibrium, compatibility and elasticity equations (36), (41) and (42), respectively:
The solving system for the assembled finite element mesh preserves the same structure. The setting up of this system requires direct allocation of the elementary matrices and vectors, without the summation operations that characterize the assemblage of the solving system of conventional elements. As the basis weighting vector, z, is strictly element dependent, matrix D is blockdiagonal. The generalized force and pressure vectors, f_{t} and f_{p}, respectively, are shared by at most two connecting elements. 4.2 Stress modelThe stress model of the hybridTrefftz finite element formulation is derived approximating directly the stress and pressure fields in the domain of the element, as stated by equations (16) and (17), respectively, and the solid displacements and the outward normal component of the displacement in the fluid on its Neumann boundary,
This dependent displacement approximation may not be associated with compatible strains, as required by equation (22). This condition, together with the elasticity relation (24), is called upon only to establish the boundary integral description of the finite element elasticity condition, the feature that distinguishes the hybrid and hybridTrefftz variants of the stress model of the finite element method. 4.2.1 Dual variables The generalized strain, defined by the dual transformation of the stress approximations (16) and (17),
Similarly, the generalized boundary forces defined by the dual transformations of approximations (50) and (51) are used to enforce on average, still in the sense of Galerkin, the Neumann boundary conditions (6) and (7):
4.2.2 Compatibility condition The stress element compatibility equation represents the weak form of the kinematic admissibility conditions (3), (4), (8) and (9),
Implementation of these operations recovers definitions (37) and (38) for the damping terms, and yields the following expressions for the boundary compatibility matrices and for the terms defining the effect of the prescribed displacements:
4.2.3 Equilibrium and elasticity conditions As conditions (1) and (2) are locally satisfied, the finite element equilibrium condition reduces to the average enforcement of the Neumann conditions (6) and (7). The following expression, the dual of the compatibility condition (55),
The finite element elasticity equation is obtained inserting the flexibility format of the local condition (5) in definition (52) for the generalised strain vector, under the incompressibility condition, and enforcing next the stress approximations (16), (17), to yield,
4.2.4 Solving system The solving system for the hybrid stress model is obtained combining the compatibility, equilibrium and elasticity equations (55), (58) and (59), respectively:
5. Numerical applicationsThe benchmark tests used in the literature, namely the unconfined indentation and the compression of cartilage specimens under confined and unconfined conditions, are solved in the frequency domain to asses the convergence characteristics of the alternative stress and displacement elements. They are implemented also in the time domain to asses the quality of the estimates they produce for the stress, pressure, displacement and velocity fields at every instant of the analysis [10].
The confined compression tests for plain strain and axisymmetric applications are represented in Figure 1. The confining chamber is assumed to be rigid, impermeable and lubricated, to yield boundary conditions (65) and (66) for twodimensional and axisymmetric problems, respectively,
5.1 Frequency domain testsThe displacement driven loading program is used in the assessment of the convergence of the finite element solutions for confined compression tests, under boundary conditions (67) and (68) for twodimensional and axisymmetric problems, respectively. They are implemented on the singleelement mesh and the (regular, unbiased) 2×2 and 4×4element meshes for an increasing refinement of the approximation bases, under the combination d_{V}=2d_{Γ}+1, where d_{V} is the order of the domain approximation and d_{Γ} is the degree of the boundary approximation.The results obtained for twodimensional and axisymmetric problems are presented in Figures 2 and 3, respectively, where N defines the dimension of the solving systems (49) and (64). The strong convergence under prefinement is typical of hybridTrefftz elements and, for coarse meshes, the rate of convergence under hrefinement is directly affected by the boundary conditions of the problem, which affect differently the alternative stress and displacement models. In the twodimensional case, the rate of convergence is also affected by the ratio between the typical length of the element and the wavelength, which in the tests shown varies from 0.7 (4×4element mesh) to 3.1 (singleelement mesh), for the forcing frequency ω=0.10 rad/s. The reference values for the strain energy norm correspond to stabilized solutions obtained with boundary approximations of order (d_{V};d_{Γ})=(19;9) for both the twodimensional and axisymmetric problems.
5.2 Time domain testsThe modelling of the response of the cartilage specimen under the confined compression loading defined by conditions (65) and (67) is used to illustrate the variation of the results in time under twodimensional testing conditions and by conditions (66) and (68) under axisymmetric testing conditions.The finite element estimates for the displacement, stress and pressure fields are determined directly from approximations (14) to (17) for each forcing frequency generated by the nonperiodic spectral decomposition method described in [18]. No smoothing is used in the representation of the (double precision) finite element solutions presented below. The displacementdriven, rampload program shown in Figure 4 is designed to reach a deformation of 5% at instant t_{0}=500s, which corresponds to a prescribed displacement . The time domain tests are implemented in a single time step: Δt=t_{max}=1'000s. The results are obtained using a nonperiodic wavelet basis defined on the (unit) interval [20,21], with dimension N=256 (family 3 and refinement 7). This highorder but rather stable basis is implemented on a single time step, as the support Δt of the time basis T_{n}(t) is identified with the full duration of the test. No use is made of the adaptive refinement capabilities of wavelet bases.
The confined compression tests on twodimensional and axisymmetric problems are modelled with meshes of 2×2 and 4×1 elements, respectively. The approximation used is (d_{V};d_{Γ})=(11;5). The solutions shown here are extracted from the animations that can be accessed on tomamil.eu. The variation in time of the results is shown in Figures 5 and 6 for twodimensional and axisymmetric problems, respectively. They recover the analytic solution [1] and the values reported in Vermilyea and Spilker [14], obtained with a mesh with 5×2 pairs of hybrid elements and a trapezoidal time integration rule, implemented with Δt=5s (200 time steps). The coordinates of the control point are measured in the system of reference defined in Figure 1. The displacement and stress models produce basically identical results. For the confined compression tests, the analytic results reported in [1] predict that the fluid should flow upwards uniformly. Indeed, the results are uniform at any radial location. 6. ClosureThe derivation of Trefftz elements is often muddled by the fact that the domain approximation basis locally satisfies all domain conditions of the problem. This source of confusion is strengthened by the existence of two alternative models and the multiple fields involved in the modelling of the response of incompressible soft tissues.The option followed here is to derive first the hybrid formulation, where the constraints on the domain basis are weaker and better focused: the only constraint set a priori in the development of the stress (displacement) model is to ensure that the domain equilibrium (compatibility) condition is locally satisfied. The need to complement the domain approximation with an independent boundary approximation follows naturally from the derivation of the finite element equations, namely boundary forces and boundary displacements for the stress and displacement models. Duality, the vector representation of the invariance of the inner product in all finite element mappings (approximations), is used consistently to support the implementation of the approximation criteria adopted in the development of the alternative stress and displacement models. This derivation shows that the alternative displacement and stress models of the hybridTrefftz formulation are based on the direct approximation of the displacements in the solid and fluid phases and the stress and pressure fields, respectively, in the domain of the element. In the displacement (stress) model, the boundary and the interelement displacement (force) continuity conditions are enforced explicitly in both the fluid and solid phases of the mixture. This approach leads to the independent definition of the weak forms of static and kinematic admissibility, which are related with the weak implementation of the constitutive relations. It is noted that the mixture incompressibility condition is locally respected in the derivation of the displacement model and enforced on average in the alternative stress model. The finite element solving systems are derived by direct combination of the basic finite element equations. These systems are sparse, wellsuited to adaptive refinement and parallel processing. This paper presents only very small part of the results for hybridTrefftz stress and displacement elements for transient analysis of incompressible saturated porous media. They are all present in [410]. As for twodimensional problems, it is shown that stable accurate solutions for each phase of the medium can be obtained using relatively coarse finite element meshes and a low number of degreesoffreedom. The tests reported here confirm, again, the experience gained in elastostatics concerning the role played by the boundary approximation. The frequency domain tests are used to illustrate the patterns and rates of convergence that can be attained under both p and hrefinement procedures. They show that the element is basically insensitive to gross distortion in shape. The tests show, also, that the element remains stable when each phase of the medium approaches the incompressibility limit. The time domain tests, that are reported here, are captured by coupling the finite element model with a highorder wavelet time approximation basis, which is implemented in single time step. Although the adaptive refinement features of this basis are not exploited in this implementation reported here, the time integration procedure performed well and proved its ability to capture the peaks in acceleration that characterize the benchmarks tests on finite element modelling of the response of hydrated soft tissues. Due to the nature of the cylindrical description of the Trefftz solutions, the modelling of axisymmetric hydrated soft tissue specimens requires the implementation of elements with low widthtoheight ratios. Consequent upon this weakness, refinement of the discretization in the zdirection may be needed to solve adequately particular applications. References
