Numerical study of temperature effects on the poro-viscoelastic behavior of articular cartilage

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

Contents lists available at ScienceDirect

Journal of the Mechanical Behavior of Biomedical Materials journal homepage: www.elsevier.com/locate/jmbbm

Numerical study of temperature effects on the poro-viscoelastic behavior of articular cartilage

T

⁎

Reza Behrou , Hamid Foroughi, Fardad Haghpanah Department of Civil Engineering, Johns Hopkins University, Baltimore, MD, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Viscoelasticity Poro-viscoelasticity Thermo-poro-viscoelasticity Articular cartilage Finite element method

This paper presents a new approach to study the effects of temperature on the poro- elastic and viscoelastic behavior of articular cartilage. Biphasic solid-fluid mixture theory is applied to study the poro-mechanical behavior of articular cartilage in a fully saturated state. The balance of linear momentum, mass, and energy are considered to describe deformation of the solid skeleton, pore fluid pressure, and temperature distribution in the mixture. The mechanical model assumes both linear elastic and viscoelastic isotropic materials, infinitesimal strain theory, and a time-dependent response. The influence of temperature on the mixture behavior is modeled through temperature dependent mass density and volumetric thermal strain. The fluid flow through the porous medium is described by the Darcy's law. The stress-strain relation for time-dependent viscoelastic deformation in the solid skeleton is described using the generalized Maxwell model. A verification example is presented to illustrate accuracy and efficiency of the developed finite element model. The influence of temperature is studied through examining the behavior of articular cartilage for confined and unconfined boundary conditions. Furthermore, articular cartilage under partial loading condition is modeled to investigate the deformation, pore fluid pressure, and temperature dissipation processes. The results suggest significant impacts of temperature on both poro- elastic and viscoelastic behavior of articular cartilage.

1. Introduction Diarthodial joints are the most common movable skeletal joints which are specifically characterized by some common structural features including: a layer of fibrocartilage that covers the opposing bony surfaces; a lubricating synovial fluid within the joint cavity; and an enclosing fibrous capsule which is lined with an active tissue (i.e. the synovium) (Mow et al., 1993). Fig. 1 illustrates the structure of knee joint as an example of diarthodial joints. Articular cartilage is a thin connective fibrocartilage (Fig. 1) that primarily consists of water, collagens, and proteoglycans. With 80% of the wet weight, water is the dominant constituent of articular cartilage. The network of collagen, which constitutes 60% of the dry weight, together with water increase the strength of the tissue to transmit high instantaneous loads. For smooth joint articulation, articular cartilage provides a lubricated surface that facilitates load transmission to the underlying subchondral bone (Sophia Fox et al., 2009). The mechanical behavior of articular cartilage has been studied using different models. Pena et al. (2005) studied the effect of meniscal tears and meniscectomies on knee joints considering cartilage as a linear elastic isotropic and homogeneous material. Viscoelastic models

⁎

are widely used to study the influences of an interstitial fluid on soft biological tissues (Mak, 1986; Shirazi et al., 2008; Gupta et al., 2009; Gu and Li, 2011). Poro-elastic and poro-viscoelastic models are developed to account for fluid flow in articular cartilage where biphasic mixture theory is used to describe deformation of the solid and the solid-fluid interactions (Huang et al., 2003; Netti et al., 2003; Hoang and Abousleiman, 2009; Tomic et al., 2014). While most current models have studied the poro- elastic and viscoelastic behavior of articular cartilage, experimental studies (June and Fyhrie, 2010) revealed that increasing temperature affects the behavior of both fluid and solid phases. Increasing temperature increases equilibrium stiffness and stress relaxation at higher temperatures, and decreases fluid viscosity. A comprehensive understanding of articular cartilage mechanics requires considering the effects of temperature. This study aims to investigate the influence of temperature on the poroelastic and viscoelastic behavior of articular cartilage. In this paper, a finite element model is developed to explore the thermal effects coupled with the poro- elastic and viscoelastic behavior of articular cartilage. The poro-mechanical behavior of a fully saturated articular cartilage is modeled through the solid-fluid mixture theory. The mechanical deformation of the solid skeleton, fluid flow, and

Corresponding author. E-mail address: [email protected] (R. Behrou).

https://doi.org/10.1016/j.jmbbm.2017.11.023 Received 19 October 2017; Received in revised form 6 November 2017; Accepted 13 November 2017 Available online 22 November 2017 1751-6161/ © 2017 Elsevier Ltd. All rights reserved.

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

nα (x α , t ) =

dv α , dv

dv = J α dV α ,

J α = det F α ,

(1)

where dV α is the initial differential volume, dv α is the current differential volume, dv is the current total differential volume of the mixture, J α is the Jacobian of the deformation matrix, and F α is the deformation gradient (Boer, 2006). In biphasic mixture theory, ns + nf = 1 where ns and nf are volume fractions of the solid and fluid phases, respectively. The real ( ρ αR ) and partial ( ρ α ) mass densities of the α phase and the total mass density, ρ , are defined as follows:

dmα , dv α α dm = ρ αR (x α , t ) nα (x α , t ), ρ α (x α , t ) = dv ρ (x α , t ) = ρ s (x α , t ) + ρ f (x α , t ),

ρ αR (x α , t ) =

(2)

where dmα is the differential mass. The material time derivative of any spatial field, ψα (x α , t ) , for the α phase can be split into the local and convected terms as follows (Boer, 2006):

Fig. 1. Schematic diagram of the knee joint.

D α ψα (x α , t ) ∂ψα (x α , t ) ∂x α ∂ψα (x α , t ) + = Dt ∂x α ∂t ∂t ∂ψα + ∇ψ α · v α , = ∂t

