**4. Mathematical model for articular cartilage**

Several authors (Mow, 1980 ; Haider & Schugart, 2006; Wilson et al., 2005; Haider & Guilak, 2007; Meng et al., 2002; Wu et al., 1997; Terada et al., 1998; Donzelli et al., 1999; Donzelli & Spilker, 1998); have conducted their research on cartilage from biphasic behavior that this tissue exhibits. This allowed analysis of the same material as a poro-elastic behavior capable of supporting loads. The mathematical model of the AC as a biphasic material, analyzes the displacement *u (t, x)* of the solid (matrix) and pressure *p (x)* of the fluid displaced by the load, thanks to its characteristic of poro-elasticity. This model is described by the equations (1) and (2):

$$-\nabla.(2\mu\_{s}\varepsilon(\mathbf{u}) + \lambda\_{s}\nabla.\mathbf{u}\,\mathbf{I}) + \nabla\mathbf{p} = 0 \text{ en } \Omega\tag{1}$$

$$\frac{\partial}{\partial t}(\nabla \cdot \mathbf{u}) - \nabla . (\mathbf{k} \nabla \mathbf{p}) = 0 \quad \text{en } \Omega \tag{2}$$

Equation (1) is derived from the law of conservation of momentum and corresponds to the linear elasticity equation (term 1a) coupled with a term that represents fluid pressure (term 1b). The term (u) corresponds to the strain tensor acting on the surface enclosed by. s s y are the "Lame elastic constants" for the solid, related to Young's modulus and Poisson's ratio (E, ν). For its part, the equation (2) refers to the change of the dilation of the solid matrix (term 2a) due to the mechanical load created by the divergence of the gradient of the pressure of fluid contained in the domain (term 2b) (Frijns, 2000). In this equation, k is a constant representing the permeability of solid module.

### **4.1. Boundary conditions**

204 Injury and Skeletal Biomechanics

(1) and (2):

1b). The term (u)

element method (FEM) (Garzón, 2007).

**3. Mechanical behavior of articular cartilage (AC):** 

**4. Mathematical model for articular cartilage** 

matrix strong function of time (flow-independent mechanism) (Wilson, 2005).

substances. Both biological and mechanical parts are combined in a computational model, which considers the application of forces, mechano-transduction, cellular expression, genetics and the transformation of the characteristics of the extracellular matrix. The typical method of numerical implementation of these mechano-biological problems is the finite

Finite element computational analysis has been used as an approach to diverse biological processes including the biomechanical behavior of articular cartilage (Wilson et al., 2004; Ateshian et al., 1997; Chan et al., 2004; Wilson, 2005; Donzelli et al, 1999; Donzelli & Spilker, 1998; Almeida & Spilker, 1998; Wu & Herzog, 2000; Levenston et al., 1998). Using the material representation of the continuum phase of cartilage, results have indicated that local intermittent hydrostatic pressure promotes cartilage maintenance (Carter & Wong, 2003).

Mechanical properties of articular cartilage (AC) are attributed to its complex structure and to composition of its ECM including a fluid phase (water with dissolved ions), and a solid matrix (collagen type II, aggregates of PGs, proteins, lipids, and cells) (Haider & Guilak, 2007). With the mechanical load, the interstitial fluid is redistributed through the pores of the permeable solid matrix, resulting in predominantly poro-elastic conduct. This behavior of the AC is mainly due to two mechanisms: (a) the frictional force due to drag flow of interstitial fluid through the porous solid matrix (flow-dependent mechanism), and (b) the deformability of the

Several authors (Mow, 1980 ; Haider & Schugart, 2006; Wilson et al., 2005; Haider & Guilak, 2007; Meng et al., 2002; Wu et al., 1997; Terada et al., 1998; Donzelli et al., 1999; Donzelli & Spilker, 1998); have conducted their research on cartilage from biphasic behavior that this tissue exhibits. This allowed analysis of the same material as a poro-elastic behavior capable of supporting loads. The mathematical model of the AC as a biphasic material, analyzes the displacement *u (t, x)* of the solid (matrix) and pressure *p (x)* of the fluid displaced by the load, thanks to its characteristic of poro-elasticity. This model is described by the equations

> .(2 (u) .u I ) p 0 en s s

> > ( .u) .(k p) 0 en

