1 Introduction

As a fluid-saturated composite tissue, the articular cartilage (AC) maintains two functions. One is to provide a bearing surface, and the other is buffering and distributing the loads between bones with its interstitial fluid. It achieves the functions through its complex structure comprised of a collagen fibril network and proteoglycans (PGs) saturated with an interstitial fluid phase^[1]. The AC is devoid of blood vessels, lymphatics, and nerves. Chondrocytes are the only found cells. Therefore, the AC can be easily damaged and has a limited capacity for intrinsic healing and repair. Surely, the interstitial fluid must play a key role in the injury mechanism.

The mechanical behavior of the AC has been first successfully described by Mow et al.^[2] where the solid matrix and the interstitial fluid were both assumed to be incompressible and nondissipative. The intrinsic biphasic material properties were determined by standard experimental procedures such as uniaxial (confined) compression^[2-4], unconfined compression^[5-6], and indentation^[7-8]. Deformation of cartilage under compression was described reasonably well by various poroelastic models^{[2, 9-11]}, in which the mechanical response is governed primarily by the viscous or rheological stresses generated during the flow of interstitial fluid, and thus intimately related to the changing pore structure of the matrix^[12]. Chondrocytes are the only cell type in the AC and depend on the pore fluid flow to complete the metabolism and mechanotransduction, which involves the multiple biophysical signals^[13], such as solid deformation, osmotic pressure, fluid seepage, shear stress, and streaming potential.

The AC (tissue) often bears the loads coming from the complicated physical environment (walking, sitting, and jumping). It would be in the fluid flow stimulus form acting on the chondrocyte when these external load signals transmitted to the AC microstructure. The cell response would induce the tissue structural remodeling to adapt to the external environment. This is the mechanism of signal transduction and feedback in the AC tissue. However, the scientific quantitative pathway of this mechanism is not very clear, especially the multi-scale transduction mechanism of the fluid stimulating signals is uncertain. These signals include fluid pressure, shear stress, and streaming potential. It is difficult to use experimental studies to quantify these signals in the AC, especially in the local injured or damaged part of the AC. Therefore, in order to approach and solve these scientific problems, the clear distributions of fluid pressure and velocity in the normal or laboratorial testing conditions should be addressed first by using theoretical modeling. However, in most analytical poroelastic models, the relaxation behavior^{[7-8, 14-15]}, friction performance^[16-18], and contact stress^[18-20] of the AC were focused and analyzed, but the completely clear analytical solutions for the interstitial fluid in the AC under the compression experiment remain unsolved. Therefore, the interstitial fluid flowing behavior especially its participating role in the injury mechanism in the AC should be noticed and assessed. In this paper, the pore fluid pressure and velocity distributions in the AC under the experimental conditions are obtained. Moreover, the factors of loading and material parameters are quantified.

2 Mathematically modeling an AC 2.1 Governing equations

According to the theory of poroelasticity, the constitutive equations for the AC can be presented as^[21]

(1)

where σ_ij is the total stress tensor for the entire medium, ε_kl stands for the strain tensor for the solid matrix, and p is the fluid's pore pressure. M_ijkl is the drained elastic stiffness tensor^[22], and α_ij is the Biot coefficient.

The variation in the fluid content can be expressed as

(2)

where ξ is the variation of fluid content per unit reference volume. M is the Biot modulus.

The momentum conservation equation in the medium is

(3)

where u_i^s is the solid displacement component. The derivative of a quantity with respect to the time t is denoted by superposed dot ^., and the second derivative is ^... ρ is the mixture density and has the relation ρ=ϕρ_f+(1-ϕ)ρ_s, in which ϕ is the porosity, and ρ_f and ρ_s are the densities of the fluid phase and the solid skeleton, respectively.

The mass conservation equation for the fluid is

(4)

Darcy's law for the fluid seepage is

(5)

where q_i is the fluid velocity component. k_ij is the permeability tensor (diagonal tensor). μ is the pore fluid dynamic viscosity. p_{, i} represents the partial derivative of pressure with respect to the coordinate.

2.2 Model description and simplification

As shown in Fig. 1, in a cylindrical coordinate system (r, θ, and z), the AC is modeled as a disk. H and R are the thickness and radius, respectively. In this work, the AC is treated as an isotropic material. Thus, M_ijkl is defined by three independent parameters: M₁₁, M₁₂, and G; α_ij and k_ij are both assumed to be diagonal tensors defined by one independent parameter representing by α and k, respectively. In this case, the constitutive equations (1) and (2) can be presented by the following two equations:

(6)

(7)

In the isotropic case, α and M can be written as^[23]

(8)

(9)

where the notations K_s and K_f are the bulk moduli for the solid and fluid phases, respectively. K stands for the drained bulk modulus. It has the relation of K=E/(3(1-2ν)).

