|Home | About | Journals | Submit | Contact Us | Français|
The finite element method is used in biomechanics to provide numerical solutions to simulations of structures having complex geometry and spatially differing material properties. Time-varying load deformation behaviors can result from solid viscoelasticity as well as viscous fluid flow through porous materials. Finite element poroelastic analysis of rapidly loaded slow-draining materials may be ill-conditioned, but this problem is not widely known in the biomechanics field. It appears as instabilities in the calculation of interstitial fluid pressures, especially near boundaries and between different materials. Accurate solutions can require impractical compromises between mesh size and time steps. This article investigates the constraints imposed by this problem on tissues representative of the intervertebral disc, subjected to moderate physiological rates of deformation. Two test cylindrical structures were found to require over 104 linear displacement-constant pressure elements to avoid serious oscillations in calculated fluid pressure. Fewer Taylor–Hood (quadratic displacement–linear pressure elements) were required, but with complementary increases in computational costs. The Vermeer–Verruijt criterion for 1D mesh size provided guidelines for 3D mesh sizes for given time steps. Pressure instabilities may impose limitations on the use of the finite element method for simulating fluid transport behaviors of biological soft tissues at moderately rapid physiological loading rates.
Biological connective tissues can be considered as materials having both solid and fluid phases. Therefore, they have time-dependent properties due in part to the viscous behavior of the fluid–solid interactions, as well as having solid elastic and viscoelastic properties. Since most biological tissues have anisotropic properties and complex geometry, the finite element method is often used to analyze the behavior of anatomical structures. They are commonly represented as poroelastic or biphasic materials.23
The field of biomechanics increasingly employs analyses that include both solid stresses and fluid flow associated with fluid pressure gradients. Such analyses attempt to link the mechanical stimuli and biological response of hydrated tissues under physiological loading. For example, Prendergast et al.27 proposed a model in which interstitial fluid flow, driven by pore pressure, was a major stimulus of mechano-regulation of tissue differentiation. This theory has a broad potential application to the understanding of, for example, fracture healing, bone remodeling, implant osseointegration, and cartilage degeneration. It has shown good correspondence to in vivo studies of tissue differentiation.15 The relationship between intervertebral disc compression, interstitial fluid motion, solute transport, and local metabolic response has been explored in the triphasic model of Huang and Gu.12
Poroelastic theory deals with materials and tissues having an elastic solid phase with fluid saturated pores. These materials can be considered to be governed by linear elastic theory and by Darcy's law for viscous flow through a porous medium. Values of the tissue properties can be estimated from experiments in which a sample of tissue is tested one-dimensionally in confined compression. The solid phase may also have viscoelastic and creep properties.13 The effective stress principle of saturated porous solids was developed for use in soil mechanics by Terzaghi and Peck,29 and subsequently generalized by Biot,3,4 and is known as Biot's Consolidation theory. The theory is expressed as a system of first-order differential equations that is integrated by the use of an appropriate time marching scheme.6 However, numerical complexities arise in the solutions because the basis functions in the u–π (solid-displacement and pressure) poroelastic formulation are of differing order. Stresses are related linearly to relative displacements (strains), while the fluid pressures associated with viscous flow are related linearly to flow velocity (first derivative of fluid displacements).1 Consequently, it has been reported by several authors5,8,10,22,31 that inaccuracies in pressure estimates can result from spurious spatial fluctuations of estimated pore pressure between adjacent elements. This is especially likely to occur at permeable boundaries in the initial time-marching steps, or when short time steps are used to represent rapidly changing or transient applied loads or pressures.1 These possible instabilities and their prevention in poroelastic analyses are not often referenced in the biomedical engineering field.
Vermeer and Verruijt30 give a physical explanation of this pressure instability problem. Small time steps and fast changes in external loads produce consolidation at the surfaces; consequently the derivative of the pore pressures at the draining surface or at the interface between strata of different permeability may be singular. Thus, stability considerations require a long-time step, but a contradictory short-time step is required to simulate the fast loading conditions. They derived an expression for the relationship between mesh size, time step, and material properties in the one-dimensional case. Also Ferranato et al.10 provided a similar empirical relation for a lower bound critical integration step length (Δtcrit) below which ill-conditioning may suddenly occur. The critical time step is larger for soft and low permeable porous media discretized on coarser grids and is similar to the criterion of Vermeer and Verruijt.30 However, the need for shorter time steps to model rapid loading and for smaller mesh dimensions to provide stability may impose unacceptable compromises, limited by computational costs.
Improved stability during fast loading simulations has been achieved by addressing the differing order of the solid strain and fluid pressure components, by use of elements employing differing degree piecewise polynomials for the displacement and pressure fields. Typically, elements employ quadratic interpolation of the solid strains and linear interpolation of the fluid pressure. Such elements (e.g., Taylor–Hood elements28) require additional nodes on their edges, hence increased degrees of freedom for their numerical solution. However, Murad et al.25 noted that, even with Taylor–Hood elements, potentially unstable results still occurred when there were discontinuous initial conditions at the initiation of a time-marching solution.
An expression was derived by Vermeer and Verruijt30 for the relationship between mesh size, time step, and material properties that was consistent with stable finite element solutions. They examined a one-dimensional u–π finite element model with linear interpolation elements and a uniform mesh of a linear porous material. They gave the required time step Δt as a function of the specific gravity of the fluid γ, the length of the elements (mesh dimension) Δh, the elastic modulus E, and the permeability of the porous material k:
Or, expressed in terms of the minimum mesh dimension Δh:
This criterion is cited in some commercial finite element software packages (e.g., Abaqus documentation manual, version 6.7.1, section 1.14.1).
Note that the solid–fluid friction coefficient fsw associated with fluid viscosity is related to the inverse permeability (k–1) by the porosity n (equal to the fluid/total volume ratio Vfluid/Vtotal):
This article provides empirical evidence of instability and associated loss of accuracy in poroelastic analysis of biological soft tissues at physiological rates of loading or deformation, with reference to tissues of the spine. In the structures of the musculoskeletal system, there are tissues with a wide range of properties. In the spine, vertebral bone and the nucleus of intervertebral discs have extremely different properties, with disc annulus possessing properties between these extremes.7,20,24 Articular cartilage and intervertebral discs are relatively soft and have low permeability.11,17 Discs are known to undergo substantial volumetric changes under physiological conditions as a result of internal fluid flow and fluid flow across tissue boundaries.21 Many relevant physiological loading rates are fast (e.g., gait at ~1 Hz) compared to the long creep and relaxation times of intervertebral discs, which are on the order of hours.15,19
Illustrative minimum time steps for the material properties of the intervertebral disc derived from Eq. (1) are given in Table 1, assuming that γ = 104 N/m3 (water) and for the case of mesh element characteristic dimension Δh = 1 mm. The values indicate that events having duration less than 17 min cannot be simulated with mesh sizes of 1 mm, since the minimum time step required for disc nucleus tissue would be over 1000 s, and for annulus over 500 s. The material properties are constrained to their actual values so the only parameter that can be altered to satisfy the Vermeer–Verruijt criterion is the characteristic element size (Δh), (see Fig. 1). The number of elements, and hence the number of degrees of freedom, increases inversely as the cube of the mesh size—hence reduction of element size will be limited by memory and computational time constraints. This limitation would have similar profound implications for modeling a wide variety of other connective tissues having similar properties (cartilage, blood vessels, etc.), and indicates that extremely fine element meshes would be required to obtain accurate solutions.
To determine the extent of the problem of stability and accuracy of poroelastic analysis of slow-draining tissues undergoing fast loading in a three-dimensional analysis, and to search for potential solutions in modeling biological materials, we examined two materials representing tissues of the intervertebral disc in a finite element analysis in an illustrative model having cylindrical shape. The analysis was run with differing boundary conditions and with two different element types—a linear displacement–constant pressure element and a quadratic displacement–linear pressure element.
The objective of this study was to determine empirically the circumstances under which pressure instability problems occur for tissue properties representing those of the intervertebral disc. The nucleus and circumferential direction annulus properties were employed, since they are the worst and best cases in Table 1. The limitations on mesh refinement and time step size reported in a one-dimensional analysis by Vermeer and Verruijt30 were evaluated for applicability to a three-dimensional analysis of intervertebral disc tissue.
A test case consisting of a cylinder with initial height and diameter of 5 mm was simulated in three-dimensional finite element analyses. Its material properties represented the nucleus and annulus of the intervertebral disc. The tissue properties were as given in Table 1, and for both tissues, the solid Poisson's ratio (ν) was set to 0.4. Two different boundary conditions at the loaded surfaces provided either free (slippery) or fixed (adhesive) displacements relative to the rigid non-porous or porous loading platens. In all cases, the fluid was allowed to flow in and out on the lateral surface.
In these illustrative simulations, a vertical compressive saw tooth (triangular) wave displacement with peak amplitude of 0.25 mm was applied to the top surface of the cylinder over 1 s, and then back to the original state in one-second. For each case, the applied force and the internal pressures were calculated for each of 20 steps (the time steps were 0.1 s).
The finite elements were hexahedral (brick) element with either 8 or 20 nodes. Four different mesh refinements were investigated, with both height and diameter each divided into 5, 10, 20, or 50 elements (Fig. 2). These meshes had characteristic mesh dimensions of 1, 0.5, 0.25, and 0.1 mm, and had 160, 960, 7680, and 122800 elements, respectively. According to the Vermeer–Verruijt criterion, the required number of elements to analyze this structure with material properties of disc nucleus, and for time steps compatible with applied loading frequency of 1 Hz would be in the order of 106 elements (Fig. 1).
Several cases (permutations of the two tissue material properties, linear, or quadratic (Taylor–Hood) displacement elements and the boundary conditions) were investigated, and results are presented for the following:
The outcome measures were:
The Abaqus finite element solution software (Abaqus/CAE Version 6.7-1, Dassault Systèmes Simulia Corp., Providence, RI) was employed. This analysis package uses the u–π method (Abaqus Users manual, Abaqus/CAE Version 6.7-1, www.simulia.com, 2008), in which the continuity equation is satisfied in the finite element model by using excess wetting liquid pressure as the nodal variable, interpolated over the elements. The equation is integrated in time by using the backward Euler approximation. The total derivative of this integrated variational statement of continuity with respect to the nodal variables is required for the Newton iterations used to solve the nonlinear, coupled, equilibrium, and continuity equations. Two pore pressure element types were compared. The C3D8RPH (8-node brick, trilinear displacement, hybrid with constant pore pressure, reduced integration) and the C3D20RPH [20-node brick, triquadratic displacement, hybrid with trilinear pore pressure, reduced integration (Taylor–Hood element)] elements were used in the simulations. Although an axisymmetric simplification would be possible for this cylindrical model, a full 3D simulation was made, to represent the case of analyses of structures having geometrically complex anatomy. The comparison for linear and quadratic coarse meshes were carried out by a personal computer (Intel® core™ 2 CPU, 2.14 GHz, with 2 GB RAM). For the very refined meshes, a computational server (DELL PowerEdge 6850, 4x—Intel® Xeon® Dual Core 7120M, 3.0 GHz, 4 MB L3 Cache, 800 MHz FSB, 16 GB RAM) was used.
The outputs of the simulations were highly dependent on the mesh refinement, except in the case of the reaction force at the upper surface of the cylinder, which was consistently calculated in all cases (max deviation 3.7% between finest and coarsest mesh). In contrast, the calculated fluid pressures were dependent on mesh size. In the case of the pressure calculated at the center of the lower cylinder end (Fig. 3), the values were substantially over-estimated when using the linear displacement–constant pressure elements for coarse-mesh models, especially in adhesive platen boundary conditions with non-porous platens. The pressure values converged with use of the Taylor–Hood elements (lowest two panels in Fig. 3). However, in nearly all simulations the calculated pressure distribution along a radial line at the cylinder's mid height (results for non-porous platens in Fig. 4) illustrated pressure instability, with oscillating values between adjacent nodes. With use of the finest mesh (0.1 mm) for nucleus and the 0.25 mm mesh for annulus the calculated pressure did not have large variations between adjacent elements, except at the outer edge, which was a draining boundary.
With use of the 20-node (Taylor–Hood) elements, pressure oscillation was evident for 1 and 0.5 mm mesh sizes. This comparison is presented only for material representing annulus (Figs. 3, ,4),4), since it was a ‘best case’ (less liable to instabilities). Thus, quadratic elements did not resolve the instability substantially better than simple mesh refinement in the less ‘demanding’ case of annulus. The computational cost of using these elements precluded simulations with any finer meshes because of the amount of addressable RAM in the computer available for these simulations (PC, Intel® Core™ 2 CPU, 2.14 GHz, 2 GB RAM). The computational cost of using the Taylor–Hood elements (compared to linear elements) was six times greater with a 1 mm element size, and 11 times greater with 0.5 mm elements, when measured as the total computational time (problem assembly and solver time). The simulations with tissue properties representing annulus converged faster with mesh refinement as expected, since the elastic modulus and permeability are both greater.
These simulations demonstrated that the use of the finite element method to obtain solutions for poroelastic 3D representations of intervertebral disc tissue loaded at a moderately rapid physiological rate with practical levels of mesh refinement is inaccurate, with spatially oscillating values of the pore pressures. Oscillation of the calculated fluid pressure, especially at the lateral borders of the cylinder, was evident for all test cases presented in Fig. 4. While the conditions tested were for tissue having material properties of intervertebral disc tissues, the same solution problems would be expected in other biological soft tissues and other anatomical structures. Pressure fluctuations are likely to occur both at free draining external boundaries, and at discontinuities between tissue boundaries (e.g., at bone-soft tissue boundaries, or between nucleus and annulus of the intervertebral disc). These boundaries are typically regions of higher cellularity of intervertebral disc and cartilage.
Spurious pressure fluctuations were reduced but not eliminated with use of the C3D20RPH element, but were also accompanied by a large increase in computational cost which was due to the large increase in degrees of freedom in the model (20-node elements vs. 8-node elements).
The calculation of the total force in each case was consistent with differing mesh sizes, indicating that the fluid pressure oscillations gave correct total values when averaged over the loading surfaces. The pressure oscillations in the solutions were evident (especially at the outer draining boundary edges—see Fig. 4) except with mesh refinement approaching that indicated by the Vermeer–Verruijt criterion corresponding to 122800 elements (nucleus) and 7680 elements (annulus). Therefore, in dynamic loadings of complex anatomical structures at higher frequency (e.g., full intervertebral disc) the number of elements required for reliable poroelastic simulations could become prohibitively large, since the memory required for the large number of degrees of freedom in the model may exceed the addressable RAM of the computer used. The calculated pressures were also sensitive to boundary conditions. With the adhesive platens as compared to the slippery platens, larger oscillations were observed with coarser meshes (Fig. 3).
By use of Eq. (2), the minimum required mesh size can be estimated based on the Vermeer–Verruijt criterion, and this was 10 microns for nucleus and 49 microns for annulus, when using linear elements. These values are somewhat smaller than those observed empirically, as the test-case simulations were observed to be close to convergent with mesh refinement with 100- and 250-micron mesh dimensions, respectively.
Various remedial actions have been proposed to prevent the pressure instabilities. Biased meshes provide an efficient version of mesh refinement that can be employed near the drained boundaries. However, implementation of such a biased mesh might be challenging for more complex simulations (e.g., intervertebral disc), which contain multiple boundaries between tissue regions with large differences in material properties, with each boundary requiring mesh refinement. Vermeer and Verruijt30 showed a small improvement by employing quadratic interpolation (rather than linear interpolation) between nodes in the 1-dimensional case, and this changed the denominator in Eq. (1) from 6 to 10. The Taylor–Hood elements28 require additional nodes on their edges, hence increased number of degrees of freedom for their numerical solution. However, in the simulations reported here, pressure instabilities were still evident with use of these elements, as also reported by Murad et al.25 when there were discontinuous initial conditions in a time-marching solution.
Aguilar et al.1 proposed a numerical stabilization by adding an extra term in the flow equations that permitted piecewise linear approximation for displacements and pressure, and this parameter appeared to be optimal at least in the one-dimensional case. Another stabilizing approach was proposed by Vermeer and Verruijt30 who suggested a “double calculation” procedure as a method to circumvent the minimum time step bound. In the first calculation, the drained boundary is assumed to be impermeable (with the fluid “somewhat” compressible). In the second calculation, the boundary pressure is set to zero. In exploratory analyses using the Abaqus finite element package, applying no permeability in the first step, and allowing the fluid drain in the next step was found to achieve improved convergence.
This study confirms that 3D poroelastic finite element analysis of biological tissues having low values of elastic modulus and permeability, loaded at a moderately rapid physiological rate, can be limited by instability in calculations of pressure. The 1D analysis published by Vermeer and Verruijt30 provided a guide to the required mesh size within an order of magnitude in the 3D examples tested here. This indicates that finite element analyses of ‘slow draining’ biological tissues at fast loading rates must be designed with great care. Extremely fine and/or locally refined meshes near external boundaries may be required in complex structures having many interfaces such as the intervertebral disc. This can be very computationally intensive. Calculated pressures should be checked carefully to ensure that numerical instabilities are not present in the solution. In a simulation having numerical instability, the magnitudes of oscillations are inherently unpredictable, producing unknown magnitude and type of errors.
This work was supported by NIH Grant AR R01 049370 and the Swiss National Science Foundation (NCCR CO-ME). Authors acknowledge helpful discussions and input from Gerard Ateshian, Weiyong Gu, Raghu Natarajan, Robert Spilker, Hai Yao, and Arthur Michalek.