temperature distribution in the solid-fluid mixture are modeled by the balances of linear momentum, mass, and energy. The mechanical behavior of the solid skeleton is described by both linear elastic and viscoelastic isotropic materials, infinitesimal strain theory, and a transient response. The motion of fluid through the porous medium is described by the Darcy's law. The influence of temperature on the mixture behavior is modeled by temperature dependent mass density and volumetric thermal strain. The governing equations are discretized in time with an implicit backward Euler scheme and in space by a standard Galerkin finite element method. The accuracy and efficiency of the finite element model are evaluated using analytical reference solutions. The remainder of the paper is organized as follows: in Section 2, the methodology for theoretical formulation of the thermo-poro- elastic and viscoelastic models is elaborated. The numerical implementation and finite element analysis are described in Section 3. In Section 4, the effects of material type, temperature, and boundary conditions on the behavior of articular cartilage are discussed. Insights gained from the results and areas for future research are summarized in Section 5.

where

vα

(3)

is the velocity vector.

2.2. Balance of mass For the α phase, the balance of mass is defined as follows:

Dαρα + ρα ∇·v α = γ α, Dt where

Dα ρα Dt

(4)

is the material time derivative of ρ α , γ α is the mass supply, ∂ (•)

and ∇ ·(•) = ∂x i is the divergence of a field. Considering volumetric i thermal expansion for the α phase, the temperature dependent real mass density is expressed as follows (Lewis and Schrefler, 1999):

1 D α ρ αR 1 D α pα Dαθα = α − 3αθα , αR ρ Dt K Dt Dt

2. Methodology

Kα

(5)

αθα

pα

where is the bulk modulus, is the pressure, is the thermal expansion coefficient, and θ α is the temperature field. From material time derivative of the fluid motion with respect to the motion of the solid phase, the following equations are derived (Boer, 2006):

The proposed thermo-poro-mechanical model predicts the coupled thermo-poro- elastic and viscoelastic behavior of articular cartilage. The model, developed based on the continuum theory of mixtures, describes the mechanical deformation, pore fluid pressure and thermal phenomena in a fully saturated articular cartilage. For more information about the basic theory of mixtures and continuum mechanics of porous media the reader is referred to Boer (2006). The displacement of solid skeleton is modeled by the balance of linear momentum. The flow of fluid though the porous medium is described by the Darcy's law. The temperature evolution in the solid-fluid mixture is modeled by the balance of energy. The model for the thermo-poro- elastic and viscoelastic behavior of the saturated articular cartilage is described below.

D f nf D s nf vf, = + ∇nf ·∼ Dt Dt D f pf D spf vf, = + ∇p f ·∼ Dt Dt Df θf D sθ f vf, = + ∇θ f ·∼ Dt Dt

(6)

v f = v f − v s is the relative velocity of the fluid phase with rewhere ∼ spect to the motion of the solid phase. Use of Eqs. (2), (5) and (6) in (4), the balance of mass for the mixture, with respect to the solid phase motion, which accounts for the volumetric thermal expansion in both solid and fluid phases, is expressed as follows:

2.1. Basic definition of mixture