Equation (1) is derived from the law of conservation of momentum and corresponds to the linear elasticity equation (term 1a) coupled with a term that represents fluid pressure (term

s s y are the "Lame elastic constants" for the solid, related to Young's modulus and

corresponds to the strain tensor acting on the surface enclosed by.

(2)

t

(1)

Boundary conditions of the model are defined in domain and may be dependent on time. The mathematical expression of these conditions is (3-6):

$$\mathbf{u} = \mathbf{g}\_u^{\mathrm{D}} \quad \text{en } \Gamma\_u^{\mathrm{D}} \tag{3}$$

$$\operatorname{Im} \left( 2\mu\_s \mathop{\rm c} \nolimits \nolimits\_\dash \right) + \lambda\_s \nabla \mathop{\rm u} \mathop{\rm l} \nolimits \mathop{\rm l} - \mathop{\rm p} \mathop{\rm l} \nolimits = \mathop{\rm g}\_u^N \mathop{\rm en} \mathop{\rm l}^N\_u \tag{4}$$

$$\mathbf{p} = \mathbf{g}\_p^{\mathcal{D}} \text{ en } \Gamma\_p^{\mathcal{D}} \tag{5}$$

$$-\mathbf{n}.\text{(k}\nabla\mathbf{p}) = \mathbf{g}\_p^N \text{ en } \Gamma\_p^N\tag{6}$$

$$\text{with } \Gamma = \Gamma\_u^P + \Gamma\_u^N + \Gamma\_P^D + \Gamma\_p^N$$

A widely used method for solving partial differential equations in complex geometries is the finite element method (Garzón, 2007). This method allows implementing the numerical model presented in equations (1) and (2) simply and with low computation cost. The method consists of using a vectorial function *W* or *weighting function* and a scalar function of *q*, which minimizes the terms of the equations (1) and (2). Multiplying (1) by *W* and (2) by *q*, and performing integration by parts in the domain, we obtain a variational of the form (Frijns, 2000):

$$\mathbf{a(u,W) + b(W,p) = (g\_u^N, W)}\tag{7}$$

$$\frac{\partial}{\partial t}\mathbf{b}(\mathbf{u},\mathbf{q}) - \mathbf{c}(\mathbf{p},\mathbf{q}) = (\mathbf{g}\_{\mathbf{p}}^{\mathcal{N}}, \mathbf{q}) \tag{8}$$

where:

$$\mathbf{a(u,W)} = \int\_{\Omega} (\mathfrak{A}\mu\_s \mathop{\mathbf{e(u)}}\_{\mathop{\mathbf{e(u)}}} \mathbf{e(W)} + \lambda\_s (\nabla.\mathbf{u})(\nabla.W) \mathop{\mathbf{V}}\_{\mathop{\mathbf{u}}} \mathbf{W}) \tag{9}$$

$$\mathbf{b}(\mathbf{W}, \mathbf{q}) = -\int\_{\Omega} \nabla. \mathbf{W} \mathbf{q} \tag{10}$$

$$\mathbf{c}(\mathbf{p}, \mathbf{q}) = \int\_{\Omega} (\mathbf{k} \nabla \mathbf{p}) \nabla \mathbf{q} \tag{11}$$

$$\mathbf{f}\left(\mathbf{g}\_{u}^{\mathcal{N}},\mathcal{W}\right) = \int\_{\Gamma\_{u}^{\mathcal{N}}} \mathbf{g}\_{u}^{\mathcal{N}} \mathcal{W} \tag{12}$$

#### 206 Injury and Skeletal Biomechanics

The solution method is based on dividing the domain (continuum) in which are defined equations (1) and (2) in their integral form in a series of non-intersecting subdomains <sup>e</sup> called *finite elements*. The set of finite elements forms a partition of the domain called *discretization* which together are described by equations (14) and (15):

$$\mathbf{a}(\mathbf{u}\_{\mathbf{u}'}^{\hbar}, \mathbf{W}\_{\mathbf{h}}^{\hbar}) + \mathbf{b}(\mathbf{W}\_{\mathbf{h}}^{\hbar}, \mathbf{p}\_{\mathbf{n}}^{\hbar}) = (\mathbf{g}\_{\mathbf{u}'}^{\hbar}, \mathbf{W}\_{\mathbf{h}}^{\hbar}) \tag{14}$$

Mechanical Behavior of Articular Cartilage 207

km f en e (18)

KM F en (19)

**Figure 5.** Representation of the *shape functions standard*. a) Case 1D, 3 nodes per element. b) Case 2D

Solution of these equations for both *u* and *p* are implemented using a routine programmed

Conditions for the simulation were based on experiments by Frijns, 2000 and Wu et al., 1999. Simulations were performed for 1D and 2D looking tissue response to the application of: (a) compression, (b) tensile and (c) oscillating or cyclical loads. The simulation time is 45 seconds of load in all cases. The loads applied in 1D and 2D simulations were performed with the parameters shown in Table 1. The tests considered the AC as a continuous and

> E K Matrix 0.40 Mpa 0.1 1.0 x10-2mmN-1s-1 Chondrocyte 1.0 kPa 0.2 1.0 x10-5mmN-1s-1

