Technical University of Lisbon, Instituto Superior Tecnico, Department of Civil Engineering and Architecture, Av. Rovisco Pais, 1049-001, Lisboa, Portugal

Abstract

The stress and displacement models of the hybrid-Trefftz 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 two-dimensional 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 two-dimensional 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 non-periodic spectral decomposition time integration procedure.

1. Introduction

The 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 quasi-static poroelasticity in the case of incompressible constituents [2,3].

Besides the hybrid-Trefftz method reported in [4-10] 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 mixed-penalty formulation of Spilker and Maxian [13], the hybrid formulation of Vermilyea and Spilker [14,15] and the displacement-pressure formulation of Wayne et al. [16] or the similar velocity-pressure formulation of Oomens et al. [17]. Penalty and mixed-penalty 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 velocity-pressure 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 degrees-of-freedom.

In hybrid-Trefftz 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 strain-displacement 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 inter-element 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 inter-element 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 equations

Let 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 inter-element 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:

		(1)
		(2)
		(3)
		(4)
		(5)
		(6)
		(7)
		(8)
		(9)

Vectors σ and ε collect the independent components of the stress and strain tensors, p and γ are the pressure and the volumetric change of the mixture, b and u are the body force and displacement vectors, respectively. Vector b_o is used to define the equivalent body forces associated with the initial condition of the problem, and parameters ζ, Φ_s and Φ_f define the diffusive drag and the solid and fluid fraction ratios, with Φ_s + Φ_f = 1. The (linear) differential operators, the divergence D and the gradient , present in the equilibrium equations (1), (2) and their conjugates present in the compatibility equation (3) and continuity equation (4), D^* and ^*, are defined elsewhere for different applications, e.g. in [10], as well as the terms present in the elasticity condition (5).

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 inter-element 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 inter-element 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 inter-element 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 bases

Combination of equations (1) to (5) leads to the mixed (Navier-Beltrami) system of differential equations,

(10)

where is the anti-gradient vector and represents its conjugate: . They are defined in [10], where the P- and S-wavenumbers k_p and k_s are given in terms of parameter .

The homogeneous system (10) has three sets of solutions, namely constant pressure modes,

(11)

harmonic pressure modes,

(12)

and Helmholtz pressure modes,

(13)

where is the Laplacian and its conjugate.

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,

		(14)
		(15)
		(16)
		(17)
		(18)
		(19)

vector z collects the weights of the strain inducing solution modes, and denotes particular solution terms. These (non-nodal) finite element unknowns, need not, and in general will not have unequivocal physical interpretations.

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:

		(20)
		(21)
		(22)
		(23)
		(24)

Similar relations can be stated for the particular solution terms:

		(25)
		(26)
		(27)
		(28)
		(29)

4. Alternative finite element models

Two 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 inter-element 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 model

The displacement model of the hybrid-Trefftz 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,

		(30)
		(31)

where matrices Z collect the boundary approximation functions, typically Legendre or Chebyshev polynomials, and the weighting vectors f collect generalized forces on the solid and fluid phases.

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,

(32)

are used to enforce on average, in the sense of Galerkin, the equilibrium conditions (1) and (2) for the Neumann conditions (6) and (7). In the notation used here, defines the (complex) conjugate of the transpose of array : .

Similarly, the generalized stress defined by the dual transformation of the dependent strain approximation (18),

(33)

is used to enforce on average the constitutive relation (5) and the generalized boundary displacements defined by the dual transformations of approximations (30) and (31) are used to enforce on average the Dirichlet conditions (8) and (9), still in the sense of Galerkin, in terms of average displacements and :

		(34)
		(35)

4.1.2 Equilibrium condition

The displacement element equilibrium equation represents the weak form of the static admissibility conditions (1), (2), (6) and (7):

(36)

This equation is derived as follows: the local equilibrium conditions (1) and (2) is inserted in definition (32) for the generalized body forces; the resulting weak form of the domain equilibrium condition is integrated by parts to mobilize the boundary terms; the displacement approximations (14) and (15), under constraints (22), (23), and definition (33) for the generalized stresses are enforced in the domain integral terms of the equation; its boundary term is uncoupled to enforce conditions (6) and (7) on the Neumann boundary of the element and approximations (30) and (31) on its Dirichlet boundary.

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:

		(37)
		(38)

The boundary equilibrium matrices and the (average) prescribed boundary force vector are defined as follows:

		(39)
		(40)

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),

		(41)

is found inserting the Dirichlet conditions in definitions (34) and (35) for the generalized boundary displacements, enforcing the displacement approximations (14) and (15), and recalling results (39) and (40).

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),

(42)

where the (Hermitian) stiffness matrix and the associated particular solution vector are defined by:

		(43)
		(44)

Integrating by parts the definition (42) for the stiffness matrix and enforcing the Trefftz constraints (20) to (24), the following alternative expression is found with the damping matrix (37):

(45)

where the (in general, non-Hermitian) dynamic matrix has the following boundary integral expression:

(46)

Definition (44) can be processed similarly, using now the particular solution relations (25) to (29). The following alternative expression for the contribution of the particular solution term, where results (32) and (38) hold, is obtained:

		(47)

It is possible, also, to replace definition (32) for the body force resultant by the following equivalent boundary integral expression:

(48)

4.1.4 Solving sytem

The solving system for the hybrid-Trefftz displacement model is obtained combining the equilibrium, compatibility and elasticity equations (36), (41) and (42), respectively:

(49)

It is noted that all terms present in system (49) have boundary integral expressions, as it is typical of Trefftz finite element formulations. Moreover, the system is Hermitian except for the damping term, as shown by definition (45).

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 block-diagonal. The generalized force and pressure vectors, f_t and f_p, respectively, are shared by at most two connecting elements.

4.2 Stress model

The stress model of the hybrid-Trefftz 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,

		(50)
		(51)

The development of hybrid stress models is based on the assumption that the equilibrium conditions are locally satisfied. In the context of time dependent problems [19], this constraint is fulfilled by assuming that the (primary) stress approximation is supported by an equilibrating (dependent) displacement approximation, as stated by equations (14), (15) and constraints (20), (21) and (25), (26) in the present application.

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 hybrid-Trefftz 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),

(52)

is used to enforce on average, in the sense of Galerkin, the compatibility conditions (3) and (4), and the constitutive relation (5). It is recalled that the condition on the null volumetric change of the mixture is locally fulfilled by the Trefftz constraint.

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):

		(53)
		(54)

4.2.2 Compatibility condition

The stress element compatibility equation represents the weak form of the kinematic admissibility conditions (3), (4), (8) and (9),

(55)

and can be recovered as follows: the local compatibility condition (3) and the continuity condition (4) is inserted in definition (52) for the generalized deformation; the resulting weak form of the domain compatibility condition is integrated by parts to mobilize the boundary terms; the dependent displacement approximation (14), (15) and constraint (20) are enforced in the domain integral terms of the equation; its boundary term is uncoupled to enforce conditions (8) and (9) on the Dirichlet boundary of the element and approximations (50) and (51) on its Neumann boundary. It is recalled that, in the notation used here, defines the complex conjugate of the forcing frequency, .

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:

		(56)
		(57)

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),

		(58)

is found inserting the Neumann conditions in definitions (53) and (54) for the generalized boundary forces, enforcing the stress approximations (16), (17), and recalling results (56), (57).

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,

(59)

where the (Hermitian) flexibility matrix and the associated particular solution vector are defined as follows, where :

		(60)
		(61)

Integrating by parts the definitions above and enforcing the Trefftz constraints, the following alternative expressions are found, in terms of results (37), (38) and (46):

		(62)
		(63)

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:

(64)

The solving system for the assembled finite element mesh presents the same structure and preserves the features identified for the displacement model.

5. Numerical applications

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

Figure 1: Confined compression tests for plain strain and axisymmetric applications.

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 two-dimensional and axisymmetric problems, respectively,

		(65)
		(66)

where subscript n identifies the displacement component normal to the boundary. A prescribed displacement is applied to the permeable, lubricated loading platen:

		(67)
		(68)

The testing data used in the time domain applications reported here is the following [13]: width w=6.35mm, height h=1.78mm, modulus of elasticity E=0.675MPa, Poisson's ratio ν=0.125, permeability κ=7.6×10^-15m⁴N^-1s^-1, and fluid fraction Φ_f=0.83. The testing data used in the frequency domain tests is taken from [1]: width w=3.0mm, height h=1.0mm, Poisson's ratio ν=0.126, permeability κ=1.16×10^-15m⁴N^-1s^-1, and fluid fraction Φ_f=0.8. The modulus of elasticity used is E=0.442855MPa and E=0.45MPa for two-dimensional and axisymmetric applications, respectively.

5.1 Frequency domain tests