1 ns D s p s nf D s p f v f) + f + f ∇p f ·(nf ∼ K s Dt K Dt K D sθ s D sθ f v f⎞ − 3nsαθs − 3nf αθf ⎛ + ∇θ f ·∼ Dt ⎝ Dt ⎠

Consider a porous medium consisting of a porous solid phase, denoted by s, with pore spaces occupied by a fluid, denoted by f . The medium, therefore, is characterized by two phases. At time t, the current position of an α phase is defined by x α where α represents either the solid or the fluid phase. Throughout the paper, any variable appearing in the form of (•)α represents the (•) property of the α phase. The volumetric fraction occupied by the α phase is defined by:

⎜

γs γf v f ) = sR + fR . + ∇ ·v s + ∇ ·(nf ∼ ρ ρ 215

⎟

(7)

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

∼ σ s = σ ′ − (B − nf ) p f I − σ f,visc.

2.3. Balance of linear momentum

For thermo-elastic problems, the total strain tensor of the solid skeleton is decomposed into the elastic and thermal components as follows:

The balance of linear momentum for the α phase is defined as:

∇ ·σ α + ρ α bα + hα = ρ α aα + γ α v α ,

σα

(8)

bα

aα

is the partial stress, is the body force vector, and is the where acceleration. The total stress is defined as σ = σ s + σ f . The internal body force due to drag, hα , on the α constituent caused by the other constituents is given by:

∑

ε s = ε skel,e + ε s,thm,

D sε s Df ε f D s ε skel, e D sθ s + σf: = σ ′: + αθs σ ′: I Dt Dt Dt Dt ∼ v f + σ f,visc: ∇∼ vf. − B p f ∇ ·v s − nf p f ∇ ·∼

(9)

ψ s = ψ s (ε skel, e, θ s ), 2.4. Balance of energy

Dαeα ρα − σ α : ε˙ α + ∇ ·qα − ρ α r α Dt 1 − γ α ⎛ v α ·v α − e α ⎞ + hα ·v α − e α̂ = 0, ⎝2 ⎠ D α ηα 1 − ρ α r α − α ∇θ α ·qα + ∇ ·qα ≥ 0, Dt θ

ε˙ α

∂ 2ψ s ∂θ s ∂ε skel, e

(21)

(11) (12)

s s skel, e ⎛σ ′ − ρ s ∂ψ ⎞: D ε skel, e ∂ε Dt ⎠ ⎝ s ψ ∂ D sθ s s s − ⎛ρ + ρ η + 3nsαθs p f − αθs σ ′: I⎞ ⎝ ∂θ ⎠ Dt

(13)

Use of the material time derivative of the Helmholtz free energy and the first law of thermodynamics in the second law of thermodynamics leads to the Clausius-Duhem inequality for the α phase as follows:

1 γ α ⎛ηα θ α + v α ·v α − e α ⎞ + σ α : ε˙ α − hα ·v α + e α̂ 2 ⎝ ⎠ D α ψα 1 Dαθα − ρα − ρ α ηα − α ∇θ α ·qα ≥ 0. Dt Dt θ

σ′ = ρs ρ s ηs

=

∂ψ s

,

∂ε skel, e ∂ψ s −ρ s ∂θ

ρ f ηf = −ρ f

∂ψ f ∂θ

− 3nsαθs p f + αθs σ ′: I , − 3nf αθf p f .

(24)

Thus, Eq. (23) reduces to:

1 − (∇p f − ρ fR g)·∼ v fD − ∇θ ·q ≥ 0, (25) θ f ∼ ∼ f f where v D = n v is the Darcy relative fluid velocity. In Eq. (25), the constitutive equation for the Darcy relative fluid velocity gives nonnegative dissipation (Coussy et al., 2004), where:

where the partial fluid pressure is related to the Cauchy pore fluid pressure through p f = nf p f . The fluid viscous stress is defined as Holzapfel (2000):

1 (∇v f + (∇v f )T ), 2

(23)

The Inequalities in Eq. (23) must hold valid for all possible thermodynamic states. Following Coleman and Noll (1963) argument for independent processes (i.e. D s ε skel, e / Dt , D s θ / Dt , and, D f θ / Dt ), the following constitutive equations must hold:

(14)

(16)

df =

(22)

∂ψ f Df θ − ⎛⎜ρ f + ρ f ηf + 3nf αθf p f ⎟⎞ ⎝ ∂θ ⎠ Dt 1 ∼ f fR f f − (∇p − ρ g)·n v − ∇θ ·q ≥ 0. θ

In the mixture theory, the total Cauchy stress, σ , is defined as a function of the effective Cauchy stress, σ ′, acting on the solid skeleton, and the Cauchy pore fluid pressure, p f : ∼ σ = σ ′ − B p f I, (15) ∼ s s where B = 1 − K skel/ K ≤ 1 is the Biot coefficient, and K skel and K are the bulk moduli of the porous medium and solid phase, respectively (Gawin et al., 1996). The fluid partial stress, σ f , is described by fluid pressure and viscous terms as Holzapfel (2000):

σ f = −p f I + σ f,visc,

Cpα ∂ 2ψ α =− α, α 2 ∂ (θ ) θ

where is the specific heat capacity per unit mass. We further assume that (1) both solid and fluid constituents are nearly incompressible, i.e. 1 ≈ 0 ; (2) the fluid flow is inviscid and the viscous effects are neKα glected, i.e. σ f,visc = 0 ; (3) there is no mass exchange between solid and fluid phases and the mass supply of the α phase is zero, i.e. γ α = 0 ; (4) the sum of the power of the solid and fluid phases is set to zero, i.e. e ŝ + e f̂ = 0 ; (5) θ s = θ f = θ ; (6) q = qs + qf ; (7) as = af = 0 ; (8) ∼ b s = b f = g ; and (9) B = 1. Thus, the Clausius-Duhem inequality of the solid-fluid mixture is expressed as:

eα

ψ α = e α − θ αηα .

= 0,

Cpα

where is the rate of the strain tensor, is the internal energy per unit mass, qα is the heat flux vector, r α is the heat input rate per unit mass, e α̂ is the power density supply by other phases on the α phase, and ηα is the entropy per unit mass. The Helmholtz free energy per unit mass is defined as Boer (2006):

kf

ψ f = ψ f (θ f ),

with the following assumptions:

The first and second laws of thermodynamics for the α phase are given by Boer (2006):

σ f,visc = k f (trdf ) I + 2μf df ,

(20)

With the assumption that both solid and fluid constituents are nearly incompressible (i.e. ρ αR is constant), the Helmholtz free energy per unit mass for the solid skeleton and the fluid phase is expressed as follows:

(10)

Furthermore, the balance of angular momentum states that the stresses are symmetric, i.e. σ α = (σ α )T .

γ αη α θ α + ρ α θ α

(19)

σ s:

By summing up the individual components of the balance of linear momentum for the solid and fluid phases, the balance of linear momentum for the mixture is expressed as:

∇ ·σ + ρ s b s + ρ f b f = ρ s as + ρ f af + γ s v s + γ f v f .

ε s,thm = αθs (θ s − θ0s ) I ,

where ε skel,e is the elastic mechanical strain, ε s,thm is the volumetric thermal strain, and θ0s is the initial temperature of the solid phase. Thus, the stress power for the solid-fluid mixture can be expressed as follows:

hα = hs + hf = 0.

α = s,f

(18)

∼ v fD = −kp (∇p f − ρ fR g),

(17)

(26)

and kp is the isotropic hydraulic permeability of the fluid. Finally, use of Eqs. (23) and (25) in (11) leads to the balance of energy for the solidfluid mixture as follows:

μf

is the partial bulk viscosity, is the partial shear viscosity, where and df is the symmetric portion of the velocity gradient for the fluid phase. Thus, the partial stress of the solid phase is derived as: 216

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

D sp f D sσ ′ D sθ :I − 3nsαθs θ + αθs θ Dt Dt Dt D sp f D sθ v f ⎞ − 3nf αθf θ ⎛⎜ v f ⎞⎟ + ∇p f ·∼ + ∇θ ·∼ + ρ f Cpf ⎛ Dt ⎠ ⎝ ⎝ Dt ⎠ v fD = 0. + ∇ ·q − ρ s r s − ρ f r f + (∇p f − ρ fR g)·∼

Δt hα (cmp) (t n + 1) = exp ⎜⎛− i ⎟⎞ hα (cmp) (t n ) ⎝ τm ⎠ ∞ (cmp) σ′ (t n + 1) − σ ′∞ (cmp) (t n ) ⎞ + ωmi τmi ⎛ Δt ⎠ ⎝

ρ s Cps

⎜

(27)

⎟

⎛ − exp ⎛− Δt ⎞ ⎞, ⎜ i ⎟⎟ ⎝ τm ⎠ ⎠ ⎝

⎜1

where i represents either the volumetric, k, or the deviatoric, μ , components of Eq. (31). For a given i, “cmp” represents the volumetric (i.e. cmp = vol for i = k ), the deviatoric (i.e. cmp = dev for i = μ ), and the thermal (i.e. cmp = thm for i = k ) components of Eq. (33). At any time t, the viscoelastic ISVs are recovered through a simple recursive update formulation and by knowing the ISVs from the previous time step.

2.5. Constitutive equations For an isotropic material, the effective stress of the solid skeleton in presence of the volumetric thermal strain can be expressed as:

σ ′ = 3K skel ε s,vol + 2μskel ε s,dev − 3K skel ε s,thm,

(28)

3. Finite element analysis

where μskel is the shear modulus of porous medium. The total mechanical strain of the solid skeleton is decomposed into the volumetric, ε s,vol , deviatoric, ε s,dev , and thermal, ε s,thm , components. For the viscoelastic behavior, a generalized Maxwell model shown in Fig. 2 is considered. The stress-strain relationship for the viscoelastic behavior of the solid skeleton can be characterized in a form of convolution integral:

A standard Galerkin finite element method is adapted to present the weak form of the governing equations. The displacement, u , the pore fluid pressure, p f , and the temperature, θ , are considered as independent state variables. The weak form of the governing equations is constructed by multiplying the strong form with a set of admissible weighting functions and integration over the domain. We seek to find f u (x, t ) ∈ L u͠ , p f (x, t ) ∈ L p , and θ (x, t ) ∈ L θ with t ∈ [0, T ] such that:

t

σ ′ = ∫−∞ G vol (t − ξ ) ε˙ s,vol (ξ ) dξ t

+ ∫−∞ G dev (t − ξ ) ε˙ s,dev (ξ ) dξ

⎧ R u: ∫Ω ∇wu·σdv − ∫Ω ρ wu·bdv ⎪ ∈Ω ⎪− ∫Γt wu σ ·nds = 0 ⎪ f ∼ s ⎪ R p f : ∫Ω w p f ∇ ·v dv − ∫Ω ∇w p f ·v D dv ⎪ D sθ ⎪− ∫Ω (3nsαθs + 3nf αθf ) w p f dv Dt ⎪ f f f ∼ f ⎪− ∫ 3α w f ∇θ ·v dv + ∫ w f ∼ ∈Ω θ p D p v D · nds = 0 Ω Γs ⎪ s ⎪ ⎪ R : ∫ w (ρ s C s + ρ f C f ) D θ dv θ θ p p , Ω Dt ⎨ sp f ⎪ D f f s s dv ⎪− ∫Ω wθ (3n αθ θ + 3n αθ θ) Dt ⎪ s ⎪+ s D σ′ ∼f fR f ⎪ ∫Ω wθ αθ θ Dt : Idv + ∫Ω wθ ρ Cp ∇θ ·v D dv ⎪ f f ∼f ⎪− ∫Ω wθ 3αθ θ ∇p ·v D dv − ∫Ω ∇wθ·qdv ⎪ f f s s ⎪− ∫Ω wθ (ρ r + ρ r ) dv ⎪ f ∼ f ⎪+ ∫Ω wθ v D (∇p − ρ fR g) dv + ∫Γθ wθ q·nds = 0 ∈ Ω ⎩

t

− ∫−∞ G thm (t − ξ ) ε˙ s,thm (ξ ) dξ .

(29)

The relaxation modulus function, G (t ) , that considers viscoelastic behavior of material across a range of time can be expressed as Simo and Hughes (2006):

⎧ thm Nv k ∞ ⎛ ⎛ t ⎞⎞ vol ⎪G (t ) = G (t ) = 3 K skel ⎜1 + ∑m = 1 ωm exp ⎜− k ⎟ ⎟ ⎪ ⎝ τm ⎠ ⎠ ⎝ , ⎨ Nv μ ∞ ⎛ ⎛ t ⎞⎞ ⎪G dev (t ) = 2μskel ⎜1 + ∑m = 1 ωm exp ⎜ − μ ⎟ ⎟ ⎪ ⎝ τm ⎠ ⎠ ⎝ ⎩

(30)

where Nv is the number of spring-dashpot devices. For the mth springdashpot device, ωm and τm are defined as follows:

ωmk =

eK m skel , ∞ K skel

ωmμ =

eμ m skel ∞ μskel

τmk =

vK m skel , eK m skel

τmμ =

vμ m skel , eμ m skel

(31)

Nv

∑ m=1

Nv

+

∑ m=1 Nv

where L u , L , L θ are the trial solution spaces defined as follows:

−

∑ m=1 Nv

+

∑ m=1

L u = {u: Ω × [0, T ] ↦ n, u ∈ H1, u (t ) = ∼ u (t ) on Γu, u (x, 0) = u 0 (x)}, f ∼f L p = {p f : Ω × [0, T ] ↦ n, p f ∈ H1, p f (t ) = p (t ) on Γs, p f (x, 0) = p0f (x)},

t n+1 ωmk exp ⎜⎛− k ⎟⎞ σ ′∞ (vol) (0) ⎝ τm ⎠

L θ = {θ : Ω × [0, T ] ↦ n, θ ∈ H1, θ (t ) = θ͠ (t ) on Γθ, θ (x, 0) = θ0 (x)},

t n+1 ωmμ exp ⎜⎛− μ ⎟⎞ σ ′∞ (dev) (0) ⎝ τm ⎠ ωmk

t exp ⎜⎛− ⎝

n+1

Nv

Nv

∑ m=1

hm (dev) (t n + 1) −

∑

(35)

and wu , w p f , and wθ are weighting functions of the independent state variables. H1 is the first Sobolev space (Hughes, 2012). The symbol Ω is the volume occupied by the mixture, Γt is the applied traction boundary, Γu is the prescribed displacement boundary, Γs is the prescribed pressure boundary, and Γθ is the prescribed temperature boundary. The prescribed and initial conditions for the state variables are denoted by ∼f u(t ) , p (t ) , θ͠ (t ) ] and [u 0 (x) , p0f (x) , θ0 (x) ], respectively. [∼ The domain, Ω is discretized into the non-overlapping sub-domains Ωe , 1 ≤ e ≤ nel , where nel is the number of elements in the domain. To avoid oscillation in mechanical stresses, pressure, and temperature, the solutions are approximated by using mixed elements such that a quadratic interpolation is used for the displacement field and linear interpolation is used

⎞ ∞ (thm) (0) ⎟ σ′

τmk ⎠

hm (vol) (t n + 1) +

(34)

pf

m m m with v μskel and v Kskel being the viscous shear and bulk moduli and e μskel e m and Kskel being the elastic shear and bulk moduli of the mth viscoelastic device, respectively. Assume all stresses before t = 0 are zero, then at time t n + 1 ≥ 0 , the effective stress of the solid skeleton is expressed as:

σ ′ (t n + 1) = σ ′∞ (t n + 1) +

(33)

hm (thm) (t n + 1),

m=1

(32) where σ ′∞ is the elastic stress defined in the Eq. (28) and hm (cmp) (t n + 1) is stress like Internal State Variables (ISVs) defined as Taylor et al. (1970):

217

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

Fig. 4. Normalized pressure at the bottom of the model.

Fig. 2. Schematic representation of the generalized Maxwell model.

Table 1 Material properties and model parameters for the laterally confined 1D column. Description Model height Fluid volume fraction Hydraulic permeability Young's modulus of solid skeleton Shear modulus of solid skeleton Applied traction Initial pressure

Symbol

Value

Unit

Hc

0.3 0.42

m –

1.0 × 10−14 101.8 47 10 10

m2/(Pa·s) MPa MPa MPa MPa

50 3696

°C

nf kp K skel μskel Fx p0f

Initial temperature Real mass density of solid

θ0

Real mass density of fluid

1000

kg/m3

Thermal conductivity of solid

ρ fR κθs

1.38

W/(m. K)

Thermal conductivity of fluid

κθf

0.6

W/(m. K)

Specific heat capacity of solid

Cps

703

J/(kg. K)

Specific heat capacity of fluid

Cpf

4.18 × 103

J/(kg. K)

Thermal expansion coefficient of solid

αθs

1.65 × 10−6

1/K

Thermal expansion coefficient of fluid

αθf

2.07 × 10−4

1/K

ρ sR

kg/m3

Fig. 5. Normalized temperature at the bottom of the model.

depends on the spatial and temporal discretization step sizes and mechanical parameters of the model. For more information, the reader is referred to Favino et al. (2013). 3.1. Time integration The weak form of the governing equations is discretized in time by an implicit backward Euler scheme. The discretized governing equations are written in the following compact form:

1 n + 1 ^n + 1 ^n n+1 − s ) + Rn + 1(^s ) = 0, M (s Δt 1 n+1 n+1 +1 = + Jn + 1 (^s ), Jndyn M Δt

+1 = R ndyn

(36)

+1 R ndyn

where is the dynamic residual of all governing equations given in Eq. (34) and written in a compact form as follows: +1 R ndyn = [R u, R p f , R θ]T .

(37)

Mn + 1 is the capacitance matrix that collects the contributions of the n+1 n inertia terms in Eq. (34). ^s and ^s are the vector of all state variables at current and previous time steps, respectively. Rn + 1 and Jn + 1 are the residual and Jacobian of the static contribution to the dynamic residual, +1 is the analytical derivatives of dynamic residual equations, and Δt Jndyn is the size of time step.

Fig. 3. Normalized displacement at the top of the model.

for approximating the pressure and temperature fields. For poro- elastic and viscoelastic problems, numerical oscillations may occur if the spatial and the temporal discretization parameters are not properly chosen (Murad and Loula, 1992; Wan, 2003; Preisig and Prévost, 2011). To ensure the stability of the system, the stability condition is computed through the relationship given by Favino et al. (2013). This stability condition

4. Results and discussion In this section, several examples are presented and discussed to illustrate the efficiency of the numerical model. For this purpose, the accuracy of the finite element method is verified through comparison 218

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

Fig. 6. Schematic of confined and unconfined compression models.

Table 2 Material properties and model parameters for the articular cartilage model. Description

Symbol

Value

Unit

Ref.

Specimen wide

a

0.25 × 10−3

m

Hoang and Abousleiman (2009)

Specimen thickness

b

m

Hoang and Abousleiman (2009)

Fluid volume fraction

nf kp K skel μskel Nv e m K skel

0.50 × 10−3 0.8 1.0 × 10−15 0.2 0.1 5 [89, 36, 36, 8.9, 2.9]

m2/(Pa·s) MPa MPa – KPa

Hydraulic permeability Bulk modulus of solid skeleton Shear modulus of solid skeleton Number of spring-dashpot devices Elastic bulk moduli of devices

– Armstrong et al. (1984) Loret and Simoes (2017) Loret and Simoes (2017)

Elastic shear moduli of devices

e

m μskel

[97, 27, 16, 9.7, 1.7]

KPa

Viscous bulk moduli of devices

v

m K skel

[0.089, 0.036, 3.6, 8.9, 29]

KPa

Viscous shear moduli of devices

v

m μskel

[0.097, 0.27, 1.6, 9.7, 17]

KPa

Applied traction

fy

5.0 × 10 4

Pa

Real mass density of solid

ρ sR

1260

kg/m3

Loret and Simoes (2017)

Real mass density of fluid

997.1

kg/m3

Loret and Simoes (2017)

Thermal conductivity of solid

ρ fR κθs

0.6

W/(m. K)

Loret and Simoes (2017)

Thermal conductivity of fluid

κθf

0.6

W/(m. K)

Loret and Simoes (2017)

Specific heat capacity of solid

Cps

3.17 × 103

J/(kg. K)

Loret and Simoes (2017)

Specific heat capacity of fluid

Cpf

4.18 × 103

J/(kg. K)

Loret and Simoes (2017)

Thermal expansion coefficient of solid

αθs

0.01 × 10−3

1/K

Loret and Simoes (2017)

Thermal expansion coefficient of fluid

αθf @24 °C

0.246 × 10−3

1/K

Thermal expansion coefficient of fluid

αθf @37.5 °C

0.376 × 10−3

1/K

Thermal expansion coefficient of fluid

αθf @55 °C

0.544 × 10−3

1/K

Fig. 7. Variations of displacement and pore fluid pressure along the side of the confined model (dashed line in Fig. 6(a)) and for different loading time. The solid lines represent the poroelastic responses without thermal coupling and the dashed lines show the thermo-poro-elastic responses. The origin of the normalized length is represented with a blue circle in Fig. 6(a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

219

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

Fig. 8. Variations of displacement and pore fluid pressure along the side of the confined model (dashed line in Fig. 6(a)) and for different loading time. The solid lines represent the poroviscoelastic responses without thermal coupling and the dashed lines show the thermo-poro-viscoelastic responses. The origin of the normalized length is represented with a blue circle in Fig. 6(a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9. Evolution of y-displacement and pore fluid pressure at the top and bottom midpoints of the confined model and for different temperatures.

extended to verify the thermo-poro-elastic response of the model. For this purpose, a laterally confined 1D column with the height Hc proposed by Bai and Abousleiman (1997) is considered. A load with magnitude of Fx is applied at the top of the column (i.e. at x = 0 ) while the bottom of the column is rigid. At the top of the column, the pressure and temperature are set to zero, which allows the fluid and the temperature to dissipate from the top. The geometric configuration and material parameters are given in Table 1. For numerical modeling, the column is discretized with 100 quadratic elements, and the simulation is continued until the pressure and temperature dissipate from the entire column. The variations of displacement at the top and pressure and temperature at the bottom of the column are presented as a function of consolidation time as shown in Figs. 3–5. The displacement, pressure, and temperature are normalized with respect to their maximum values. The results are compared against analytical solutions (Bai and Abousleiman, 1997). The comparison demonstrates good agreement between numerical model and analytical solutions.

against analytical reference solutions. To gain insight into the effects of different parameters on the performance of articular cartilage, confined and unconfined models are considered. In particular, the effects of material type, temperature, and boundary conditions are investigated. To determine the physical response appropriate for 3D domain, the plain strain assumption is used in 2D models. For all numerical examples, mesh refinement studies are performed and only results for sufficiently fine meshes with negligible discretization errors are shown. Newton-Raphson method is used to solve nonlinear equations using analytically derived Jacobians. The discretized linear problem is considered converged if the relative residuals are less than 10−9 . 4.1. Verification example The accuracy of the implemented model discussed previously is verified through comparison against analytical reference solutions. Originally verified for a poro-elastic model (Behrou, 2017), this study is 220

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

Fig. 10. Evolution of y-displacement and pore fluid pressure at the top and bottom midpoints of the unconfined model and for different temperatures.

recovered. For the poro-elastic behavior, models with thermal effects show a significant reduction in both pore fluid pressure and displacement at early times. However, for the poro-viscoelastic behavior, this reduction is not considerable due to creep recovery and constant stress over the loading period. To explore the effect of temperature on the behavior of articular cartilage, different temperatures (24, 37.5, and 55 °C) are considered for both confined and unconfined models with thermo-poro-viscoelastic behavior. For the unconfined model, as shown in Fig. 6 (b), the fluid pore pressure can dissipate along the sides, and the temperature can dissipate from the top of the model. At the bottom, the displacements in both x and y directions are constrained to zero. Impermeable solid plates at the top and bottom of the model are used. The top plate uniformly transfers the applied load to the articular cartilage specimen. For both confined and unconfined models, the evolution of y-displacement at the top midpoint and pore fluid pressure at the bottom midpoint are presented for different temperatures, as shown in Figs. 9 and 10, respectively. The results reveal the influence of temperature on the articular cartilage behavior. The results show a significant influence of temperature at very early times. After dissipation of pore fluid pressure and temperature, elastic response is recovered in both confined and unconfined models.

Fig. 11. Confined compression articular cartilage model with partial loading. Dissipation of pore fluid pressure and temperature is not allowed through the solid blue line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.2. Confined and unconfined articular cartilage models To explore the characteristics of the articular cartilage model, both confined and unconfined boundary conditions are considered, as shown in Fig. 6. For the confined model (Fig. 6 (a)), articular cartilage is compressed into a perfectly smooth and impermeable box. The fluid pore pressure and temperature can dissipate from the top boundary. The model is fixed at the bottom and displacements in the x direction are constrained to zero along the sides. A porous plate at the top of the model is used to uniformly transfer the applied load, f y , to the articular cartilage specimen. For numerical modeling, the domain is discretized with 55 × 110 elements; the geometric configuration and material parameters, used in the test, are shown in Table 2. The behavior of the articular cartilage model is studied considering poro- elastic and viscoelastic material behaviors. The results are compared with those that consider the influence of temperature on both poro- elastic and viscoelastic behavior through thermal coupling. To this end, the confined model is analyzed at θ = 37.5 °C . The results are presented in terms of variations of displacement and pore fluid pressure along the side of the confined model and for different loading times. The final loading time, t f , that allows full dissipation of pore fluid pressure and temperature in the model is set to 4 × 103 seconds. For poro- elastic and viscoelastic materials, the displacement and pore fluid pressure profiles are shown in Figs. 7 and 8, respectively, where solid lines represent responses without thermal coupling and dashed lines show those with thermal coupling. The results of poro- elastic and viscoelastic models show that at very early times, large fluid pore pressure results in small displacements. As the pore fluid pressure dissipates, the deformation increases until a purely elastic response is

Fig. 12. Evolution of displacement along the top of the model and for different temperatures and loading stages.

221

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

Fig. 11. All model and material parameters are adopted from Table 2, except for the width of the specimen which is set to a = 5.0 × 10−3 m . The fluid pore pressure and temperature can dissipate only from the bottom of the porous plate. The displacement in the x direction are constrained to zero along the sides, and all displacements are fixed at the bottom. Again, a porous plate at the top of the model is used to uniformly transfer the applied load, f y , to the articular cartilage specimen. The domain is discretized with 125 × 25 quadratic elements and the final loading time, t f , that allows full dissipation of pore fluid pressure and temperature is set to 2 × 105 seconds. The viscoelastic behavior is studied at different temperatures and the evolution of displacement along the top and pore fluid pressure along the bottom of the model are given at different loading stages (Fig. 12 and Fig. 13). Initially, high pore fluid pressure is established in the half of the model where the pressure can dissipate. However, as the simulation proceeds the vertical displacement increases and the pore fluid pressure dissipates from the bottom of the porous plate. The results also show that at very early times, temperature has a significant influence on the articular cartilage behavior. High temperature results in low pressures and large vertical displacements. Furthermore, the results show fast dissipation of temperature at early times. Again, once the temperature dissipates, the linear elastic response is recovered. The inset snapshots given in Figs. 14 and 15 at the early loading stages depict a significant pressure drop and large displacement in the half of the model. Due to free displacement boundary conditions on the right-half, expansion occurs at the top middle of the model.

Fig. 13. Evolution of pore fluid pressure along the bottom of the model and for different temperatures and loading stages.

4.3. Partially loaded confined articular cartilage The effects of partial loading are an extremely important facet of articular cartilage mechanics on account of the fact that medical injuries in articular cartilage are mostly due to partial loads. The following example extends the confined model (Fig. 6 (a)) to study the behavior of articular cartilage under partial loading condition shown in

Fig. 14. Displacement distribution in the articular cartilage specimen at θ = 37.5 °C and for different loading stages.

Fig. 15. Pore fluid pressure distribution in the articular cartilage specimen at θ = 37.5 °C and for different loading stages.

222

Journal of the Mechanical Behavior of Biomedical Materials 78 (2018) 214–223

R. Behrou et al.

The model introduced in this paper can have applications for studies in ablating diseased cartilage and bone cavity treatment (Lu et al., 2000; Preininger et al., 2012; Alambeigi et al., 2017). For instance, Radiofrequency (RF) energy has gained an extensive attention for articular cartilage shrinking and cutting. The RF treatment smooths and stabilizes the articular surface in a partial-thickness cartilaginous defect, and consequently, interrupts the cascade that could lead to osteoarthritis and joint dysfunction. In RF treatment, an alternating current is used between the application probe and the grounding plate, which produces molecular friction in tissues resulting in tissue heating. As another example, Poly-Methylmethacrylate (PMMA) is a common cavity filling material for bone cavity treatment following tumor excision. During this treatment, a great degree of heat is set free which affects articular cartilage and its behavior.

Lett. 2 (3), 1480–1487. Armstrong, C., Lai, W., Mow, V., 1984. An analysis of the unconfined compression of articular cartilage. J. Biomech. Eng. 106 (2), 165–173. Bai, M., Abousleiman, Y., 1997. Thermoporoelastic coupling with application to consolidation. Int. J. Numer. Anal. Methods Geomech. 21 (2), 121–132. Behrou, R., 2017. A thermo-poro-viscoelastic model for the behavior of polymeric porous separators. ECS Trans. 80 (8), 583–591. Boer, R.De, 2006. Trends in Continuum Mechanics of Porous Media. 18 Springer Science & Business Media. Coleman, B.D., Noll, W., 1963. The thermodynamics of elastic materials with heat conduction and viscosity. Arch. Ration. Mech. Anal. 13 (1), 167–178. Coussy, O., 2004. Poromechanics. John Wiley & Sons. Favino, M., Grillo, A., Krause, R., 2013. A stability condition for the numerical simulation of poroelastic systems. In: Poromechanics V: Proceedings of the Fifth Biot Conference on Poromechanics, pp. 919–928. Gawin, D., Schrefler, B.A., Galindo, M., 1996. Thermo-hydro-mechanical analysis of partially saturated porous materials. Eng. Comput. 13 (7), 113–143. Gu, K., Li, L., 2011. A human knee joint model considering fluid pressure and fiber orientation in cartilages and menisci. Med. Eng. Phys. 33 (4), 497–503. Gupta, S., Lin, J., Ashby, P., Pruitt, L., 2009. A fiber reinforced poroelastic model of nanoindentation of porcine costal cartilage: a combined experimental and finite element approach. J. Mech. Behav. Biomed. Mater. 2 (4), 326–338. Hoang, S.K., Abousleiman, Y.N., 2009. Poroviscoelastic two-dimensional anisotropic solution with application to articular cartilage testing. J. Eng. Mech. 135 (5), 367–374. Holzapfel, G.A., 2000. Nonlinear Solid Mechanics. 24 Wiley, Chichester. Huang, C.-Y., Soltz, M.A., Kopacz, M., Mow, V.C., Ateshian, G.A., et al., 2003. Experimental verification of the roles of intrinsic matrix viscoelasticity and tensioncompression nonlinearity in the biphasic response of cartilage. Trans.-Am. Soc. Mech. Eng. J. Biomech. Eng. 125 (1), 84–93. Hughes, T.J., 2012. The finite element method: linear static and dynamic finite element analysis. Cour. Corp. June, R.K., Fyhrie, D.P., 2010. Temperature effects in articular cartilage biomechanics. J. Exp. Biol. 213 (22), 3934–3940. Lewis, R., Schrefler, B., 1999. The finite element method in the static and dynamic deformation and consolidation of porous media. Meccanica 34 (3), 231–232. Loret, B., Simoes, F.M.F., 2017. Biomechanical Aspects of Soft Tissues. CRC Press. Lu, Y., Hayashi, K., Hecht, P., Fanton, G.S., Thabit, G., Cooley, A., Edwards, R.B., Markel, M.D., 2000. The effect of monopolar radiofrequency energy on partial-thickness defects of articular cartilage. Arthrosc.: J. Arthrosc. Relat. Surg. 16 (5), 527–536. Mak, A., 1986. The apparent viscoelastic behavior of articular cartilagethe contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J. Biomech. Eng. 108 (2), 123–130. Mow, V.C., Ateshian, G.A., Spilker, R.L., 1993. Biomechanics of diarthrodial joints: a review of twenty years of progress. J. Biomech. Eng. 115 (4B), 460–467. Murad, M.A., Loula, A.F., 1992. Improved accuracy in finite element analysis of biot's consolidation problem. Comput. Methods Appl. Mech. Eng. 95 (3), 359–382. Netti, P., Travascio, F., Jain, R., 2003. Coupled macromolecular transport and gel mechanics: poroviscoelastic approach. AIChE J. 49 (6), 1580–1596. Pena, E., Calvo, B., Martinez, M., Palanca, D., Doblaré, M., 2005. Finite element analysis of the effect of meniscal tears and meniscectomies on human knee biomechanics. Clin. Biomech. 20 (5), 498–507. Preininger, B., Matziolis, G., Pfitzner, T., Hardt, S., Perka, C., Röhner, E., 2012. In situ tele-thermographic measurements during pmma spacer augmentation in temporary arthrodesis after periprosthetic knee joint infection. Technol. Health Care 20 (4), 337–341. Preisig, M., Prévost, J.H., 2011. Stabilization procedures in coupled poromechanics problems: A critical assessment. Int. J. Numer. Anal. Methods Geomech. 35 (11), 1207–1225. Shirazi, R., Shirazi-Adl, A., Hurtig, M., 2008. Role of cartilage collagen fibrils networks in knee joint biomechanics under compression. J. Biomech. 41 (16), 3340–3348. Simo, J.C., Hughes, T.J., 2006. Computational Inelasticity. 7 Springer Science & Business Media. Sophia Fox, A.J., Bedi, A., Rodeo, S.A., 2009. The basic science of articular cartilage: structure, composition, and function. Sports Health 1 (6), 461–468. Taylor, R.L., Pister, K.S., Goudreau, G.L., 1970. Thermomechanical analysis of viscoelastic solids. Int. J. Numer. Methods Eng. 2 (1), 45–59. Tomic, A., Grillo, A., Federico, S., 2014. Poroelastic materials reinforced by statistically oriented fibresnumerical implementation and application to articular cartilage. IMA J. Appl. Math. 79 (5), 1027–1059. Wan, J., 2003. Stabilized Finite Element Methods for Coupled Geomechanics and Multiphase Flow. (Ph.D. Thesis), Stanford University.

5. Conclusions The presented finite element framework proved reliable for a variety of 2D articular cartilage models. The accuracy of the implemented framework is verified through comparison against analytical reference solutions. To gain insight into the characteristics of articular cartilage, both confined and unconfined models are considered. The poro-mechanical behavior of articular cartilage is investigated for both elastic and viscoelastic materials, with and without thermal coupling. The influence of temperature is studied through examining the behavior of articular cartilage for different temperatures. Partially loaded confined articular cartilage model is considered to investigate deformation, pore fluid pressure and temperature dissipation processes. The results show a significant influence of temperature on the behavior of articular cartilage at early loading stages. As the temperature and pore fluid dissipate from the system, the steady-state elastic behavior is recovered. The model presented in this paper can have applications for medical treatments in which the temperature of articular cartilage is changed beyond the normal body temperature, e.g. ablating diseased cartilage and bone cavity treatment. While only 2D numerical models are considered in this study, the developed method allows for the convenient extension of the framework to 3D behavior of articular cartilage. For future studies, this method needs to be extended to explore the 3D characteristics of articular cartilage in the finite strain regime. Conflict of interest statement The authors declare that there is no conflict of interest. Funding sources This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. References Alambeigi, F., Wang, Y., Sefati, S., Gao, C., Murphy, R.J., Iordachita, I., Taylor, R.H., Khanuja, H., Armand, M., 2017. A curved-drilling approach in core decompression of the femoral head osteonecrosis using a continuum manipulator. IEEE Robot. Autom.

223