J. Biomechanrcs Vol. 25, No. 9, pp. 1027-1045. Printed m Great Britain

MI21 9290/92 $5.00 + .oO Pergamon Press Ltd

1992.

A TRANSVERSELY ISOTROPIC BIPHASIC FINITE ELEMENT MODEL OF THE MENISCUS ROBERTL. SPILKER,* PETER S. DONZELLI* and VAN C. Mow? *Department of Mechanical Engineering, Aeronautical Engineering, and Mechanics, Rensselaer Polytechnic Institute, Troy, NY 12180-3590,U.S.A. and TDepartments of Orthopaedic Surgery and Mechanical Engineering, Columbia University, New York, NY 10032,U.S.A. Abstract-A finite element model of the meniscus is presented, based on an axisymmetric geometric approximation of the menisci and a biphasic description of the tissue as a mixture of solid and fluid components. The highly fibrous nature of the meniscal tissue is accounted for by using a fiber-reinforced, transversely isotropic description of the solid phase. This model is used to study the response of a meniscus resting on a perfectly lubricated tibia1 surface and subjected to distributed loads applied to the femoral surface, and to examine the effects of changes in loading conditions at the femoral and tibia1 interfaces. Quantities of interest include the stress, pressure and strain distributions at discrete times early in the meniscal response, and the flow of the fluid phase relative to the solid phase. Of particular interest are regions of large tensile strains which could lead to meniscal failures such as the bucket-handle tear. We show that all components of strain are positive in regions of the outer third of the meniscus, and that the maximum tensile strain perpendicular to the circumferentially arranged fibers (largest principal strain in the axisymmetric cross section) is positive throughout most of the cross section. Changing the partition of the load on the femoral surface and the permeability at the tibia1 surface changes the time-dependent response, but has little effect on the strain distributions at times of the order of 5 s considered in this study. The inclusion of a transversely isotropic, fibrous representation of the solid phases is shown to be essential to proper meniscal simulation. The results demonstrate the importance of the biphasic representation since the fluid phase is shown to carry a significant part of the applied load.

INTRODUCTION The menisci in the human knee are crescent-shaped structures interposed between the relatively flat tibia1 surface and the convex femoral condyles (Fig. 1). The cross section of the meniscus is nearly triangular, with a relatively straight distal (tibial) surface and concave proximal (femoral) surface, to conform to their adjacent surfaces. The medial and lateral menisci differ somewhat in their size with the more circular-shaped lateral meniscus covering a larger portion of the tibia1 articular surface. The primary attachments of the menisci are to the tibia at the horns (Arnoczky et al., 1988). Meniscal tissue is hydrated and highly fibrous, with dominant fiber bundles oriented in the circumferential direction. Some smaller and less dense radial fibers running perpendicular to these large bundles are also observed. This fibrous architecture produces a tissue which is much stiffer in the fiber direction than in directions perpendicular to the fibers. As described initially by Shrive et al. (1978), this fiber arrangement follows the function of the tissue. Assuming that the interactions between the meniscal and adjacent artitular cartilage tissues are frictionless, then loads transmitted to the meniscus must be perpendicular to the meniscal surface. On the femoral surface in particular, the load will have both an axial and a radial component. The radial component of force, which

Received injnalfirm

13 November 1991.

tends to expand the crescent-shaped meniscus, is resisted primarily by the internal forces in the circumferential direction, and, thus, the stiff circumferential fibers provide the required resistance. The inclusion of this unique fibrous structure is, therefore, essential in the simulation of a meniscal mechanical response to loading. Like articular cartilage, meniscal tissue exhibits a time-dependent, viscoelastic behavior under applied load. To represent this behavior we use the biphasic model of soft hydrated tissues, which is based on the theory of mixtures, and was first developed for soft tissue as the linear KLM biphasic model by Mow et al. (1980). In this theory, the tissue is modeled as a two-phase mixture of fluid (interstitial fluid) and solid (collagen and proteoglycan) forming a fluidfilled, porous permeable domain. The resistance of the fluid as it flows through the porous permeable solid phase gives rise to the time-dependent response of the tissue. While each of the constituents is assumed to be intrinsically incompressible, the tissue as a whole is compressible due to fluid exudation. Nonlinearities which are important at physiologically relevant levels of loading can also be included, corresponding to (i) the change in local permeability and/or fluid-solid content with changes in the compressive strain (Lai and Mow, 1980; Mow et al., 1984; Holmes et al., 1985), and (ii) the intrinsic nonlinear behavior of the solid phase of the tissue (Mak, 1980; Holmes, 1986; Mow et al., 1986; Kwan et al., 1990). This biphasic continuum description of soft hydrated tissues has been applied to a number of tissues, including articular cartilage

1027

1028

R. L. SPILKER

et al.

Quadriceps Tendon

apatellar Pad of Fat

Cancellous Bone

Tibia

(a) Human Knee Joint

Cancellous Bone

ubchondral Cortex

(b) Load Bearing Region Fig. 1. Schematic of the menisci in the knee joint.

and meniscus (Armstrong and Mow, 1982; Armstrong et al., 1984; Mow et al., 1984, 1989; Whipple et al., 1985; Brown and Singerman, 1986; Mak et al., 1987; Salzstein and Pollack, 1987; Salzstein et al., 1987; Chern et al., 1989, 1990; Fithian et al., 1989, 1990; Proctor et al., 1989). Recently, Lai et al. (1989, 1991)

have developed a triphasic theory to include the dissolved solute concentration which is known to be responsible for the Donnan osmotic pressure and the swelling effect of soft tissue. A comprehensive solution of the meniscus problem under compression will provide details of the deformation, flow fields, and internal stresses, pressure and strains as a function of time. The finite element method is a natural choice for numerical solution of partial differential equations on such complex geometries. A number of finite element models of the meniscus have been reported. Initially, axisymmetric toroidal models of the meniscus with frictionless contact between the meniscus and deformable femur and

tibia were developed in which the materials were represented as single-phase linear elastic materials (Sauren et al., 1984; Hefzy et al., 1987; Hefzy and Zoghi, 1988) or transversely isotropic materials (Schreppers et al., 1990). Nonlinearities were introduced into an axisymmetric model resting on a rigid frictionless tibia1 surface with loads applied to the femoral surface (Aspden, 1985; Aspden et al., 1986); tensile properties assumed in the circumferential direction were nonlinear, and compressive properties in the cross section were chosen to match the assumed elastic modulus of cartilage in this model. Axisymmetric contact analyses with transversely isotropic nonlinear single-phase materials have been presented more recently (Tissakht et al., 1991) and demonstrate the importance of the fiber-reinforced description. Three-dimensional models of the meniscus including contact and single-phase material representations allowing for different properties in the circumferential direction have also been developed (Tissakht et al.,

A model of the meniscus 1989; Tissakht and Ahmed, 1990). Important insights into meniscal behavior have been obtained with these models, However, each of these analyses assumes a single-phase solid meniscal material and, therefore, cannot capture the important fluid flow characteristics within the tissue which provide the time-dependent behaviors. A finite element method which solves the equations of the biphasic theory has been developed by Spilker and co-workers (Spilker and Maxian, 1990; Spilker and Suh, 1990; Vermilyea and Spilker, 1992). The extended formulations also include the nonlinearities due to strain-dependent permeability (Spilker et al., 1988) and finite deformation (Suh, 1989; Suh et al., 1991). In this study we will use a version of the finite element codes which is based on the linear KLM biphasic model for soft hydrated tissues and incorporates a transversely isotropic fibrous representation of the solid phase of meniscal tissue (Spilker and Maxian, 1990). We will maintain loading in our simulations at levels for which the linear theory can be reasonably applied. A correlation of the theoretical and experimental studies on the mechanical testing of soft hydrated tissues such as articular cartilage suggests that the linear theory is adequate for local compressive strains as high as 20-30% (Mow et al., 1980; Kwan et al., 1990). The experimentally measured uniaxial tension stress-strain response of bovine (Whipple et al., 1985; Proctor et al., 1989) and human (Fithian et al., 1989) meniscal tissue oriented in the circumferential direction is nonlinear. However, reasonable linear fits can be assumed in the range of 2-8% strain. The measured uniaxial stress-strain response in the radial direction (Whipple et al., 1985; Proctor et al., 1989) (perpendicular to the predominant fiber direction) is reasonably approximated as linear up to engineering strains as high as 50%. Strain levels within these magnitudes will be maintained with physiologically relevant load levels approximately equivalent to the body weight. These linear models will be used to calculate the deformation and fluid flow, pressure, stress and strain fields within the meniscal tissue under applied loading, and to identify possible pathologies of meniscal failure. Attention in this study is focused primarily on the tensile strains which arise in the tissue. It is likely that combinations of tensile strains are responsible for pathologies such as the bucket-handle tear. In future studies, nonlinearities can be added in order to extend the simulations to higher levels of loading.