The 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 two-dimensional and axisymmetric problems, respectively. They are implemented on the single-element mesh and the (regular, unbiased) 2×2- and 4×4-element 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 two-dimensional 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 p-refinement is typical of hybrid-Trefftz elements and, for coarse meshes, the rate of convergence under h-refinement is directly affected by the boundary conditions of the problem, which affect differently the alternative stress and displacement models. In the two-dimensional 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×4-element mesh) to 3.1 (single-element 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 two-dimensional and axisymmetric problems.

Figure 2: Convergence patterns under p- and h-refinement of the solution of plane strain problems.

Figure 3: Convergence patterns under p- and h-refinement of the solution of axisymmetric problems.

5.2 Time domain tests

The 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 two-dimensional 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 non-periodic spectral decomposition method described in [18]. No smoothing is used in the representation of the (double precision) finite element solutions presented below.

The displacement-driven, ramp-load program shown in Figure 4 is designed to reach a deformation of 5% at instant t₀=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 non-periodic wavelet basis defined on the (unit) interval [20,21], with dimension N=256 (family 3 and refinement 7). This high-order 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.

Figure 4: Loading program for the compression tests.

Figure 5: Evolution in time of the confined compression test solution of two-dimensional problems.

Figure 6: Evolution in time of the confined compression test solution of axisymmetric problems.

The confined compression tests on two-dimensional 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 two-dimensional 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 co-ordinates 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. Closure

The 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 hybrid-Trefftz 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 inter-element 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, well-suited to adaptive refinement and parallel processing.

This paper presents only very small part of the results for hybrid-Trefftz stress and displacement elements for transient analysis of incompressible saturated porous media. They are all present in [4-10].

As for two-dimensional 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 degrees-of-freedom. 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 h-refinement 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 high-order 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 width-to-height ratios. Consequent upon this weakness, refinement of the discretization in the z-direction may be needed to solve adequately particular applications.

References

	[1]	V. C. Mow, S. C. Kuei, W. M. Lai, C. G. Armstrong, Biphasic creep and stress relaxation of articular cartilage in compression, Theory and experiments, Journal of Biomechanical Engineering (1980), Vol. 102, No. 73.
	[2]	M. A. Biot, General theory of three-dimensional consolidation, Journal of Applied Physics (1941), Vol. 12, pp. 155-164.
	[3]	M. A. Biot, Mechanics of deformation and acoustic propagation in porous media, Journal of Applied Physics (1962); Vol. 33, pp. 1482-1498.
	[4]	J. A. T. Freitas, M. Toma, Hybrid-Trefftz stress and displacement elements for axisymmetric incompressible biphasic media, Part I, Formulation, To be submited for publication (2008).
	[5]	M. Toma, J. A. T. Freitas, Hybrid-Trefftz stress and displacement elements for axisymmetric incompressible biphasic media, Part II, Assessments, To be submited for publication (2008).
	[6]	J. A. T. Freitas, M. Toma, Hybrid-Trefftz displacement elements for incompressible biphasic media, Part I, Formulation, To be submited for publication (2008).
	[7]	M. Toma, J. A. T. Freitas, Hybrid-Trefftz displacement elements for incompressible biphasic media, Part II, Assessments, To be submited for publication (2008).
	[8]	J. A. T. Freitas, M. Toma, Hybrid-Trefftz stress elements for incompressible biphasic media, Part I, Formulation, To be submited for publication (2008).
	[9]	M. Toma, J. A. T. Freitas, Hybrid-Trefftz stress elements for incompressible biphasic media, Part II, Assessments, To be submited for publication (2008).
	[10]	M. Toma, Modelling of hydrated soft tissues using hybrid-Trefftz finite elements, Lisbon, Portugal [PhD thesis]: Technical University of Lisbon, Instituto Superior Tecnico (2007), 228 p. [PDF, 9.5 MB]
	[11]	R. L. Spilker, J-K. Suh. Formulation and evaluation of a finite element model for the biphasic model of hydrated soft tissues. Computers and Structures (1990), Vol. 35, No. 4, pp. 425.
	[12]	J-K. Suh, R. L. Spilker, M. H. Holmes. A penalty finite element analysis for nonlinear mechanics of biphasic hydrated soft tissue under large deformation. International journal for numerical methods in engineering, 32:1411, 1991.
	[13]	R. L. Spilker, T. A. Maxian. A mixed-penalty finite element formulation of the linear biphasic theory for soft tissues. International journal for numerical methods in engineering, 30:1063, 1990.
	[14]	M. E. Vermilyea, R. L. Spilker. A hybrid finite element formulation of the linear biphasic equations for hydrated soft tissue. International Journal for Numerical Methods in Engineering, 33:567, 1992.
	[15]	M. E. Vermilyea, R. L. Spilker. Hybrid and mixed-penalty finite elements for 3-D analysis of soft hydrated tissue. International journal for numerical methods in engineering, 36:4223, 1993.
	[16]	J. S. Wayne, S. L. Y. Woo, M. K. Kwan. Application of the u-p finite element method to the study of articular cartilage. Journal of Biomechanical Engineering, 113:397, 1991.
	[17]	C. W. J. Oomens, D. H. Van Campen, H. J. Grootenboer. A mixture approach to the mechanics of skin, Journal of Biomechanics, 20:877-885, 1987.
	[18]	J. A. T. Freitas, Mixed finite element formulation for the solution of parabolic problems, Comput Methods Appl Mech Eng (2002), Vol. 191, pp. 3425-3457.
	[19]	J. A. T. Freitas, Hybrid finite element formulations for elastodynamic analysis in the frequency domain, Int J Solids and Structures (1999), Vol. 36, pp. 1883-1923.
	[20]	I. Daubechies, Orthonormal bases of compactly supported wavelets, Communications Pure and Applied Mathematics 1988; 41:909-996.
	[21]	P. Monasse, V. Perrier, Orthonormal wavelet bases adapted for partial differential equations with boundary conditions, SIAM J Numerical Analysis, 1998; 29:1040-1065.