For the 1D simulation a mesh was performed that represents a fragment of 0.3 mm AC. From this mesh 641 nodes and 320 elements of 3 nodes were obtained. In the 2D simulation a mesh was performed that represents a fragment of articular cartilage 0.09 mm2 (0.3 mm x 0.3 mm). In this case 10,201 nodes were obtained, equivalent to 10,000 elements of 4

in Fortran and a desktop PC with 2.4 GHz AMD processor and 1.0 GB of RAM

homogeneous material with the chondrocytes being part of the continuum.

**Table 1.** Parameters for the homogeneous material characteristics (Wu et al., 1999).

elements with four nodes (Radi, 2007).

**4.2. Model implementation** 

**4.3. Computer simulation** 

*4.3.1. Meshing* 

nodes.

$$\mathbf{b}(\mathbf{u}\_{\text{n}}^{\text{h}}, \mathbf{q}^{\text{h}}) - \Delta \mathbf{t} \mathbf{c}(\mathbf{p}^{\text{h}}, \mathbf{q}\_{\text{n}}^{\text{h}}) = \mathbf{b}(\mathbf{u}\_{\text{n-1}}, \mathbf{q}^{\text{h}}) + \Delta \mathbf{t} (\mathbf{g}\_{\text{b}}^{\text{N}}, \mathbf{q}^{\text{h}}) \tag{15}$$

In (15) the n 1 subscript indicates the value of a parameter at a t 1 time, while the subscript n indicates the value in the next time, n n1 tt t (Frijns, 2000). Using the Galerkin method, the functions *W* and *q* are approximated by the matrix expression (16) where the *N* vector contains *shape functions standard (,)* which allow interpolation within each element of the domain.

$$\mathbf{q}^{\text{h}} = \mathbf{W}\_{\text{h}}^{\text{h}} = \begin{pmatrix} \text{N}\_{1} \\\\ \text{N}\_{2} \\\\ \text{N}\_{3} \\\\ \text{N}\_{3} \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \\ (1 + \xi)(1 - \xi) \\\\ 1 + \xi(1 + \eta) \\\\ \frac{1}{2}\xi(\xi + 1) \end{pmatrix} \tag{16}$$

$$\mathbf{W}\_{\text{h}}^{\text{h}} = \begin{pmatrix} \text{N}\_{1} \\\\ \text{N}\_{2} \\\\ \text{N}\_{3} \\\\ \text{N}\_{3} \\\\ \text{N}\_{4} \end{pmatrix} = \begin{pmatrix} 1 \\ 4 \\\\ \frac{1}{4}(1 + \xi)(1 - \eta) \\\\ \frac{1}{4}(1 + \xi)(1 + \eta) \\\\ \frac{1}{4}(1 - \xi)(1 + \eta) \end{pmatrix} \tag{17}$$

In (16), the shape functions correspond to the case of one-dimensional element with three nodes (See Fig. 5a). For two-dimensional case four node elements are used, whose standard shape functions are shown in matrix form (17) and in Figure 5b.

If the integration space *(x,y)* is changed to *(, )* by the *Jacobian* of the *transformation* and vector notation is adjusted, the equations (1) and (2) can be reduced to a matrix system of elementary type of the form (18) which corresponds to the algebraic discretization in e domain of an element, with *k* elementary stiffness matrix, *m* the unknowns and *f* the independent term. Joining the result of (18) for total elements in , a general matrix system is obtained defined as (19) where *K* is the global matrix of stiffness, *M* is the vector of unknowns and *F* is the global vector entries.

Mechanical Behavior of Articular Cartilage 207

$$\mathbf{k} \cdot \mathbf{m} = \mathbf{f} \text{ en } \Omega^{\mathbf{e}} \tag{18}$$

$$\mathbf{K} \cdot \mathbf{M} = \mathbf{F} \text{ en } \Omega \tag{19}$$

**Figure 5.** Representation of the *shape functions standard*. a) Case 1D, 3 nodes per element. b) Case 2D elements with four nodes (Radi, 2007).

### **4.2. Model implementation**

206 Injury and Skeletal Biomechanics

each element of the domain.

The solution method is based on dividing the domain (continuum) in which are defined equations (1) and (2) in their integral form in a series of non-intersecting subdomains <sup>e</sup> called *finite elements*. The set of finite elements forms a partition of the domain called

h h hh N h

hh hh h Nh

In (15) the n 1 subscript indicates the value of a parameter at a t 1 time, while the subscript n indicates the value in the next time, n n1 tt t (Frijns, 2000). Using the Galerkin method, the functions *W* and *q* are approximated by the matrix expression (16)

1

qW N 1 1

<sup>N</sup> <sup>1</sup> ( 1) <sup>2</sup>

<sup>1</sup> ( 1) <sup>N</sup> <sup>2</sup>

<sup>N</sup> <sup>1</sup> 1 1 4