1029

negligible. In the equations which follow, boldface denotes a tensor, superscripts s and f denote solid and fluid phases, respectively, and a bar over a quantity denotes a prescribed value of that quantity. The symbol V is the gradient operator, Vsym refers to the symmetric part of the gradient, C$represents a constituent volume fraction (note that & + #’ = 1 for the fully saturated mixture assumed in the biphasic theory), v velocity, u stress, II the body force caused by diffusive drag, p the pore pressure, I the identity tensor, D the fourth-order tensor of material constants of the solid phase, E strain, u displacement, t traction, and K a diffusive drag coefficient which, for slow flow, is inversely related to the permeability K via the relation (Lai and Mow, 1980) E(

FINITE ELEMENT

AND

MODEL

The linear biphasic equations for soft tissue (Mow that: the biphasic mixture is chemically inert and immiscible; each phase is intrinsically incompressible; the solid phase is linearly elastic; the strains are infinitesimal; the fluid phase is inviscid; and the inertia forces are et al., 1980) are based on fhe assumptions

(“)’ K

The governing equations and boundary conditions for domain R with boundary I consist of a continuity equation for the mixture, v *(@vS+ #v’) = 0, linear momentum phases,

equations

(24

for the solid and fluid

V*a”+l-IJ=O,

(2b)

v~a’+I-I’=o,

(2c)

constitutive relations for each phase, tr’ = - @pI + DE’,

(24

ur = - f#JfpI,

(2e)

and the body force corresponding ance

to diffusive resist-

I-I’ = - I-I’ = K(v’ - v’) - ~VI$‘.

(2f)

It is assumed that the solid strains are related to the solid-phase displacements through the infinitesimal strain-displacement relation &”= (VU’)BYrn.

(2g)

Boundary conditions correspond to tractions applied to the solid and fluid phases: tS= f”

on I,s,

(3a)

p = fi on Ilr,

(3b)

and to prescribed values of solid displacement or fluid velocities as a function of position and time: us = I’

DESCRIF’TION OF THE MATHEMATICAL

_

on I..,

vf = ii’ on

rvr.

(3c)

(34

Initial conditions complete the problem description. Note that for isotropic materials, D contains only two independent constants. For the transversely isotropic description used in this study, five independent material parameters and the direction of the fibers are required to define D, details are presented later. Furthermore, in the present study we assume that 4’ is

R. L. SPILKER

1030

constant throughout the domain and, thus, the last term in equation (2f) is eliminated. A number of finite element formulations have been developed for these biphasic equations, including those based on penalty methods (Spilker and Suh, 1990), mixed-penalty methods (Spilker and Maxian, 1990), and hybrid methods (Vermilyea and Spilker, 1992). These and other formulations have also been extended to nonliner biphasic equations (Oomens et al., 1987; Spilker et al., 1988; Suh, 1989, Suh ef al., 1989, 1991). For the present study, we consider the linear biphasic model, and use the mixed-penalty formulation (Spilker and Maxian, 1990). In the mixedpenalty formulation, the continuity equation [equation (2a)] is replaced by the following penalty form:

v -(fp

+

q5’v’)+ ; = 0,

where /I is the user-specified penalty number. For large fl, equation (4) reduces to the continuity equation, and p as calculated from equation (4) is the pressure. Guidelines for the selection of/l (Spilker and Suh, 1990) indicate that the order of magnitude of fl should be

where 4 = @ or of, 1’ and $ are the Lame constants of the solid phase for an isotropic material (we will use appropriate substitutions for the transversely isotropic material), to is a characteristic time scale (e.g. ramp time), and c is set to the order of 106-1010. Experience with penalty methods suggests that variations by the four orders of magnitude implied by the range of c have negligible effect on the computed solution. The finite element formulation utilizes the Galerkin weighted residual method for the solid and fluid phases. The momentum equations are expressed in terms of the solid- and fluid-phase kinematic variables by substituting equations (2d)-(2g) into equations (2b) and (2~). The resulting momentum equations, the penalty form [equation (4)] of the continuity equation [equation (2a)] and the natural boundary conditions (traction conditions) [equations (3a) and (3b)] are substituted into a weighted residual statement. Equations (3~) and (3d) are treated as the essential boundary conditions and are, therefore, exactly satisfied by all approximations. The resulting Galerkin finite element formulation is expressed in terms of unknown variables corresponding to the displacement of the solid phase, velocity of the fluid phase and pressure. Standard families of continuous interpolation functions are applied to the kinematic variables and simple polynomials (discontinuous between elements) are used to interpolate the pressure. In the resulting ‘biphasic’ element the pressure variables are eliminated at the element level. Using conventional

et al.

assembly and constraint operations a system of firstorder differential equations is obtained in the form: Cv + Kd = F,

(6)

where C, K, and F are the assembled conductivity and stiffness matrices and force vector, respectively, and d and v are the nodal displacement and velocity degrees of freedom, respectively, at the nodes for the complete (assembled) continuum for the solid and fluid phases. These equations are solved using the finite difference method for the response at discrete time steps, subject to a set of initial and boundary conditions. Once the nodal displacements and velocities of the phases are known at a time step, the intra-element stress and strain components and the pressure can be calculated, producing a complete solution at that time step. We use a global leastsquares smoothing technique (Niu and Shephard, 1990) to produce the stress, strain, and pressure contours presented in this study. In this approach, nodal values of the quantity of interest (e.g. a stress component, strain component, or pressure) are calculated as the least-squares best fit to the values of the quantity at the optimal sampling points within the element. These smoothed nodal values are then used directly in polygonally based contouring algorithms in he MOVIE.BYU post-processor (Engineering Computer Graphics Laboratory, Brigham Young University, Provo, Utah). DESCRIPTION

OF THE MENISCAL GEOMETRY,

BOUNDARY

CONDITIONS

AND LOAD

In order to obtain a precise description of the meniscal geometry and subsequent response, a threedimensional model requiring significant computational resources would undoubtedly be needed. We expect to eventually pursue such three-dimensional studies; however, as a first step in understanding meniscal response, and in order to perform more exhaustive parametric studies, a geometrically simplified representation will be employed. The actual meniscal geometry will be approximated by an axisymmetric toroidal geometry in which the cross section takes on a realistic wedge shape as indicated in Fig. 2. The problem is best described in cylindrical coordinates corresponding to the radial direction, r, the axial direction, z, and the circumferential direction, 8. In this axisymmetric description, the geometry, material properties, boundary conditions and loading do not change with 8. As a result of the axisymmetric assumption, the deformation, flow fields, stress and strain will be calculated in the cross-sectional r-z plane. Thus, the three-dimensional problem is reduced to a two-dimensional analysis which can be solved in the plane of the cross section. In the present description of the menisci, we assume that the distal surface is horizontal to conform to an assumed horizontal tibia1 surface. The proximal surface is concave to conform to the condyles of the distal

A model

of the meniscus

Z

,.-I ;i .

,

:

* :

I’

\I

(a)

‘.

m

t ___----__-- -____

Circumferential

Fibers

I

-..>_E:__-----.----____

-_

.---------mm__

___---

__-’ *-__ ----____________-----

I

r

Tibia1 Surface

(b) Fig. 2. Schematic of the (a) axisymmetric representation and (b) definition of geometry for the meniscal simulations.

femur. The peripheral surface of the meniscus is assumed to be convex. The meniscal geometry (Fig. 2) is, therefore, represented by the inner radius, ri, and outer radius, ro, of the distal (tibial) surface; the height, h, at the outer radius; the height, h,, at the mid-radius of the meniscus, expressed as a fraction of the maximum height, h; and the bulge, b, at the peripheral surface, expressed as a fraction of the radial dimension of the meniscus. The proximal and peripheral surfaces are formed by passing a second-order curve through the endpoint and midpoint coordinates. In our simulations, we examine the response of a meniscus which is resting on the tibia and with prescribed loads imposed by the femur. We assume that the interface between the meniscus and the tibia is perfectly lubricated (i.e. frictionless) and impermeable. Thus, the meniscus is free to expand radially at the tibia1 surface, but the displacement and fluid flow in the axial direction at the tibia1 surface are required to be zero. We will also consider a case when the fluid phase is free to flow out normal to the tibia1 surface, and show that this has a minimal effect on the early time response in which we are most interested. The traction along the tibia1 surface is always compressive for the applied loads used in this study and, thus, the meniscus would not separate from the tibia1 surface even if the model allowed for this possibility. The effects of the femur are simulated by applying a load distribution to the femoral surface of the meniscus (Fig. 3). It is reasonable to assume that the interface between the meniscus and the femoral surface is

1031

also perfectly lubricated. The load will, therefore, be assumed to act normal to the femoral surface, since the tangential component, corresponding to frictional forces, must be zero. The total applied load must also be partitioned between the solid and fluid phases. In theory, two extreme cases can be identified. If the load from the femur onto the meniscus is passed directly to the solid phase, in which case the fluid is free to flow through the loaded surface and the pressure on the fluid is zero, then the load would be partitioned 100% to the solid phase and no load on the fluid phase. In the other extreme, if the load is imparted to the meniscus through the layer of fluid, then the theoretical load partition must be in accordance with the solid and fluid contents (Hou et al., 1989). For example, if the fluid content is 80%, then 80% of the load is applied to the fluid phase, and the remaining 20% is applied to the solid phase. In the present simulations, we will examine a range of partitions for the normal distribution of force applied to the femoral surface, but will concentrate on the case when the load is split 75% to the solid phase and 25% to the fluid phase to simulate a ‘mixed’ loading condition. For the present biphasic analysis, a triangular, axisymmetric element with curved sides, and nodes at the vertices and mid-side locations (i.e. six nodes per element) is used, as depicted in Fig. 4(a). The element is, therefore, a ring whose cross section in the r-z plane is triangular with curved sides. To construct a finite element model of this meniscal model, the cross section is divided in a regular fashion into a series of these triangular elements. Finite element meshes are generated using Finite Quadtree, developed by Shephard and co-workers (Shephard et a[., 1988); a sample mesh ofelements is shown in Fig. 4(b). Details of the formulation and implementation of this six-node element can be found in Spilker and Maxian (1990).

DEFINITION

OF TRANSVERSELY

SOLID-PHASE

MATERIAL

ISOTROPIC

PROPERTIES

The solid phase of the meniscal tissue is highly fibrous. With the exception of a thin superficial layer, the fibers are oriented primarily in the circumferential direction. As a result, it should be expected that the solid phase of the tissue will have very different properties (stiffness, for example) in directions parallel to the fibers in comparison to directions perpendicular to the fibers; this has in fact been verified in a number of experimental studies (Bullough et al., 1970; Chern et al., 1989, 1990, Proctor et al., 1989; Skaggs and Mow, 1990). Based on these known directional variations, we assume the solid phase to be transversely isotropic. In a transversely isotropic material, the mechanical properties parallel and perpendicular to the fiber direction are different, but all directions perpendicular to the fiber direction have the same properties.

1032

R. L.

SPILKER et al.

a. Distribution of Normal Parabolic Load Over the Perfectly Lubricated Femoral Surface .

b. Parabolic Pressure Distribution Applied via Fluid Fib at the Tibial Plateau

Fig. 3. Load distributions applied to the meniscus. (a) The load is applied normal to the femoral surface and partitioned between the solid and fluid phases. The total vertical component of the applied load is equal to 667 N (approximately body weight). (b) The tibia1 surface is made permeable, and a parabolic pressure distribution is applied to the fluid phase.

,_--

_____-*---

-----------.

. Nodal Points *.

= Film Dire&ion

(

I ‘.. :..:--- __-- -__-_______ --m____ ___--i,.6Node ,----------‘.._

--____

----______

___---

_____------

*.’

Triangular Element

t-----I

= Perpendicular to Fiber I a

=FiberAngle

(4 Fig. 5. Material coordinates used for the definition of material constants in the fibrous, transversely isotropic solid phase of the meniscal tissue.

(b) Fig. 4. {a) Geometry of an axisymmetric six-noded triangular ring element. (b) A typical finite element mesh for a cross section of the meniscus.

We define the material properties in a material coordinate system (xi, x2, x3) in which xi is in the fiber direction, and x1 and x3 are mutually perpendicular and perpendicular to the fiber direction, as depicted for an annulus of fiber-reinforced material in Fig. 5. With this choice of coordinates the stressstrain relation for the material can be defined with respect to stress and strain components in material coordinates. If the stress-strain relation is required in a different coordinate system, typically the global coordinate system, it can be obtained by standard tensor transformation laws. Subscripts 1, 2, or 3 are used to indicate the material coordinate directions associated with a particular material property.

Because the fibers in the meniscus are predominantly arranged in the circumferential direction (i.e the fiber angle in Fig. 5 is zero), therefore the material coordinates (xi, x2, x3) and global coordinates (0, z, r) are coincident, and the material constants could also be identified with respect to the global cylindrical coordinates. In summary, seven parameters are required (we use the material coordinates in this discussion): the Young’s moduli, El and &; the Poisson’s ratios, vi2 and v2s; the shear moduli, Giz and Gz3; and the angle(s), CC, describing the direction of the fibers with respect to a global reference coordinate system. The modulus El parallel to the fibers represents the ratio of the stress in the xi direction to the strain in the xi direction under the action of a load in the x1 direction. Similarly, the modulus Et perpendicular to the fibers represents the ratio of the stress in the x2 direction to the strain in the x2 direction under the action of a load in the x2 direction. The Poisson’s ratio vi2 is calculated as the ratio of the contractile strain in the x2 direction to the tensile strain in the x1 direction under the action of a load in the xi direction.

A model of the meniscus The Poisson’s ratio vz3 represents the ratio of the contractile strain in the xg direction to the tensile strain in the xz direction under the action of a load in the x2 direction. The shear moduli, G12 and Gz3, correspond to the ratio of shear stress to shear strain in a plane including the fibers and perpendicular to the fibers, respectively. The modulus Gz3 corresponds to the shear modulus of the material surrounding the fibers, for a transversely isotropic material, and is usually assumed to be related to Ez and vz3 by the relation G23 =

E2

2(I +

v23)’

These material properties have been measured experimentally for normal bovine and human meniscus (Whipple et al., 1985; Chern et al., 1989; Proctor et al., 1989; Fithian et al., 1990; Skaggs and Mow, 1990); these values form the basis for the range of material properties considered in this study. We will use El (Fithian et al., 1989; Proctor et al., 1989) and v12 (Fithian et al., 1989) values obtained from uniaxial tension tests of circumferentially oriented meniscal tissue specimens. The shear modulus Gi2 is determined from shear tests of meniscal cylindrical disk specimens having fiber orientation along the longitudinal direction of the cylinder (Chern et al., 1989, 1990). The value of E2 in tension can be obtained from uniaxial tension tests of radially directed specimens (Proctor et al., 1989; Skaggs and Mow, 1990). The transverse modulus, E2, used in this study is obtained from the unconfined compression response of cylindrical disks of meniscal tissue prior to cyclic shear testing (Chern et al., 1989, 1990). The value of v23 could be obtained from either of these tests. Details of the material property matrix for a transversely isotropic material, and the implementation of this material description into the biphasic element are found in the Appendix.

RESULTS

Geometric, material, and loading parameters

This study focuses on the time-dependent response of the meniscus as a result of loads applied to the femoral surface. To describe this response, we consider the distributions of stress, strain and pressure at discrete times and the time-dependent deformation, as characterized by the radial displacement at the outer radius along the tibia1 surface (point A in Fig. 2) and the axial and radial displacements at the midpoint of the femoral surface (point B in Fig. 2). The components of strain are of particular interest. It is reasonable to assume that large tensile normal strains can lead to failure of the meniscus. For example, large tensile axial strain or radial strain could lead to the commonly observed bucket-handle tear (Arnoczky et al., 1988) in the meniscus. The presence of large shear strains in conjunction with tensile normal strains

1033

could also contribute to tissue failure. We will, therefore, calculate the principal strains and, thereby, determine the maximum tensile strain in the cross section and in the circumferential (fiber) direction. Following an extensive examination of what we will term the ‘normal’ meniscus, we will consider some variations in the modeling of the loading conditions at the tibia1 and femoral surfaces and determine their effects on meniscal response. The geometry assigned to the normal meniscus is chosen on the basis of an examination of human meniscal specimens. The material properties are chosen as described earlier to be within the range measured experimentally (Whipple et al., 1985; Chern et al., 1989, 1990; Fithian et al., 1989; Proctor et al., 1989; Skaggs and Mow, 1990). These geometric and material parameters are summarized in Table 1. It should be noted that the effects of variations in the geometric and material parameters on meniscal response are of interest and will be the subject of a future publication. A distributed load which is parabolically shaped in the plane of the cross section and uniform in the circumferential direction is applied to the toroidal meniscus model, to represent the action of the femur on the meniscus [Fig. 3(a)]. This canonical distribution, used in the absence of experimental data or analytic studies suggesting a more precise load distribution, will be varied in future studies. The load is applied instantaneously at time t = 0 s (i.e. a step load) and held at a constant value; the initial load partition is also held constant. As a result, a timedependent creep response of the meniscus will be found in which the displacement and strain increase at a given position as a function of time, eventually reaching an equilibrium value. We will focus most of our attention on the response at early times (e.g. t = 5 s).

At each point on the femoral surface, the applied load can be decomposed into an axial and a radial component, and these components summed up over the complete femoral surface to determine the net axial force, F,, and net radial force, F,, applied to the meniscus. Both F, and F, will depend on the geometry of the femoral surface and the magnitude of the distributed load. For all cases in the present studies, we adjust the magnitude of the distributed load so that the net axial force is equal to the body weight (i.e. 675 N). This load level is at least twice the resting load on the meniscus, and could represent on the order of four times the resting value, depending on the share of total load carried by the meniscus. Previous singlephase finite element contact analyses of the meniscus suggest that the meniscus carries approximately 50% of the applied load; the actual load share depends on the mechanical properties of the meniscal tissue and the adjacent tissue (Sauren et al., 1984; Hefzy and Zoghi, 1988; Schreppers et al., 1990). The resulting peak value of the parabolic distribution of load is on the order of 0.3 MPa, which is small compared to peak values measured in uitro (Ahmed and Burke,

1034

R. L.

SPILKER

et al.

Table 1. Geometric and material parameters for normal meniscal tissue modeled as a biphasic material with a transversely’ isotropic fiber-reinforced solid phase Geometric parameters Inner radius Outer radius Height Peripheral bulge Mid-radius height Material parameters Solid and fluid content Permeability Young’s modulus Poisson’s ratio Shear modulus

ri = 25 mm r,=4Omm h=6mm b = & (r, - ri) = 1.5 mm h, = $ h = 2.25 mm & = 0.25, 4’ = 0.15 K = 1.26 x lo-l5 m4 N-’ S-I (Proctor et al., 1989) E, = EB = 200 MPa (Fithian et al., 1989; Proctor et al., 1989) E, = E, = E, = 0.055 MPa (Chem et al., 1989, 1990; Proctor et al., 1989; Skaggs and Mow, 1990) vi2 = 2.0 (Chern et al., 1989, 1990) v*J = 0.05 G12 = 0.025 MPa (Chem et al., 1989, 1990) G,, from equation (7)

1983; Ahmed, 1991). As will be demonstrated, the strain values predicted in this case approach the limit of the linear theory, i.e. 30% strain. Larger loads can be analyzed only with the nonlinear finite deformation biphasic theories, and are not considered in the present study. Results for the normal meniscus

The radial displacement at the outer radius of the tibia1 surface (point A in Fig. 2), and the axial and radial components of displacement at the middle of the femoral surface (point B), as a function of time, are shown in Fig. 6. The displacements jump to an initial value slightly after the imposition of the load (i.e. at t = O+) and continue to increase (creep) with time. At much larger times, the displacement will reach an equilibrium value. A similarly increasing, time-dependent response would be observed for all quantities of interest. While the response for long times may be of interest in some applications, we will focus much of our attention on the short-time response, since these times are of greatest interest clinically in the response to suddenly applied loads, and consider the distributions of stress, strain, pressure, and fluid velocity over the cross section at discrete times. Figures 7(a)-(e) show the distribution of stress components and pressure at time t = 5 s. The inset in each figure serves as a reminder of the direction of the component of stress being considered. In these figures, lines of constant stress or pressure are indicated and correspond to the values given in the legend. The contour bands represent regions with values in a certain range; hence, the maximum and minimum present in the cross section will, in general, exceed the largest and smallest values shown. The axial stress

[Fig. 7(a)] and radial stress [Fig. 7(b)] are compressive throughout most of the cross section, with small tensile values in the extreme peripheral region. The values range from about 5 kPa tension to about 80 kPa compression and the maximum compressive

values are found in the central region of the cross section. Circumferential stress [Fig. 7(c)] is tensile to balance the radial component of the load and ranges from - 1800 kPa near the inner radius to over 2700 kPa near the outer radius. Circumferential fibers, therefore, are present to carry the significant circumferential loads, as observed in the original study of form and function by Shrive et al. (1978). Shear stress [Fig. 7(d)] is negative in the inner twothirds of the meniscus and increasingly positive in the outer third, reaching a peak value of 9 kPa near the apex of the femoral surface. Pressure [Fig. 7(e)] is also positive throughout the cross section and reaches a maximum value in excess of 300 kPa in the central region of the cross section. It is important to note that the pressure is much larger than the axial and radial stress throughout the cross section, and is similar in magnitude to the circumferential stress in the central region. Therefore, the pressure carries a significant portion of the applied load in these earlier times even though 75% of the applied load was placed on the solid phase. Such behaviors could not be captured by single-phase elastic or viscoelastic models. With time this pressure will decrease as the fluid flow slows and eventually ceases; in this process, the.biphasic model will also provide the time-dependent response. To demonstrate the flow of fluid within the meniscus, Fig. 8 shows the velocity of the fluid phase relative to the solid phase at the femoral surface for time t = 1.5 s. The magnitude of the relative velocity and of the individual solid and fluid phases (not shown here) is of the order of pm s - ‘. The data shown in Fig. 8 are a least-squares fifth-order fit (correlation coefficient, r > 0.96) of the nodal data, which shows some numerical oscillations about this curve fit. While the magnitude is small, the results in Fig. 8 indicate that the fluid is being exuded from the meniscus on the femoral surface (i.e. relative motion in the negative radial and positive axial directions). Similar results are found for t = 5 s (the time used in the previous contour plots). With only 25% of the load on the fluid, there is

A model of the meniscus

0.53,....,..,.,..,.,....,...., 0 1

2

1035

3

4

5

Time (set)

-0.21 . . 0

. , . . . . , . 1 2

. . , . 3

. . , . . 4

. ( 5

Time (set) Fig. 6. (a) Radial displacement at the outer radius of the tibia1 surface (point A) and (b) displacements at the midpoint of the femoral surface (point B) as a function of time for the normal meniscus simulation.

moderate resistance to fluid flow at the femoral surface. This fluid motion influences the solid-phase deformations as will be shown below and, as observed in various other biphasic problems, both analytically (Armstrong and Mow, 1982; Armstrong et al., 1984; Mow et al., 1984, 1986; Holmes et al., 1985) and numerically (Spilker et al., 1988, 1990; Suh, 1989; Spilker and Suh, 1990; Suh et al., 1991). Figures 9(a)-(d) show the distribution of each component of strain throughout the meniscus cross section at t = 5 s. For the present simulations, we focus attention on those regions of the cross section which contain the highest tensile strains. Recall that the magnitude of the strain can be expected to increase with increasing time, and to increase roughly in proportion to the magnitude of the load. The increase is precisely linear with increasing load within the context of linear theory, but is not linear for the more

realistic nonlinear theories which would be required at these increased load levels. The axial strain, sz, is compressive throughout most of the cross section [Fig. 9(a)], as would be expected for this compressive applied load. Tensile axial strains are found, however, throughout most of the bulged peripheral section, with a maximum value of 0.14. The radial strain, E, [Fig. 9(b)], is predominantly positive, although negative values are found in the bulged region. Except for a very small region at the outer radius of the tibia1 surface, where the strain reaches 0.15, the largest area of tensile strain is in the central portion of the meniscus, where values range from 0.05 to 0.1. The circumferential strain, sO, shown in Fig. 9(c), is tensile throughout the cross section, and is relatively uniform through the thickness. The maximum value of 0.013 is reached at a radial position corresponding approximately to the outer radius of

1036

R. L. SPILKER et al.

(a)

-3.5 -19.5 -35.4 -51.4 -67.3 -83.3 kPa

(bf

-3.12 -17.0 -30.8 -44.7 -58.5 -72.3 kPa

(cl

2300 1560 816 74 -668 -1410 kPa

kPa Fig. 7(a-d).

A model of the meniscus

1037

(e) n 302 247

Fig. 7. Distribution of the (a) axial, (b) radial, (c) circumferential and (d) shear stress components, and (e) pressure for the normal (average)meniscal simulation at time t = 5 s after application of the load. Lines of constant stress/pressure are indicated.

_41”‘I 0

.‘.l.. 0.2

7,. . .I.. 0.4 0.6 0.8 Edge Parameter, 5

., 1

Fig. 8. Radial and axial components of the velocity of the fluid relative to the solid at the femoral surface for the normal meniscus at t = 1.5 s. Values are plotted vs nondimensional length, 5, from the inner radius (5 = 0) to the outer radius (5 = I).

the tibia1 surface. These relatively smaller values of circumferential strain, compared with the other components of strain, are a result of the stiffening effect of the circumferential fibers. The shear strain, y,. [Fig. 9(d)], is positive over a large region of the meniscus. Maximum positive values between 0.26 and 0.37 are found near the outer radius just below the apex. These results provide insight into possible mechanisms for meniscal failure. In particular, we hypothesize that the bucket-handle tear could be caused by large tensile strain, via a single large component or a combination of tensile components. Experimental results (Proctor et al., 1989) for meniscal tissue suggest that specimens aligned with the fiber direction fail at the order of 0.05 tensile strain, whereas specimens oriented perpendicular to the fiber direction fail at the order of 0.3 strain. In the present simulation of a normal meniscus, the region of maximum normal strain BM25:9-G

for each of the strain components is in the mid to outer third of the meniscal radius, where the shear strain is also high. In view of the observed combination of maximum tensile strains, and given that the strain values presented in Figs 9(a)-(d) will increase with time and/or with increased load level, the present results suggest that the mid to outer third of the meniscal geometry is most susceptible to the buckethandle mode of failure. The largest (and smallest) normal strains at a point in a continuum correspond to the largest (and smallest) principal strains at that point. For an axisymmetric representation, one of the principal strains is in the circumferential direction and the other two principal strains are in the plane of the cross section (i.e. in the r-z plane). Because the material response parallel to the fibers (circumferential) and perpendicular to the fibers (in the plane of the cross

1038

R. L. SPILKERet al.

-0.008 -0.064 -0.12

0.0115 0.0078 0.0041 0.00044 -0.0033 -0.0069

(d)

0.31 0.19 0.07 -0.05 -0.17 -0.29

Fig. 9. Distribution of the (a) axial, (b) radial, (c) circumferential and (d) shear strain components for the normal (average) meniscal simulation at time t = 5 s after application of the load. Lines of constant strain are indicated.

A model of the meniscus section) is quite different,

we are, therefore,

interested

in the location of maximum tensile circumferential strain, and the maximum tensile value of the principal strains in the plane of the cross section. Any tensile strain in the cross section, regardless of its direction in the cross section, could be the source of a circumferential bucket-handle tear, for example. The circumferential strain has already been plotted [Fig. 9(c)]. Shown in Fig. 10 is the distribution of the maximum principal strain in the cross section. Note that the maximum normal strain is positive at all locations in the cross section, with a maximum value of 0.37. This distribution is quite similar to the shear strain distribution shown in Fig. 9(d), demonstrating the important contribution of shear strain to the maximum normal strain. This result suggests that when the load is distributed over the entire femoral surface, the entire cross section is susceptible to failure, but that the region near the femoral surface in the outer third is the most susceptible. Effects of boundary conditions at the tibia1 ana’ femoral surfaces

The loading conditions at the femoral surface on displacement, velocity, and load at the tibia1 and femoral surfaces represent approximations of the more realistic three-body contact problem involving the meniscus, femur, and tibia. It is of interest, therefore, to determine to what extent changes in the load sharing on the femoral surface and fluid flow at the tibia1 surface affect the most significant quantities of interest in the present analysis. The boundary condition on the rigid tibia1 surface can be set to allow for either an impermeable interface (as used in the previous simulations) or a free-draining interface, in which case the fluid is free to flow into and out of the meniscus at this interface. Alternatively, a pressure can be applied to the fluid phase at the tibia1 surface, to represent the effect of the thin layer of fluid between the meniscus and tibia, and the fluid movement at the tibia1 surface left unspecified. The load applied to the femoral surface can be

1039

partitioned between the solid and fluid phases. If the load is applied entirely to the solid phase, similar to the situation when a porous indenter is used in testing soft hydrated tissue, the fluid pressure at the femoral surface is zero, and the fluid is free to flow across the boundary. This loading condition is most closely matched with the free-draining tibia1 boundary. If the load is transmitted to the biphasic material entirely through the fluid layer, then the load is partitioned according to & on the solid phase and 4’ on the fluid phase (for the case considered in this study, 25% of the load would be applied to the solid phase, and the remaining 75% to the fluid phase) (Hou et al., 1989). The two extreme cases of femoral/tibia1 loading and boundary conditions could, thus, be defined as (i) a free-draining tibia1 interface with the load on the femoral surface applied completely to the solid phase, and (ii) an impermeable tibia1 interface with femoral surface loading split according to 4” and 4’. A number of intermediate cases are also considered, as summarized in Table 2. Case 5, whichinvolves a back pressure at the tibia1 surface [Fig. 3(b)], requires some additional explanation. To determine the magnitude and distribution of the pressure on the tibia1 surface, the distribution of the total axial stress (sum of the solid and fluid stresses) was calculated along the tibia1 surface for Case 4 (note that the integral of this distribution over the tibia1 surface is equal to the net axial force applied to the femoral surface). We found that the distribution of total axial stress was reasonably approximated by a parabolic distribution and, as expected, that the distribution was essentially constant in time. We further found that the total axial stress is compressive at all locations on the tibia1 surface, suggesting that there is no tendency in the present analysis for the meniscus to leave the tibia1 surface. From theoretical developments (Hou et al., 1989), if this total load is transmitted to the meniscus through the fluid layer, it should be partitioned according to 4’ and 6. We, therefore, multiplied the parabolic distribution of total axial stress by 4’ and applied this pressure to the fluid phase at the tibia1 surface. No load was applied to the

0.34 0.272 0.204 0.136 0.068 0.0

Fig. 10. Principal strain for the normal meniscus at t = 5 s.

1040

R. L. SPILKERet al. Table 2. Description of the boundary conditions for several cases in the meniscal simulation

Case number

Case name

Femoral surface loading on the solid phase (%)

Normal 25% impermeable 25% permeable 100% impermeable 25% back pressure 75% permeable 100% permeable

15 25 25 100

25 75 100

Tibia1 surface boundary condition* Impermeable Impermeable Permeable Impermeable Fluid pressure prescribed Permeable Permeable

* Axial solid displacement is zero for all cases.

solid phase. The solid-phase displacement perpendicular to the tibia1 surface was set to zero, but no constraints were imposed on the IIuid-phase velocities at the tibia1 surface (as in all analyses thus far). To demonstrate the effect of these boundary conditions, the displacement at the middle of the femoral surface (both axial and radial components) and the radial displacement at the outer radius of the tibia1 surface are plotted vs time for the first 10 s of response in Figs 1l(a)-(c). The extreme cases, the ‘normal’ meniscus case considered previously, and several intermediate cases are considered. Of the quantities examined, the largest displacement is found at the outer radius of the tibia1 surface, but this quantity is minimally affected by the changes in boundary conditions. Those cases with an impermeable boundary condition at the tibia1 surface showed an increase in displacement with time, while the permeable cases exhibited a decrease with time. The manner in which the load on the femoral surface was partitioned had little effect on this particular quantity. Results such as these are indicative of the role of interface permeability at times soon after the application of step loadings. Radial displacements at the midpoint of the femoral surface are more significantly affected by the proportion of the load applied to the solid phase, and hence the permeability of that surface, than by the tibia1 surface permeability. The amount of creep observed increases as the percentage of load on the solid phase increases. Here again, local motion of the fluid effects the solid displacement. The largest effect, in both an absolute and relative sense, is found for the axial displacement at the middle of the femoral surface. While the radial displacements examined above were effected only by the permeability of the boundary in question, the axial displacement is influenced by the boundary conditions at both the tibia1 and femoral surfaces. Within the @st 10 s there is almost no change in the axial displacement when the load is partitioned according to 9’ and 4’ and the tibial interface is impermeable. As more load is transferred to the solid phase and/or when the tibia1 interface is made permeable, a larger increase in the displacement with time is noted.

In general, the appearance of the stress and strain distributions for the cases listed in Table 2 are similar, but there are some differences in the relative magnitudes observed. The extreme cases (Cases 2 and 7) and that in which a fluid pressure was applied at the tibia1 surface (Case 5) are of most interest. With 100% of the load on the solid phase and a permeable tibia1 surface, the magnitudes of compressive strain are larger than in the normal case; there were no increases in the tensile strains. This is not a significant elect, since tensile strains are thought to be responsible for failures. Figure 12 shows that principal strains in the cross section are compressive in the inner half of the sample, where they were tensile at all locations in the normal case. Both the shear stress and shear strain have local regions at the center of the femoral surface, where values are near the maximum tensile value observed in the cross section. The other extreme case (Case 2), when the tibia1 surface is impermeable and only 25% of the load is applied to the solid phase at the femoral surface, provides very little change in either the appearance or magnitudes of the distributions. Principal strains are all tensile and range up to a maximum of 0.38 at the outer third of the femoral surface. Peak values of the circumferential stress are nearly 10% higher in compression than those observed for the normal case, though. Case 5 has a permeable tibia1 surface, but imposes a fluid pressure there to model a fluid film interface. Some compressive principal strains exist in the inner third of the sample, but the peak tensile values are still near 0.38. Very high compressive axial strains exist at the center of the tibia1 surface and lower pressures are observed in this region, as compared to the normal case. The locations and peak values of the maximum tensile stress and strain components are unchanged.

SUMMARY AND CONCLUDING REMARKS Axisymmetric finite element simulations of meniscal response, based on a biphasic description of the tissue and a transversely isotropic representation of

1041

A model of the meniscus

-. 0

A-. --.. 0

0

0

.a.._ A’....*

0

-... ‘b’...

0

Oh

ii

0

9

0

0.54 -_,I.,,..,.,.,...,...,

0

0.535

O

a

2

Time 4 (set)

6

8

10

0.11, ..: .... 0 ...'

- '-26% /Impermerble 0 26%IPenne&le 0 1OOU Impermeabla o 26% IBuk Pm~~um A 76% IPermeable

*..-

.tj’

*... /ti ...*

3

..6’

!l 0*09

.*ei

.d

/

(d

-E 3

0.08

b

0.07 0

,,, I,,,(, ,,,... ,,..# 2

4

6

8

10

Time (set)

z

-0.25

s

c

O

2

4

Time (set)

6

8

lo

Fig. 11. Comparisons of the (a) radial displacement at the outer radius of the tibia1 surface (point A) and the (b) radial and (c) axial displacements at the midpoint of the femoral surface (point B) for several different tibia1 and femoral surface boundary conditions. The displacement at time t = 0 s is zero. Case names refer to Table 2.

1042

R. L.

SPILKER et al.

n 0.32 0.23 0.14 0.044 -0.048 -0.14

Fig. 12. Principal strains at t = 5 s for Case 7, 100% permeable.

the solid phase, have been presented. The response of a ‘normal’ meniscus of average geometry and material properties, resting on a smooth and impermeable tibia1 interface and loaded through the femoral surface, has been examined. The effects of changes to the modeling of the boundary condition at the tibia1 surface, and the partitioning of the load on the femoral surface have also been examined. The response of the meniscus within the first 5 s has been investigated, with emphasis on the distributions of stress, pressure, strain, and fluid flow throughout the cross section. A primary objective was the identification of regions of the cross section having large tensile strains. These large tensile strains, individually or in combination, can lead to failures of the meniscus such as the bucket-handle tear. On the basis of the present studies, we make the following observations. -The circumferential stresses are dominant, while the circumferential strains are relatively small. The transversely isotropic representation of the solid phase of the meniscal tissue is, therefore, essential in order to capture this fundamental, anisotropic behavior. -The pressure of the fluid phase is larger than either axial or radial normal stress and carries a significant part of the load at early times after the onset of loading. This behavior could not be captured by single-phase models and, thus, the inclusion of the biphasic model of tissue is essential. -In general, the axial strain is compressive, except for the region near the bulged peripheral surface. The radial strain is positive in a large central region of the cross section, and the circumferential strain is positive throughout the cross section. Large shear strains are found, particularly near the apex of the cross section. -The largest principal strains, representing the largest normal strains, are positive throughout the cross section and are maximum in regions of maximum shear strain. In the present simulations, the largest values are found in the region near the femoral surface and in the outer third of the cross section. -Because of the stiffening effect of the circumferential fibers, the circumferential strains are the smallest

of the strains; however, these contributions may still be important since experimental studies (Proctor et al., 1989) show that tissue specimens parallel to the fiber directions can fail in extension at relatively low strain levels (e.g. 0.05). -When the fluid is allowed to flow through the tibia1 surface and/or the load along the femoral surface is partitioned more to the solid phase, permitting increased fluid flow at the boundary, the creep response of the tissue increases. As expected, the flow fields at these interfaces are affected; however, the general observations made above are essentially unaffected. -At loading levels exceeding body weight (twice the resting weight), very large deformations and strains occur in the tissue. It will, therefore, be important to perform future meniscal studies using a finite deformation biphasic theory which includes the transversely isotropic fibrous solid phase and probably strain-dependent permeability. Acknowledgements-The authors gratefully acknowledge the helpful discussions of meniscal geometry and material properties with Drs D. C. Fithian (now with the Department of Orthopaedics, University of California at San Diego) and K. Y. Chem of the Department of Orthopedic Surgery at Columbia-Presbyterian Medical Center. The support of the National Science Foundation (EET 87-20314), the National Institute of Arthritis, Musculoskeletal and Skin Diseases (AR 38728), and the Surdna Foundation are gratefully acknowledged.

REFERENCES

Ahmed, A. M. (1991) The load-bearing role of the knee menisci. In Knee Meniscus: Basic and Clinical Foundation (Edited by Mow, V. C., Arnoczky, S. P. and Jackson, b. W.), Ch. 6. Raven Press, New York. Ahmed, A. M. and Burke, D. L. (1983) In-vitro measurement of static pressure distribution in synovial joints-part I: tibia1 surface. of the knee. J. biomech. Enma 105, 216-225. Armstrong, C. G., Lai, W. M. and Mow; v. C. i1984) An analysis of the unconfined compression of articular cartilage. J. biomech. Engng 106, 165-173.

A model of the meniscus Armstrong, C. G. and Mow, V. C. (1982) Variations in the intrinsic mechanical properties of human articular cartilage with age. degeneration, and water content. J. Bone Jt &rg. 64A; 88-G. Arnoczky, S. A., Adams, M., DeHaven, K., Eyre, D. R. and Mow, V. C. (1988) Meniscus. In Injury and Repair of the Musculoskeletal Soji Tissue (Edited by Woo, S. L-Y. and Buckwalter, J. A.), pp. 483-537. American Academy of Orthopaedic Surgeons, Chicago. Aspden, R. M. (1985) A model for the function and failure of the meniscus. Engng Med. 14, 119-122. Aspden, R. M., Yarker, Y. E. and Hukins, D. W. L. (1986) The structure of the meniscus and its relationship to stability and failure. In Miniposters, 15th Symposium of the European Society of Osteoarthrology, Kuopio. Brown, T. D. and Singerman, R. J. (1986) Experimental determination of the linear biphasic constitutive coefficients of human fetal proximal femoral chondroepiphysis. J. Biomechanics 10, 597-605. Bullough, P. G., Munuera, L., Murphy, J. and Weinstein, A. M. (1970) The strength of the menisci of the knee as it relates to their fine structure. J. Bone Jt Surg. 52B, 564. Chern, K. Y., Zhu, W. B., Kelly, M. and Mow, V. C. (1990) Anisotropic shear properties of bovine meniscus. In Transactions, 36th Orthopaedic Research Society, p. 246. Chern, K. Y., Zhu, W. B. and Mow, V. C. (1989) Anisotropic viscoelastic shear properties of meniscus. In 1989 Aduances in Bioenaineerina (Edited by Rubinsky, B.), BED_ 15, pp. 105-106: ASME, New York. Fithian, D. C., Kelly, M. A. and Mow, V. C. (1990) Material properties and structure-function relationships in the menisci. Clin. Orthop. Rel. Res. 252, 19-31. Fithian, D. C., Zhu, W. B., Ratcliffe, A., Kelly, M. and Mow, V. C. (1989) Exponential law representation of tensile properties of human meniscus. Proc. Inst. Mech. Engng in Bioengng, pp. 85-90. Hefzy, M. S., Grood, E. S. and Zoghi, M. (1987) An axisymmetric finite element model of the meniscus. In 1987 Advances in Bioengineering (Edited by Erdman, A.), BED-3, pp. 51-52. ASME, New York. Hefzy, M. S. and Zoghi, M. (1988) The meniscus: a finite element model. In Transactions, 34th Orthopaedic Research Society, p. 149. Holmes, M. H. (1986) Finite deformation of soft tissue: analysis of a mixture model in uni-axial compression. J. biomech. Engng 108, 372-381. Holmes, M. H., Lai, W. M. and Mow, V. C. (1985) Singular perturbation analysis of the nonlinear, flow-dependent compressive stress relaxation behavior of articular cartilage. J. biomech. Engng 107, 206-218. Hou, J. S., Holmes, M. H., Lai, W. M. and Mow, V. C. (1989) Boundary conditions at the cartilage-synovial fluid interface for joint lubrication and theoretical verifications. J. biomech. Engng 111, 78-87. Kwan, M. K., Lai, W. M. and Mow, V. C. (1990) A finite deformation theory for cartilage and other soft hydrated connective tissues. J. Biomechanics 23, 145-155. Lai, W. M., Hou, J. S. and Mow, V. C. (1989) A triphasic theory for articular cartilage swelling and Donnan osmotic pressure. In Biomechanics Symposium, University of California, San Diego. ASME, New York. Lai, W. M., Hou, J. S. and Mow, V. C. (1991) A triphasic theory for the swelling and deformational behaviors of articular cartilage. J. biomech. Engng 113(3), 245-258. Lai, W. M. and Mow, V. C. (1980) Drag-induced compression of articular cartilage during a permeation experiment. Biorheology 17, 111-123. Mak, A. F. (1980) The apparent viscoelastic behavior of articular cartilage-the contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J. biomech. Engng 108, 123-130. Mak, A. F., Lai. W. M. and Mow, V. C. (1987) Biphasic indentation of articular cartilage: Part I, theoretical analysis. J. Biomechanics 20, 703-714.

1043

Mow, V. C., Gibbs, M. C., Lai, W. M., Zhu, W. B. and Athanasiou, K. A. (1989) Biphasic indentation of articular cartilage+II. A numerical algorithm and an experimental study. J. Biomechanics 22, 853-861. Mow, V. C., Holmes, M. H. and Lai, W. M. (1984) Fluid transport and mechanical properties of articular cartilage: a review. J. Biomechanics 17, 377-394. Mow, V. C., Kuei, S. C., Lai, W. M. and Armstrong, C. G. (1980) Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments. J. biomech. Engng 102, 73-84.

Mow, V. C., Kwan, M. K., Lai, W. M. and Holmes, M. H. (1986) A finite deformation theory for nonlinearly permeable cartilage and other soft hydrated connective tissues. In Frontiers in Biomechanics (Edited by Woo, S. L-Y., Schmid-Schonbein, G. and Zweifach, B.), pp. 153-179. Springer, New York. Niu, Q. and Shephard, M. S. (1990) Transfer of Solution Variables for Finite Element Meshes. Scientific Comoutation Reseirch Center, Rensselaer Polytechnic Ins&ute, Troy, NY. Oomens, C. W. J., Van Campen, D. H. and Grootenboer, H. J. (1987) A mixture approach to the mechanics of skin. J. Biomechanics 20, 877-885.

Proctor, C., Schmidt, M. B., Whipple, R. R., Kelly, M. A. and Mow, V. C. (1989) Material properties of the normal medial bovine meniscus. J. orthop. Res. 7, 71-82. Salzstein, R. A. and Pollack, S. R. (1987) Electromechanical potentials in cortical bone--II. Experimental analysis. J. Biomechanics 20, 271-280.

Salzstein, R. A., Pollack, S. R., Mak, A. F. T. and Petrov, N. (19871 Electra-mechanical ootentials in cortical bone--I. k continuum theory. J. Biomechanics 20, 261-270. Sauren, A. A. H. J., Huson, A. and Schouten, R. Y. (1984) An axisymmetric finite element analysis of the mechanical function of the meniscus. Int. J. Sports Med. 5, 93-95. Schreppers, G. J. M. A., Sauren, A. A. H. J. and Huson, A. (1990) A numerical model of the load transmission in the tibio-femoral contact area. J. Engng Med., Proc. Instn Mech. Engrs 204, 53-59.

Shephard, M. S., Baehmann, P. L. and Grice, K. R. (1988) The versatility of automatic mesh generators based on tree structures and advanced geometric constructs. Comm. appl. num. Meth. 4, 379-392.

Shrive, N. G., O’Connor, J. J. and Goodfellow, J. W. (1978) Load bearing in the knee joint. Clin. Orthop. 131, 279. Skaggs, D. C. and Mow, V. C. (1990) Function of radial tie fibers in meniscus. In Transactions, 36th Orthopaedic Research Society, p. 248. Spilker, R. L. and Daugirda, D. M. (1981) Analysis of axisymmetric structures under arbitrary loading using the hybrid stress model. Int. J. num. Meth. Engng 17,801-823. Spilker, R. L. and Maxian, T. A. (1990) A mixed-penalty finite element formulation of the linear biphasic theory for soft tissues. Int. J. num. Meth. Engng 30, 1063-1082. Spilker, R. L. and Suh, J.-K. (1990) Formulation and evaluation of a finite element model of soft hydrated tissue. Camp. Struct. 35, 425439.

Spilker, R. L., Suh, J.-K. and Mow, V. C. (1988) A finite element formulation of the nonlinear biphasic model for articular cartilage and hydrated soft tissues including strain-dependent permeability. In Computational Methods in Bioengineering (Edited by Spilker, R. L. and Simon. B. R.), BED-9, pp. 81-92. ASME, New York. Soilker, R. L., Suh. J.-K. and Mow. V. C. 119901Effects of -friction on the unconfined compressive reHpon& of articular cartilage: a finite element analysis. J. biomech. Engng 112, 138-146. Suh, J.-K. (1989) A finite element formulation for nonlinear behavior of biphasic hydrated soft tissue under finite deformation, Ph.D. thesis, Rensselaer Polytechnic Institute, Troy, NY. Suh, J.-K., Spilker, R. L. and Holmes, M. H. (1991) A penalty finite element analysis for nonlinear mechanics of biphasic

1044

R. L.

SPILKER et al.

hydrated soft tissue under large deformation. ht. J. mm. Meth. Engng 32, 1411-1439. Suh, J. K., Spilker, R. L., Holmes, M. H. and Mow, V. C. (1989) A nonlinear biphasic finite element formulation for soft hydrated tissues under finite deformation. In 1989 Advances in Bioengineering, ASME

where E~~,si2, etc., and ull, u12, etc., are the strains and stresses, respectively, with respect to the material coordinate system and the material constants have been defined earlier. Due to symmetry of the stress and strain tensors, the following relations must be satisfied:

Winter Annual Meet-

ing (Edited by Rubinsky, B.). ASME, New York. Tissakht, M. and Ahmed, A. M. (1990) Effect of tibia1 axial rotation on knee meniscal stress: a finite element study. In Transactions, 36th Orthopaedic Research Society, p. 243. Tissakht, M., Farinaccio, R. and Ahmed, A. M. (1989) Effect of ligament attachments on the response of the knee menisci in compression and tdrsion: a finite element study. In

(A2a) “13

“31

E3

El

“32 -=-.

“23

Ez

E3

-=-

Transactions, 35th Orthopaedic Research Society Meeting,

WW

WC)

Since the solid of interest is transversely isotropic, we can make the following standard simplifying substitutions:

p. 203. Tissakht, M., Marchand, F. and Ahmed, A. M. (1991) Nonlinear finite element analysis of the knee menisci: a composite fiber-reinforced model. In Transactions, 37th Orthopaedic Research Society Meeting, 16, p, 294. Vermilyea, M. E. and Spilker, R. L. (1992) A hybrid finite element formulation of the linear biphasic equations for soft hydrated tissues. Int. J. num. Meth. Engng 33,

E3

=

E2r

(A34

“13

=

“129

WV

=

Glzt

G3

(A3c)

EZ

--

G23 - 2(1 + vz3)’

567-594.

Whipple, R. R., Wirth, C. R. and Mow, V. C. (1985) Anisotropic and zonal variations in the tensile properties of the meniscus. In Transactions, Orthopaedic Research Society, 10, p. 367.

For our axisymmetric problems, the nonzero components of strain (stress) are: E,, (a,), E,, (u,,), E@ (a#), and y,= (a,,). Thus, we do not need the full complement of stress-strain relations developed by Spilker and Daugirda (1981). By eliminating terms that we do not need for our problems, and by using the proper transformation relations, equation (Al) takes the followina form as the stress-strain relation in global coordinates?or our problems:

APPENDIX MATERIAL PROPERTY MATRIX FOR THE TRANSVERSELY ISOTROPIC SOLID PHASE OF A BIPHASIC MATERIAL

(A4) The properties of a transversely isotropic material are most appropriately defined in a material coordinate system. We, therefore, choose (see Fig. 5) x1 in the direction of the fiber orientation, x2 in the plane of the fibers but perpendicular to the fiber orientation, and xj in the ‘transverse’ direction perpendicular to the (x,, x2) plane of the fibers (note that xJ is in the radial direction); a is the fiber angle, measured counterclockwise positive from the Bdirection to the x1 direction. This local coordinate system (x1, x2, xp) can be mapped into the global coordinate system, (r, 0, z) through proper transformation laws. In order for the behavior of the meniscus to remain axisymmetric, G(can be 0 or 90”. Additionally, we include the case of a layer having a sufficient and equal number of + K and - K fibers so that the overall behavior is an ‘average’ of the individual + K and - a behaviors. In this context, K can take on any value. These assumptions were used by Spilker and Daugirda (1981), and eliminate all coupling terms which would produce nonaxisymmetric behaviors. We can define the stress-strain law for a transversely isotropic material, in material coordinates, as follows:

1

-

-

-

-

E2

1

“2.1 E,

-

"12

-

E,

“31

E, -

where s1, = $,

(A5a)

&=;+g+

_$Z+$ 1

(

s,3L+$+

12 >

1

12 >

++$(

s2c2,

s12

s

SZl

=

+

c2 +

$=

1 El

0

s=,

(A5e)

“23

2

(A5f)

E2

0

0

0

0

0

0

0

0

E3

-

-0

“23

1

E,

&

E,

0

0

0

0 I

0

0

o&o

0

0

0

00;

I 9

Glo 0

(A5c)

2

&3~s31=~s2+-c,

“13

(A5b)

(A5d)

s,,=&+&c2,

E3 “32

&2,

1

1

(Al)

A model of the meniscus

sz, = &

=

2

(s4+ c4)+ 1

S,d =

& +; - &

(

where B is the elastic straindisplacement matrix, and D1 and D2 are defined for an isotropic material so that

s2c2,

12 >

s,, = s,, = s,, = s34 = s*j = 0,

D = 1*D, + $D2

WV

is the material property matrix. For the transversely isotropic solid phase, D is given by the inverse of S and the sohd-phase element stiffness is given by

-flpI+D~~.

(2d) For a transversely isotropic solid phase, D is given by the inverse of the S matrix just defined [equation (A4)]. This requires a minor modification to the matrix definitions derived for the biphasic mixed-penalty method (Spilker and Maxian, 1990) and used in this study. The solid-phase element stiffness matrix defined in this reference [equation (18b) of the reference] is (I’BrD,B + /PBTD2B)dQ,

k= sn

iA

Wg)

where c = cos a and s = sin a. The transversely isotropic extension of the biphasic representation requires a simple modification that affects only the solid phase. The constitutive equation for an isotropic solid was given by equation (2d): u’=

1045

(At3

(Ag) This transversely isotropic formulation requires that we supply a, El, El, G12, v12, vz3, @, and K. Our finite element program constructs the compliance matrix, S, using equations (A4) and (AS). Then, the program analytically inverts the compliance matrix to obtain the needed anisotropic material property matrix, D = S-i. Since there are several zeros in the compliance matrix [equation (A4)], the program can calculate this inverse with minimal computational effort.