For (8), M₁₁, M₁₂, and G can be expressed by the elastic modulus E and Poisson's ratio ν,

(10)

Substituting (5) into (4) and neglecting high order terms lead to the continuity equation for the fluid seepage^[22],

(11)

where ∇² represents the Laplace operator.

In this proposed model, the applied load is a time-harmonic longitudinal strain in the form of ε_zz=ε_z0e^iωt, where ε_z0 is the strain amplitude, and ω is the angular frequency linked to the frequency f by the relation of ω=2πf.

2.3 Boundary cases

In this paper, two boundary cases are proposed. One is unconfined, and the other is confined, as shown in Fig. 2. As shown in Fig. 2(a), the sample is under the conventional unconfined compression test, with a general time-dependent axial load (load control) or an apparent axial strain (stroke control) applied through the perfectly rigid, frictionless, and impermeable end caps. The lateral surface is stress and pressure free for the fluid permeating. As shown in Fig. 2(b), the sample is confined by rigid, frictionless, and impermeable bottom plate and lateral ring while being loaded on top using a porous loading plate. In this case, the lateral displacement is confined and the fluid can flow only through the top porous indenter.

For Case Ⅰ, at r=0,

(12)

(13)

at r=R,

(14)

and at z=0, H,

(15)

For Case Ⅱ, at z=0,

(16)

at z=H,

(17)

and at r=R,

(18)

where σ^Ⅱ(t) is the average reaction stress acting on the top of the material.

3 Problem solving 3.1 Case Ⅰ

For Case Ⅰ, this is an axisymmetric problem. Thus, u_θ=0, and τ_θr=τ_θz=0. Both the radial displacement u_r and the pressure p depend only on the radius r and time t (which leads to γ_zr=0). Moreover, u_r and p are assumed to have the solution form of u_r=u_r0(r)e^iωt and p=p₀(r)e^iωt. Hence, the strain-displacement relation can be written as

(19)

(20)

Only low cyclic loading frequencies are considered (a few Hz). Hence, inertia terms can be neglected, and (3) can be written as

(21)

Substituting (6) into (21) leads to

(22)

Substituting (19) and (20) into (22) yields

(23)

Substituting (19) and (20) into (7) yields

(24)

Substituting (24) into (11) leads to

(25)

Consider u_r=u_r0(r)e^iωt and p=p₀(r)e^iωt. Then, (23) can be expressed as

(26)

where c is an integral constant. Thus, the following differential equation of p₀^Ⅰ(r) is derived from (26) by considering (25):

(27)

The general solution to the above differential equation is

(28)

where I_n and K_n are the first kind and the second kind modified Bessel functions of order n, respectively. The constant C is given by

(29)

Then, the constants A and B can be obtained by taking into account the above boundary conditions (12) and (14),

(30)

(31)

Therefore, the solution to p^Ⅰ(r, t) can be written as

(32)

Next, the last constant c will be determined.

Substituting (28) into (26), u_r0(r) can be expressed as

(33)

The general solution to (25) is

(34)

where s is another integral constant.

According to the boundary equation (13). u_r0|_r=0=0 and (34), s can be obtained to be zero,

(35)

According to (19), (20), and (35), the following two equations are obtained from (34):

(36)

(37)

According to the boundary equation (13), we have

(38)

According to the first equation of (6), (38) can be written as

(39)

Substituting (36) and (37) into (39) leads to (noting r=R)

(40)

Thus, the constant c can be obtained by solving (40) and (30),

(41)

Moreover, the velocity q is computed by modified (4),

(42)

which develops into the following expression:

(43)

The average reaction stress σ^Ⅰ(t) is equal to the integral of the normal stress σ_zz over the sample area,

(44)

According to the third equation of (6), we can get the average reaction stress σ (t),

(45)

3.2 Case Ⅱ

Case Ⅱ is also an axisymmetric problem. Thus, u_θ=0, and τ_θr=τ_θz=0. Both the axial displacement u_z and the pressure p depend only on z and t. Moreover, u_z and p are assumed to have the solution form of u_z=u_z0(z)e^iωt and p=p₀(z)e^iωt. According to u_r=0, the strain-displacement relation can be written as

(46)

(47)

(48)

Thus, the constitutive equations (6) and (7) can be presented by

(49)

(50)

Substituting (50) into (11) leads to

(51)

Consider ε_z=ε_z0e^iωt and p=p₀(z)e^iωt. The above equation (51) becomes

(52)

The general solution to the above differential equation is

(53)

where the constant C' is given by

(54)

By assuming σ_zz=σ₀(z)e^iωt, the constants A' and B' can be obtained by taking into account the above boundary conditions (16) and (17) ((49) is also needed),

(55)

Finally, the pressure solution for Case Ⅱ can be obtained,