<sup>N</sup> <sup>1</sup> 1 1 4

1 N 1 1 4

<sup>1</sup> 1 1 <sup>N</sup> <sup>4</sup>

In (16), the shape functions correspond to the case of one-dimensional element with three nodes (See Fig. 5a). For two-dimensional case four node elements are used, whose standard

> *,*

vector notation is adjusted, the equations (1) and (2) can be reduced to a matrix system of elementary type of the form (18) which corresponds to the algebraic discretization in e domain of an element, with *k* elementary stiffness matrix, *m* the unknowns and *f* the independent term. Joining the result of (18) for total elements in , a general matrix system is obtained defined as (19) where *K* is the global matrix of stiffness, *M* is the vector of

3

h 2

1

2

3

4

u h hn u h a(u ,W ) b(W ,p ) (g ,W ) (14)

*)* which allow interpolation within

(16)

(17)

*)* by the *Jacobian* of the *transformation* and

<sup>n</sup> n n1 <sup>p</sup> b(u ,q ) t c(p ,q ) b(u ,q ) t(g ,q ) (15)

*,*

*discretization* which together are described by equations (14) and (15):

h h

h h

W

shape functions are shown in matrix form (17) and in Figure 5b.

If the integration space *(x,y)* is changed to *(*

unknowns and *F* is the global vector entries.

where the *N* vector contains *shape functions standard (*

Solution of these equations for both *u* and *p* are implemented using a routine programmed in Fortran and a desktop PC with 2.4 GHz AMD processor and 1.0 GB of RAM

### **4.3. Computer simulation**

Conditions for the simulation were based on experiments by Frijns, 2000 and Wu et al., 1999. Simulations were performed for 1D and 2D looking tissue response to the application of: (a) compression, (b) tensile and (c) oscillating or cyclical loads. The simulation time is 45 seconds of load in all cases. The loads applied in 1D and 2D simulations were performed with the parameters shown in Table 1. The tests considered the AC as a continuous and homogeneous material with the chondrocytes being part of the continuum.


**Table 1.** Parameters for the homogeneous material characteristics (Wu et al., 1999).

### *4.3.1. Meshing*

For the 1D simulation a mesh was performed that represents a fragment of 0.3 mm AC. From this mesh 641 nodes and 320 elements of 3 nodes were obtained. In the 2D simulation a mesh was performed that represents a fragment of articular cartilage 0.09 mm2 (0.3 mm x 0.3 mm). In this case 10,201 nodes were obtained, equivalent to 10,000 elements of 4 nodes.

### *4.3.2. Boundary conditions*

Simulation was performed so as not to allow displacement at the bottom. Load was applied on the upper edge, allowing the fluid outlet only at the bottom of the tissue fragment, as shown in Figure 6a. For 2D, a condition of tissue confinement was simulated, as shown the Figure 6b, so as to present only flow at the bottom. The burden was placed at the top and lateral and bottom movement were restricted, similar conditions to those reported in several experimental studies, including that of Ateshian et al, 1997; Frijns, 2000 and Wu et al, 1999.

Mechanical Behavior of Articular Cartilage 209

**Figure 7.** Response of tissue to compression forces. a. Displacement of solid matrix. b. Changes in the

Figure 8 summarizes the results obtained for tensile strength. Figure 8a shows the positive displacement of the solid component in response to tensile load imposed on tissue. The figure shows how the displacements are positive in presence of tension maintained over time. From an initial value of 0 mmN-1s-1 at t = 0.67s, displacement increases up to a

**Figure 8.** Response of tissue to tensile forces. a. Displacement of solid component tissue upward. b.

Figure 8b shows the increase of fluid pressure generated by the entry thereof into the tissue due to the pressure difference between the internal and external environment (fluid suction). This pressure difference is caused by displacement of the solid component in the positive direction in response to imposed tensile load. It is shown that this progression of the increase of pressure manifests over time, creating a redistribution of the fluid in the

Change in the fluid pressure generated by the redistribution of fluid in the tissue.

fluid pressure in the tissue in presence of the displacement of the same.

maximum of 16x10-4 mmN-1s-1 at a maximum load time t = 45s.

*5.1.2. Tensil loads response* 

tissue.

**Figure 6.** AC scheme for the confinement conditions in (a) 1D and (b) 2D.

### *4.3.3. Loading conditions*

Calculations to simulate the applied load were performed from the data for the AC reported by Wu et al. For the 1D simulations, a load was applied of -2.033 N/m in compression, 2.033 N/m in tension and cyclic loading with frequency equal to 0.1 Hz. In 2D simulations a load was applied on the upper face of cartilage fragment corresponding to the value -4.0397e-6 N/m in compression, 4.0397e-6 N/m in tension. For cyclic loading, it was applied at a frequency equivalent to 0.1 Hz