(56)

The velocity q_z^Ⅱ is

(57)

The average reaction stress σ^Ⅱ(t) is

(58)

4 Numerical parameters

As shown in Table 1, the isotropic poroelastic and geometrical constants for the AC experimental model are grouped. The results are presented in terms of the amplitudes of fluid pressure and velocity, which have been carried out as a function of (ⅰ) the strain amplitude, (ⅱ) frequency, (ⅲ) strain rate amplitude, (ⅳ) radii/height, and (ⅴ) permeability.

5 Results 5.1 Strain amplitude

Relations between the pore pressure amplitude (p^*=| p(r, t)|), the velocity amplitude (q^*=|q(r, t)|) and the amplitude of strain are plotted in Fig. 3 at a constant loading frequency of 1 Hz. As shown in Fig. 3, both the pressure and velocity amplitudes are proportional to the loading amplitude, as has been displayed in (41), (43), (56), and (57).

5.2 Frequency

Figure 4 is plotted for a fixed longitudinal strain amplitude of ε_z0=0.001. As shown in Fig. 4, regularities of the pressure and frequency amplitudes vs. frequency are not obvious. As shown in Fig. 4(a), the pressure amplitude in Case Ⅰ seems changing little when the frequency goes higher. The pressure amplitude at the AC top (z=2 mm) in Case Ⅱ presents a constant value approximately the smallest value 7.7 MPa (see Fig. 4(c)), while the velocity amplitude is the largest (see Fig. 4(d)) relative to the other places. As shown in Fig. 4(b), the velocity amplitudes in Case Ⅰ decrease gradually as the frequency goes higher.

5.3 Strain rate

The strain rate is defined by the relation of^[5-7]

Thus, (32) and (43) can be rewritten by

(59)

(60)

(61)

(62)

The pressure and frequency amplitudes are proportional to the amplitude of strain rate as shown in (59)-(62). As shown in Fig. 5, both the pressure and velocity amplitudes in Cases Ⅰ and Ⅱ decrease with the frequency increasing from 1 Hz to 21 Hz. It seems that both the pore pressure and fluid velocity amplitudes in the AC depend more on the strain than on the loading frequency.

5.4 Pressure and velocity distributions along the coordinate

Figure 6 shows the variation of pore pressure and fluid velocity amplitudes along the radial ((a) and (b)) and height coordinate at ε_z0=0.01 and f=1 Hz, respectively. According to the imposed boundary conditions, the maximum value of pore pressure amplitude in Case Ⅰ presents (r=0 mm) at the center, and the pore pressure amplitude is zero at the edge of the AC (r=5 mm). On the contrary, the fluid velocity amplitude gets the maximum value at the edge (see Fig. 6(b)). The pressure amplitude in Case Ⅱ (see Fig. 6(c)) mainly presents a decrease along the height of specimen except a turning point at z=1.75 mm. As shown in Fig. 6(d), the fluid velocity amplitude presents a linear relationship with the coordinate.

5.5 Permeability

Finally, the role of permeability is examined in Fig. 7. It is plotted with the strain loading rate fixed at z₀=0.01 and f=1 Hz. In the paper, the permeability coefficient of the AC is on the order of 10^-18 m^2[24]. In our study, we choose 10^-18 m² as a case to obtain the above results (Subsections 5.1-5.5). Besides, this parameter is set between 10^-17 m² and 10^-19 m² to evaluate its influence. As shown in Figs. 7(a) and 7(c), the pressure amplitude in both cases decreases as the permeability increases to 10^-17 m², while the velocity amplitude presents an opposite variation (see Figs. 7(b) and 7(d)) as the permeability increases.

6 Discussion

The injury mechanism research of AC often focuses on the solid friction or mechanical strength, but rarely on the fluid permeate behavior. As a preliminary study, in this paper, the fluid behaviors especially the distributions of pressure and velocity in the AC are obtained under a laboratorial test. The interstitial fluid flowing behavior associated with permeate capacity would be changed when the AC has been injured. Moreover, quantifying the fluid pressure and velocity is difficult in the experiment, and thus in this work it might be done. Some other parameters, such as Poisson's ratio and permeability, can also be obtained.

As shown in Subsection 5.1, a larger loading amplitude can lead to a larger pressure and velocity amplitude in the AC. Especially, in this study, their relationships are linear. This result shows that the external load can make the interstitial fluid with a high pressure to gather and transfuse in the AC tissue. It has been shown that the high loads transmitted across joints during weight bearing are predominantly supported via the hydrostatic pressure developed within the interstitial fluid of the AC^[25-26]. High loads might induce a larger fluid velocity that benefits metabolism of cartilage cells. It should be noted that the effects of loading frequency on the fluid pressure and velocity amplitudes are not as straightforward as the loading amplitude. In the unconfined case Ⅰ, the pressure amplitude changes little as the frequency goes higher, while the velocity amplitude decreases gradually. Therefore, the loading intensity and frequency, which one is the most impact factor that affects the AC fluid flowing behavior? In this study, the strain rate is defined and evaluated. Both the fluid pressure and velocity amplitudes decrease as the frequency grows when the strain rate is fixed, which means that the poroelastic behavior of AC depends more on the loading amplitude than on the frequency. To a certain extent, in order to obtain the AC cell responses, the most effective way is to increase the loading strains. In the same tissue strain level, the increase of the frequency might not be effective to stimulate cells responding associated with the fluid shear stress.

According to the imposed boundary conditions, the maximal fluid pressure in the unconfined case Ⅰ test is at the sample center r=0, while the fluid velocity is minimum of zero. The pressure amplitude at the AC top (z=2 mm) in the confined case Ⅱ test presents constant (approximately 7.7 MPa), while the velocity amplitude is the largest. Especially, the velocity presents a linear increase along the AC specimen height. This is the one-dimensional confined compression test. The deformation and fluid flow move in the direction of the applied load. The top rigid porous-permeable indenter permits the interstitial fluid to flow out, as it exudes from the tissue. It should be noted that the minimum pressure is nearly at 87.5{%} of the specimen height. This would be a pressure release area and would be beneficial to the natural articular. This buffer zone allows the fluid pressure to release and also makes a protection for the whole AC structure from a fracture under an external load. In this study, the boundary conditions are set for the experiment, almost not the same with the physiological conditions. The unconfined case Ⅰ might simulate the condition of AC with a edge fracture. In Case Ⅱ associated with physiological conditions, the permeable boundary condition at top (or bottom) can be seen as the fluid permeating from the subchondral bone to the growth edition. Thus, the interface could likely be the ``porous indenter''. In the natural AC with the complex micro porous structure, the fluid flowing behavior is tanglesome but orderly from the high pressure region to the low one, along the load induced pressure gradient directions.

The important material parameters of the AC, i.e., the solid phase and the fluid phase, have been shown to be effectively incompressible^[27]. Therefore, Poisson's ratio is particularly set to be 0.499. For the cortical bone, the solid matrix (14 GPa) has a bulk modulus 6 times that of water (2.3 GPa), which leads to a compressibility coefficient B

(63)

of 0.40^[24]. In the case of soft tissues, the bulk moduli of the soft tissue matrix and the pore water are almost the same and the compressibility B is about 1.0. In our study, it is about 0.9. The other important material parameter would be the permeability. Its order is about 10^-18 m². In this paper, this parameter is fixed at 2.74×10^-18 m². This parameter would not be a constant with the AC's nonlinear material properties. In an experiment study^[28], the normal permeability of the AC superficial, middle, and deep layers were tested with the value of 1.03×10^-18 m², 1.10×10^-18 m², and 0.86×10^-18 m², respectively. The permeability plays an important role in the fluid flowing behavior. As shown in Figs. 6(b) and 6(d), both the fluid velocity amplitudes in two cases increase as the permeability value increases, while the pressure amplitude decreases as the permeability increases approach 10^-17 m². In most studies, this parameter is set on the order of 10^-18 m². The frictional property of AC is also influenced by the permeability. The reduced permeability results in decreasing friction coefficient in the AC^[28]. Though the permeability is difficult to quantify, several estimations of the permeability have been proposed, ranging from 10^-19 m² to 10^-17 m². For the cortical bone associated with the lacuno-canalicular system (LCS), this parameter lays between 10^-25 m² and 10^-20 m^2[29-31].

7 Conclusions

This analytical poroelastic solution links the external loads to the interstitial fluid pressure and velocity, which may be a stimulus to the AC injury mechanism involved fluid flow. However, simplifications and limitations in this ideal AC model are obvious. First, the geometrical disc shape and the boundary conditions are just for the experimental tests, and it is far away from the natural AC. Second, the nonlinear permeability and elastic properties are not considered in this model. However, some basic material parameters can be obtained by combining this theoretical solution with the experiment. A more detailed model taking these disadvantages into account may contribute to a better understanding of mechanotransduction associated with the fluid stimulate, as well as the mechanisms in AC injury or repairing processes. Though there are limitations, some main conclusions may be drawn below.

(ⅰ) Both the pore pressure and fluid velocity amplitudes are proportional to the strain loading amplitude.

(ⅱ) In a physiological loading state, it seems that both the amplitudes of pore fluid pressure and velocity in the AC depend more on the loading amplitude than on the frequency. Thus, in order to obtain the considerable fluid stimulus for the AC cell responses, the most effective way is to increase the loading amplitude rather than the frequency.

(ⅲ) Both the pressure and fluid velocity are strongly affected by the permeability variations.