I I Argatov
email: argatov@home.ru

G S Mishuris

Axisymmetric contact problem for a biphasic cartilage layer with allowance for tangential displacements on the contact surface

Keywords: Contact problem, cartilage layer, tangential displacements, analytical solution

A unilateral axisymmetric contact problem for a biphasic cartilage layer indented by a rigid punch is considered. The refined linearized kinetic relationship which takes into account the both radial and tangential displacements of the boundary points of the biphasis cartilage layer is imposed. The obtained analytical solution is valid over long time periods and can be used for increasing loading conditions. It can be used as it is and as a benchmark for verification of FEM model accuracy.

Introduction

Biomechanical contact problems involving transmission of forces across biological joints are of considerable practical interest. Many solutions to the problem of contact interaction of articular cartilage surfaces in joints are available. [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF] obtained an asymptotic solution for the contact problem of two identical biphasic cartilage layers attached to two rigid impermeable spherical bones of equal radii. [START_REF] Wu | Modeling axi-symmetrical joint contact with biphasic cartilage layers-An asymptotic solution[END_REF] extended this solution to a more general model by combining the assumption of the kinetic relationship from classical contact mechanics [START_REF] Johnson | Contact Mechanics[END_REF] with the joint contact model for the contact of two biphasic cartilage [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. An improved solution for the contact of two biphasic cartilage layers which can be used for dynamic loading was obtained by Wu et al. (1997a).

It is known [START_REF] Wu | Joint contact mechanics in the early stages of osteoarthitis[END_REF] that in real joints, the structural and biomechanical properties of the cartilage may change dramatically in the pathologic stage. Moreover, the geometrical complexity of the problem cannot give a chance to solve it analytically in the general setting and the numerical approach (mostly with FEM) moves to the frontier of the research and becomes recently the main tool for dealing with wider spectrum of medical application [START_REF] Wilson | The role of computational models in the search for the mechanical behaviour and damage mechanisms of articular cartilage[END_REF]. It is a common knowledge, however, that the numerical simulations should be carefully verified. As in any practical problem, experimental results could be the most efficient verification. Unfortunately, in most types of biomechanical contact problems, where experiments cannot be easily performed in vivo, other methods to verify FEM models become increasingly important. And if comparison of the numerical results obtained with different FEM packages seems to be the way to do so, the most efficient way still remains a comparison with the exact analytical results. In their article, Wu et al. (1997b) compared the solutions obtained using ABAQUS with other finite element models and with those obtained using analytical solutions.

Moreover, in the case of a thin contrast interface (with properties essentially different from those in the adherents) the FEM code often fails, that is its reliable application requires so a fine mesh inside the interface and in its neighborhood (see, for example, [START_REF] Mishuris | 2D modelling of a thin elasto-plastic interphase between two different materials: Plane strain case[END_REF]) that it turns out to be unrealistic to obtain any credible results especially if time dependent problems are considered. For this reason, asymptotic methods recently come back to the game, and many researchers considered various revised problems which previously have been tackled in the most simple formulations [START_REF] Benveniste | An interface model of a graded three-dimensional anisotropic curved interphase[END_REF][START_REF] Benveniste | An interface model for a three-dimensional curved thin piezoelectric interphase between two piezoelectric media[END_REF].

Relatively recently, the tangential displacements effect was found to be quite important in the nanoindentation technique for determining the elastic properties. Observing that the so-called Bulychev -Alekhin -Shorshorov relation [START_REF] Bulychev | Mechanical properties of materials studied from kineticdiagrams of load versus depth of impression during microimpression[END_REF] for analyzing nanoindentation load-displacement data does not take into account the radial displacements of material under the indenter, [START_REF] Hay | Critical examination of the fundamental relations used in the analysis of nano-indentation data[END_REF] introduced the corresponding correction factor, which turned out to be Poisson's ratio dependent. For spherical indenters, [START_REF] Hay | Small correction required when applying the Hertzian contact model to instrumented indentation data[END_REF] have shown that the additive correction is proportional to the ratio a/R, where a is the contact radius and R the indenter radius. The refined nonlinear boundary conditions for the actual contact were discussed by [START_REF] Fu | On Sneddon's boundary conditions used in the analysis of nanoindentation data[END_REF]. Therefore, the refined kinetic relationship which takes into account the radial tangential displacements of the boundary points of a biphasic cartilage layer should be imposed instead of the classical kinetic relationship for correct predictions of biomechanical parameters in the contact region of the articular cartilage.

The purpose of this study is to revisit the existing analytical solutions for the axisymmetric analytical contact problem of a biphasic cartilage layer with a rigid punch and to construct a new more accurate solution of the model contact problem which is also capable to represent the limit situation when the biomechanical properties of the contacting cartilage layers differ drastically. The obtained solution is valid over long time periods and under increasing loading conditions. It can be used as it is and as a benchmark for verification of FEM model accuracy.

The analytical approach proposed here is directly based on the methods previously developed by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF], Wu et al. (1997a), andArgatov (2004). Note that analytical approaches to various contact problems related to biomechanics were developed previously by [START_REF] Eberhardt | Normal contact of elastic spheres with two elastic layers as a model of joint articulation[END_REF] and [START_REF] Haider | Indentation of a thin compressible elastic layer: Approximate analytic and numerical solutions for rigid flat indenters[END_REF]. Recent contributions relevant to the problem under consideration approach are those of [START_REF] Pavliotis | A perturbation-based numerical method for solving a three-dimensional axisymmetric indentation problem[END_REF] and [START_REF] Batra | Analytical solution of the contact problem of a rigid indenter and an anisotropic linear elastic layer[END_REF]. However, to the best of the authors' knowledge, this is the first time that the biomechanical contact problem has been solved in the refined formulation with the tangential displacements taken into account.

In solving contact problems for a punch indenting an elastic body which are formulated in the frame of the linear theory of elasticity, it is a standard assumption that the tangential displacements of the boundary points of the elastic body that come in contact with the punch are negligible (see, e. g., [START_REF] Gladwell | Contact Problems in the Classical Theory of Elasticity[END_REF]; [START_REF] Johnson | Contact Mechanics[END_REF]). This simplified assumption is employed to uncouple the normal and tangential displacements along the contact interface, thus considerably simplifying the resulting system of integral equations which has to be solved. The contact problems with allowance for the tangential displacements at the contact zone were considered by [START_REF] Brock | Frictionless indentation by a rigid wedge: The effect of tangential displacements in the contact zone[END_REF], [START_REF] Georgiadis | Tangential-displacement effects in the wedge indentation of an elastic half-space -an integral-equation approach[END_REF], Argatov (2004). It was shown that the refined contact problem formulation implies new qualitative features in singular contact problems with geometrical discontinuities of the indenter [START_REF] Brock | Frictionless indentation by a rigid wedge: The effect of tangential displacements in the contact zone[END_REF][START_REF] Georgiadis | Tangential-displacement effects in the wedge indentation of an elastic half-space -an integral-equation approach[END_REF], whereas the tangential-displacement effects in the contact problem for two elastic bodies are proportional to the relative difference between the elastic moduli of the contacting bodies (Argatov, 2004). To the best of the authors' knowledge, this paper presents the first investigation of the tangential-displacement effects in the contact problem for a biphasic cartilage layer. The obtained results were partially presented at the 1st International Conference on Mathematical and Computational Biomedical Engineering1 .

Note that, as a rule, simplifying hypotheses should allow to obtain results at least qualitatively comparable with experimental results. The simplified shape of the punch and the mechanical properties of the cartilage are justified by the aim of the study. In fact, in contrast to the simplifying geometry we introduce the refined contact condition. The numerous simplifying hypotheses which appear throughout the derivation of the approximate analytical solution are verified by the numerical calculations on the basis of the exact solution.

Refined formulation of the contact problem

We consider a thin linear biphasic cartilage layer indented by a punch shaped like a body of revolution. Introducing the cylindrical coordinate system (r, ϕ, z), we write the equation of the punch surface (before loading) in the form z = h + Φ(r), where h is the thickness of cartilage layer. For simplicity, we assume that the punch occupies a convex domain z ≥ h + Φ(r) whereas it is in contact with the plane z = h at a single point chosen as the coordinate origin. In the particular case of a punch shaped like a paraboloid of revolution, we have

Φ(r) = r 2 2R 0 , (1)
where R 0 is the curvature radius of the punch surface at its apex.

We denote the vertical displacement of the punch by -δ 0 (t). The unilateral contact condition that the boundary points of the cartilage layer do not penetrate into the punch can be written as follows [START_REF] Kravchuk | On the Hertz problem for linearly and nonlinearly elastic bodies of finite dimensions[END_REF]:

w(r, t) + δ 0 (t) -Φ r + u(r, t) ≤ 0.
(2)

According to [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF], the horizontal and vertical displacements of the boundary points of the tissue, u(r, t) and w(r, t), are expressed in the following dimensionless form:

u(r ′ , t ′ ) = -h ε 2µ s ∂p(r ′ , t ′ ) ∂r ′ , (3)
w(r ′ , t ′ ) = h ε 2 µ s 1 3 1 r ′ ∂ ∂r ′ r ′ ∂p(r ′ , t ′ ) ∂r ′ + 1 H t ′ 0 1 r ′ ∂ ∂r ′ r ′ ∂p(r ′ , τ ) ∂r ′ dτ . (4)
Here, ε = h/ā 0 is a dimensionless small parameter with ā0 denoting a characteristic measure of the contact radius, H = (λ s + 2µ s )/µ s is a material parameter of the cartilage, where λ s and µ s represent the first Lamé coefficient and the shear modulus of the solid phase of the cartilage tissue, p(r ′ , t ′ ) is the contact pressure. The dimensionless radial coordinate and time, r ′ and t ′ , are defined by the formulas

r ′ = r ā0 , t ′ = t τ ′ , τ ′ = h 2 (λ s + 2µ s )k , (5)
where k denotes the cartilage permeability.

The equality in relation (2) determines the contact radius a(t). In other words, the following equation holds within the contact area: w(r, t) + δ 0 (t) -Φ r + u(r, t) = 0 for r ≤ a(t). We assume that the displacement u(r, t) is small compared to the contact radius a(t). In this case, the written-out nonlinear contact condition can be replaced by the following linearized kinetic relationship:

w(r, t) + δ 0 (t) -Φ(r) -Φ ′ (r)u(r, t) = 0, r ≤ a(t). (6)
Here the prime denotes differentiation.

In view of formula (1), Eq. ( 6) takes the form

-w(r, t) + r R 0 u(r, t) = δ 0 (t) - r 2 2R 0 , r ≤ a(t). (7)
Expressing the displacements w(r, t) and u(r, t) in the dimensional form by means of formulas (3) -( 5), we write the displacement compatibility condition (7) in the following form (we assume that r ≤ a(t)):

1 r ∂ ∂r r ∂p(r, t) ∂r + χ t 0 1 r ∂ ∂r r ∂p(r, τ ) ∂r dτ + κr ∂p(r, t) ∂r = m Cr 2 -δ 0 (t) . ( 8
)
Here we introduced the notation

χ = 3µ s k h 2 , κ = 3 2hR 0 , m = 3µ s h 3 , C = 1 2R 0 . (9)
Eq. ( 8) is used to find the contact pressure density p(r, t). The radius of the contact area a(t) is determined from the condition that the contact pressure is positive and vanishes at the contour of the contact area:

p(r, t) ≥ 0, r ≤ a(t); p(a(t), t) = 0. (10)
Moreover, the contact pressure should satisfy the boundary conditions ∂p(r, t) ∂r r=0

= 0, ∂p(r, t) ∂r r=a(t) = 0,

where the first one comes from the standard idea for an axisymmetric blunt punch, whereas the second one appears in the case of contact problem for a biphasis cartilage layer, where it is additionally assumed a smooth transition of the surface normal stresses from the contact region r < a(t) to the outside region r > a(t) [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. Finally, the equilibrium equation for the whole system is 2π

a(t) 0 p(ρ, t)ρ dρ = F (t), (12)
where F (t) denotes the external load.

For non-decreasing loads when dF (t)/dt ≥ 0, the contact radius increases monotonously, i. e., the following condition holds (Wu et al., 1997a): da(t)/dt > 0.

It is readily seen from formulas (3) and (4) that the tangential displacement w(r ′ , t ′ ) is an order of magnitude greater than the normal displacement u(r ′ , t ′ ). Thus, the formulation of the contact problem excluding the tangential displacement comes in contradiction with neglecting terms that could contribute significantly in certain configurations. This consideration implies the question whether and how the tangential displacements can influence the parameters of the contact pressure distribution. It is clear that the inclusion of the tangential displacements makes the formulation of contact problem more precise though the obtained refined problem becomes much more difficult, and requires a much more elaborate approach for analytical solution.

The aim of this study is to derive a solution for the refined contact problem for a biphasic cartilage layer in the case of increasing loading condition handling accurately both radial and tangential displacements of the cartilage boundary. The solution proposed here is directly based on the approach developed by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF], Wu et al. (1997a), andArgatov (2004).

In the next section we formulate the most important obtained results concerning asymptotic relationships between the main biomechanical parameters. In Section 3, we derive and solve equations for the punch displacement and the radius of the contact area. In Section 4, we obtain analytical representations for the contact pressure density. In Section 5, we show that the contact problem under consideration coincides in the case κ = 0 with that studied in detail by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF], Wu et al. (1997a). Overall, this paper addresses the question when accounting of the additional term in Eq. ( 8) is important for applications.

Some of the technical derivations useful for understanding and further applications are included in Appendices.

Outline of main results

Based on the governing integral equation ( 8), we have established two a priori relationships between the three main parameters of contact. Namely, the following equations connect the unknown radius of the contact area a(t), the unknown punch displacement δ 0 (t) and the contact load F (t):

δ 0 (t) = C a 2 (t) 2 + 2κ m F (t) πa 2 (t) , ( 13
)
πm 48 Ca 6 (t) = KF (t) + κ      a 2 (t) 4 F (t) + π 2 a(t) 0 ρ 4 ∂p(ρ, t) ∂ρ dρ      . ( 14
)
Here we have introduced the following short hand notation for the Volterra integral operator K on the right-hand side of Eq. ( 14) and its inverse K -1 :

KF (t) = F (t) + χ t 0 F (τ ) dτ, K -1 f (t) = f (t) -χ t 0 f (τ )e -χ(t-τ ) dτ. ( 15
)
In a specific case κ = 0 Eqs. ( 13) and ( 14) coincide with the results obtained by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]; Wu et al. (1997a). We underline that these equations are exact consequences of the governing integral equation ( 8).

In order to solve Eq. ( 14) and derive an approximation for the radius of the contact area, we have invented a new asymptotic method based on the recursive relation for the integral characteristics of the contact pressure density. The following equation takes into account the first order correction to the effect of tangentional displacements on the contact radius:

πm 48 C a 6 (t) + κ 4 K -1 a 8 (t) = KF (t) + κa 2 (t) 4 F (t). ( 16
)
In the domain 0 < r < a 0 , the following approximate formula for the derivative of the contact pressure distribution takes place:

∂p - ∂r (r, t) = mr 8 2Cr 2 e -χt -4K -1 δ 0 (t) (17) + κr 2 C 3 r 2 e -χt (χt -1) + K -2 δ 0 (t) + O(κ 2 ), κ → 0.
In the interval a 0 < r < a(t), we have the following asymptotic representation of the derivative of the contact pressure distribution:

2r m ∂p + ∂r (r, t) = a(t) -r C 2 a 2 0 (2r 2 -a 2 0 )e -χt -r 2 K -1 δ 0 (t) (18) + κ 4 r 4 K -2 δ 0 (t) + 1 3 Ca 2 0 (χt -1)(3r 4 -3r 2 a 2 0 + a 4 0 )e -χt (19) + R(r, t, κ) + O(κ 2 ), κ → 0.
Here < • > is the Heaviside step function defined in the standard way: x = 1 if x ≥ 0 and x = 0 if x < 0 and the auxiliary integral R(r, t, κ) is determined by formulas ( 61) and ( 62).

The main difficulty in calculating the contact pressure is the variability of the contact area that complicates the application of the Laplace transform. A sophisticated numerical algorithm have been elaborated for dealing with integral equations with variable domain.

3 Equations for the punch displacement and the radius of the contact area

Let us integrate Eq. ( 8) with respect to r and transform the obtained integral using integration by parts. Since the constant of integration D 1 (t) vanishes in accordance with the regularization condition ( 11) 1 , we get

K ∂p(r, t) ∂r + κ rp(r, t) - 2 r r 0 p(ρ, t)ρ dρ = m C r 3 4 -δ 0 (t) r 2 , (20)
which, after substituting the value r = a(t) and taking into account the boundary conditions ( 10), ( 11) as well as the monotonicity condition da(t)/dt > 0, transforms into

- 2κ a(t) a(t) 0 p(ρ, t)ρ dρ = m C a 3 (t) 4 -δ 0 (t) a(t) 2 . (21)
Finally, in view of formula (12), Eq. ( 21) takes the form of Eq. ( 13).

We now return to Eq. ( 20) which after integration with respect to variable r take the forms:

Kp(r, t) + κ    r 0 p(ρ, t)ρ dρ -2 r 0 p(ρ, t)ρ ln r ρ dρ    = m C r 4 16 -δ 0 (t) r 2 4 + D 2 (t). ( 22
)
By using the boundary condition (10), the constant of integration can be obtained as a function of the contact radius a(t):

D 2 (t) = m 4 δ 0 (t)a 2 (t) - m 16 Ca 4 (t) + κ 2π F (t) -2κ a(t) 0 p(ρ, t)ρ ln a(t) ρ dρ. ( 23
)
Multiplying both sides of Eq. ( 22) by r and integrating over the contact radius, we obtain

1 2π KF (t) + κ      a(t) 0 ρ ρ 0 p(ξ, t)ξ dξdρ -2 a(t) 0 ρ ρ 0 p(ξ, t)ξ ln ρ ξ dξdρ      = m C a 6 (t) 96 -δ 0 (t) a 4 (t) 16 + a 2 (t) 2 D 2 (t). ( 24
)
We change the order of integration in the second and third integrals in (24). After that, taking into account formulas ( 12), ( 23), and (13), we transform Eq. ( 24) into the following equation:

πm 48 Ca 6 (t) = KF (t) + κ      a 2 (t) 4 F (t) -2π a(t) 0 p(ρ, t)ρ 3 dρ      . ( 25
)
For further convenience we rewrite condition (25) in the equivalent form ( 14). Notice that Eq. ( 14) was obtained by integrating by parts in the integral with taking into account the condition (10) 2 .

Contact pressure

Returning to Eq. ( 8), we introduce a new dependent variable y(r, t) such that y(r, t) = r ∂p(r, t) ∂r .

(26)

Using the substitution (26), we represent Eq. ( 8) in the form 1 r ∂y(r, t) ∂r

+ χ t 0 1 r ∂y(r, τ ) ∂r dτ + κy(r, t) = mf (r, t) . ( 27
)
where

f (r, t) = a(t) -r Cr 2 -δ 0 (t) . ( 28
)
This allows us to rewrite the equation in the form:

1 r 1 + χ s dỹ dr + κ ỹ = m f (r, s), (29)
where the tilde denotes the time Laplace transform and s is the Laplace transform variable. Note that all these functions should be analytic in the half-plane ℜs > 0.

Taking into account the regularity condition (11) 1 , which now reads ỹ(r, s) = O(r 2 ) as r → 0, we obtain the corresponding particular solution of Eq. ( 29) for arbitrary r > 0:

ỹ(r, s) = m κ e -K(s) r 2 2 r 0 e K(s) ξ 2 2 ξK(s) f (ξ, s)dξ. ( 30
)
Here we introduced the notation

K(s) = κs s + χ . ( 31
)
Under the assumption utilized in this paper that the loading is monotonic (as a result, a(t) is an increasing function) one can compute f (r, s) to have:

f (r, s) =      Cr 2 s -δ0 (s), 0 ≤ r ≤ a 0 , Cr 2 e -sb(r) s - ∞ b(r) e -st δ 0 (t)dt, r ≥ a 0 . b(r) =      0, 0 ≤ r ≤ a 0 , a -1 (r), a 0 ≤ r. ( 32
)
where we denote by a -1 (r) the inverse to the function a(t).

Please include Figure 1 somewhere here.

Not that the solution to Eq. ( 27) is different from zero only in the domain Ω -Ω + drawn on Fig. 1. We shall consider this problem in the both domains separately.

4.1 Solution of Eq. ( 27) in domain Ω - Let us denote by y -(r, t) the solution of Eq. ( 27) in the domain Ω -(0 ≤ r ≤ a 0 (0)). After some algebra, we can rewrite Eq. ( 30) in the interval 0 ≤ r ≤ a 0 as follows:

κ m ỹ-(r, s) = Cr 2 2s + 2C sK(s) + δ0 (s) exp -K(s) r 2 2 -1 . (33)
Taking into account the asymptotic estimate ỹ(r, s) = O(r 2 ) as r → 0, we obtain

κ m ∂ p-(r, s) ∂r = Cr s + 2C sK(s) + δ0 (s) exp - κr 2 2 1 + L(s, r) r - 1 r . ( 34
)
Here we have introduced the following notations:

L(s, r) = exp {G(s, r)} -1, G(s, r) = κχr 2 2(s + χ) . ( 35
)
Then the direct inverse Laplace transform formula gives

κ m ∂p -(r, t) ∂r = Cr + exp -κr 2 2 -1 r 2C κ (1 + χt) + δ 0 (t) (36) + 1 r exp - κr 2 2 t 0 L(t -τ, r) 2C κ (1 + χτ ) + δ 0 (τ ) dτ.
The inverse to the function L(s, r) is calculated in the following form:

L(t, r) = 1 2πi i∞+η -i∞+η exp {G(s, r)} -1 e st ds, (37)
where -χ < η ≤ 0 is an arbitrary number. Properties of this function and discussion of the methods for its effective computation can be found in Appendix 1. Using (70), one can check that the condition ( 11) is valid for any solution determined by Eq. ( 34). Now, we can compute pressure distribution at any particular time t = t * employing formula (36). However, for large time this equation is not the best representation from the numerical point of view. In fact, as it was shown in ( 53), the contact area increases with time to be estimated as a(t) = O(t 1/6 ) as t → ∞. Then, from Eq. ( 13) we get a(t) = O(t 1/3 ) as t → ∞ for the fixed total load F (t) = F 0 . Then the second and the third terms in the solution (36) tend separately to infinity as t → ∞, whereas the pressure distribution, p(r, t), has to decrease with time, because the total contact load F 0 is fixed. To confirm this we integrate in the integral in (36) by parts with taking into account identity (40) to obtain

κ m ∂p -(r, t) ∂r = Cr + exp -κr 2 2 -1 r - 1 r M (0, r) exp - κr 2 2 2C κ (1 + χt) + δ 0 (t) + 1 r M (t, r) exp - κr 2 2 2C κ + δ 0 (0) (38) + 1 r exp - κr 2 2 t 0 M (t -τ, r) 2Cχ κ + δ ′ 0 (τ ) dτ.
Here the function M (t, r) is defined similarly to the kernel function L(t, r):

M (t, r) = 1 2πi i∞-η -i∞-η 1 s exp G(s, r) -1 e st ds. (39)
In fact, there is the close relationship between them:

∂ ∂t M (r, t) = L(r, t). ( 40
)
Properties of the function M (t, r) are discussed in Appendix 2. With use of (75), Eq. ( 38) can be simplified as follows:

κ m ∂p -(r, t) ∂r = Cr + M (t, r) exp - κr 2 2 2C κ + δ 0 (0) + 1 r exp - κr 2 2 t 0 M (t -τ, r) 2Cχ κ + δ ′ 0 (τ ) dτ. (41)
One can see that, in contrast with the representation (36), all the terms in Eq. ( 41) are bounded as t → ∞.

Note that Eqs. ( 36) and ( 41) have formally sense not only in the interval 0 ≤ r ≤ a 0 but for any r ∈ (0, a(t)). However, this fact does not mean that they necessarily represent the solution of Eq. ( 27) for any given function δ 0 (t) in the whole domain 0 ≤ r ≤ a(t). For this reason, we have to return back to Eq. ( 27) and consider it in the interval a 0 ≤ r ≤ a(t) independently.

4.2 Solution of Eq. ( 27) in domain Ω + If r ≥ a 0 , Eq. ( 30) can be rewritten in the following form using integrating by parts:

κ m ỹ+ (r, s) = f (r, s) -e -κ r 2 2 e G(s,r) f (0, s) -e -K(s) r 2 2 a 0 0 e K(s) ξ 2 2 ∂ ∂ξ f (ξ, s)dξ -e -K(s) r 2 2 r a 0 e K(s) ξ 2 2 ∂ ∂ξ f (ξ, s)dξ. (42) Note that ∂ ∂ξ f (ξ, s) = 2Cξ s , 0 ≤ ξ ≤ a 0 , (43)
and the first integral in ( 42) can be computed to obtain

κ m ỹ+ (r, s) = f (r, s) -e -κ r 2 2 e G(s,r) f (0, s) - 2C sK(s) e -κ r 2 -a 2 0 2 e G(s, √ r 2 -a 2 0 ) -e -κ r 2 2 e G(s,r) - r a 0 e -κ r 2 -ξ 2 2 e G(s, √ r 2 -ξ 2 ) ∂ ∂ξ f (ξ, s)dξ. ( 44
)
To calculate the remaining integral in ( 44), we need first to find the derivative f ′ ξ (ξ, s) in the interval r > a 0 . From Eq. ( 32) it follows that

∂ ∂ξ f (ξ, s) = e -sb(ξ) 2Cξ s -b ′ (ξ) Cξ 2 -δ 0 (b(ξ)) .
The following identities can be established by direct integration

1 2πi i∞+η -i∞+η e st ∂ ∂ξ f (ξ, s)ds =      2Cξ, t > b(ξ) 0, 0 ≤ t ≤ b(ξ)      -b ′ (ξ) Cξ 2 -δ 0 (b(ξ)) δ(t -b(ξ)) (45) =      2Cξ, ξ < a(t) 0, ξ ≥ a(t)      -b ′ (ξ) Cξ 2 -δ 0 (t) δ(ξ -a(t)). ( 46
)
Here δ(x) is the Dirac delta function. In the domain t > b(ξ), we have

1 2πi i∞+η -i∞+η e st e G s, √ r 2 -ξ 2 -1 ∂ ∂ξ f (ξ, s)ds = 2Cξ e κ r 2 -ξ 2 2 -1 + M t -b(ξ), r 2 -ξ 2 (47) -b ′ (ξ) Cξ 2 -δ 0 (b(ξ)) L t -b(ξ), r 2 -ξ 2 .
Notice that this function is equal to zero when t ≤ b(ξ) or that is equivalent when ξ ≥ a(t).

It remains now to invert Eq. ( 44) with taking into account ( 45) and ( 47). For a 0 < r < a(t), we have

κr m ∂p + ∂r (r, t) = Cr 2 -δ 0 (t) + e -κ r 2 2 2C κ (1 + χt) (48) +δ 0 (t) + t 0 L(t -τ, r) 2C κ (1 + χτ ) + δ 0 (τ ) dτ - 2C κ e -κ r 2 -a 2 0 2 (1 + χt)+ t 0 L(t -τ, r 2 -a 2 0 )(1 + χτ )dτ -2C r a 0 e -κ r 2 -ξ 2 2 ξM (t -b(ξ), r 2 -ξ 2 )dξ + r a 0 e -κ r 2 -ξ 2 2 L(t -b(ξ), r 2 -ξ 2 ) × b ′ (ξ) Cξ 2 -δ 0 (b(ξ)) dξ -C(r 2 -a 2 0 ).
Taking advantage from asymptotic estimate (70), one can check that

y -(a 0 -, t) = y + (a 0 +, t), ∂y - ∂r (a 0 -, t) = ∂y + ∂r (a 0 +, t). ( 49
)
We also can rewrite the solution representation (48) for r > a 0 in the same way as it has been done in (41) for 0 < r < a 0 . We have

κ m ∂p + ∂r = C r a 2 0 + 1 r e -κ r 2 2 M (t, r) 2C κ + δ 0 (0) (50) + 1 r t 0 M (t -τ, r) 2Cχ κ + δ ′ 0 (τ ) dτ - 2C κr e -κ r 2 -a 2 0 2 M (t, r 2 -a 2 0 ) + χ t 0 M (t -τ, r 2 -a 2 0 )dτ - 2C r r a 0 e -κ r 2 -ξ 2 2 ξM (t -b(ξ), r 2 -ξ 2 )dξ + 1 r r a 0 e -κ r 2 -ξ 2 2 L(t -b(ξ), r 2 -ξ 2 )b ′ (ξ) Cξ 2 -δ 0 (b(ξ)) dξ.

Asymptotics of the initial contact problem and some qualitative conclusions

Expanding the expression for the initial contact pressure p(r, 0) ≡ p -(r, 0), where p -(r, 0) is given by formula (87), into the power series of the parameter κ, we obtain

1 m p(r, 0) = δ 0 (0) 4 (a 2 0 -r 2 ) - C 16 (a 4 0 -r 4 ) -κ δ 0 (0) 32 (a 4 0 -r 4 ) - C 144 (a 6 0 -r 6 ) + O(κ 2 ).
In view of formula (90), we can rewrite the obtained asymptotics as follows:

1 m p(r, 0) = C 16 (a 2 0 -r 2 ) 2 1 + κa 2 0 36 1 -4 r 2 a 2 0 + O(κ 2 ). ( 51
)
It is easy to show that under the non-burdensome condition κa 2 0 ≤ 12 the initial contact pressure given by formula ( 51) is positive. Observe that the contact pressure p(r, 0) in the central region tends to increase with the parameter κ, whereas the accounting for the tangential displacements leads to a decrease in the contact pressure in the boundary region of contact.

Further, expanding the right-hand side of Eq. ( 93) into the power series of the parameter κ, we find 2π m

a 0 0 p(ρ, 0)ρ 3 dρ = πCa 8 0 192 1 - 7κa 2 0 20 + κF 0 a 4 0 12m + O(κ 2 ).
Substituting this asymptotic expansion into Eq. ( 92), after some algebra, we obtain

πm 48 Ca 6 0 = F 0 + κ 2 F 0 a 4 0 240 + O(κ 3 ). ( 52
)
We stress that the first-order term is absent in Eq. ( 52). This explains the paradoxical fact that formula (52) predicts a small increase in the initial radius of contact while at the same time the contact pressure at the central region of contact also increases with the increase in the parameter κ.

Since we consider a monotonic loading (the function F (t) does not decrease with time), a(t) will increases up to infinity what proves by contradiction. Assuming reasonable growth of the function F (t) that is F (t) → F ∞ as t → ∞, the following estimate holds true:

a(t) ∼ 48χF ∞ t πmC 1/6 , t → ∞. ( 53
)
As a result, one can write asymptotic expansion for the other sought-for function δ 0 (t):

δ 0 (t) ∼ 48χF ∞ t πmC 1/3 , t → ∞, (54)
Finally, from formula (90) it is easy to see that the accounting for the tangential displacements results in an increase in the punch displacement δ 0 (0).

Approximate solution to the problem for arbitrary time

In the previous three sections we have found basic equations solving the problem under consideration. They respectively are: (13), ( 25), ( 41), (50). Unfortunately, they require substantial computational efforts to reconstruct the exact numerical solution. However, in a particular case when the parameter κ in its dimensionless form is reasonably small, an approximate solution to the problem will provide acceptable answer. Notice that equation ( 13) defining the punch indentation δ 0 (t) does not need any simplification.

Approximate equation for the radius of the contact area

In order to obtain an approximation for the radius of the contact area, we have derived the recursive relation for the integral characteristics of the contact pressure density (see Appendix 3). In particular, for n = 1, Eq. ( 84) takes the form

KF 1 (t) = πmCa 8 (t) 192 + κ - a 4 (t) 24 F 0 (t) + 3 8 F 2 (t) , (55)
where F n (t), n = 0, 1, 2, are integral characteristics of the contact pressure density determined by formula (78).

Neglecting the last term in the right-hand side of Eq. ( 55), we derive

F 1 (t) = πmC 192 K -1 a 8 (t). ( 56
)
Now, substituting the expression (56) into the right-hand side of Eq. ( 79), we arrive at the approximate equation for the radius of the contact area (16).

In principle, the resolving equation for a(t) can be made more accurate by employing the interactive process based upon formula (84).

For example, using formulas ( 55) and ( 84), n = 2, we have

F 1 (t) = πmC 192 K -1 a 8 (t) + κ 24 K -1 -a 4 (t)F (t) + 3πmC 160 K -1 a 10 (t) . ( 57
)
It is readily seen that formula (57) refines formula (56) by means of the second order terms and leads to

πm 48 C a 6 (t) + κ 4 K -1 a 8 (t) + 3κ 2 80 K -2 a 10 (t) (58) = KF (t) + κ 4 a 2 (t)F (t) + κ 2 24 K -1 a 4 (t)F (t) .
Note that for t = 0 the condition ( 16) transforms with taking into account evident identities Kf (0) = f (0) and

K -1 f (0) = f (0) into the equation πmC 48F (0) a 6 0 = 1, (59)
which is consistent with Eq. ( 25).

Approximate formula for the contact pressure

Let us consider the limit value (as κ → 0) of Eq. ( 36) in the domain 0 < r < a 0 . We can show that

κ m ∂p -(r, 0) ∂r = Cr - κr 2 - κ 2 r 3 8 2C κ (1 + χt) + δ 0 (t) + 1 r 1 - κr 2 2 t 0 L(t -τ, r) 2C κ (1 + χτ ) + δ 0 (τ ) dτ + O(κ 2 ), κ → 0.
After some algebra with use of (70), we obtain

2 mr ∂p -(r, t) ∂r = -δ 0 (t) + χ t 0 e (t-τ )χ δ 0 (τ )dτ (60) + 1 2 Cr 2 e -χt + κr 2 4 C 3 r 2 e -χt (χt -1) + δ 0 (t) -χ t 0 e (t-τ )χ 2 -χ(t -τ ) δ 0 (τ )dτ + O(κ 2 ).
In the interval a 0 < r < a(t), the asymptotic representation of the pressure should be found from Eq. ( 50). After some algebra, we arrive at formula (18) with R(r, t, κ) = R 1 (r, t, κ) + R 2 (r, t, κ) and

R 1 (r, t, κ) = C r a 0 a(t) -ξ (r 2 -ξ 2 )e -χ(t-b(ξ)) × (61) 2ξ + 1 2 κ χ(t -b(ξ)) -1 ξ(r 2 -ξ 2 ) dξ, R 2 (r, t, κ) = χ b(r) 0 a(t) -a(τ ) (r 2 -a 2 (τ ))e -χ(t-τ ) Ca 2 (τ ) -δ 0 (τ ) × 1 - κ 4 (r 2 -a 2 (τ )) 2 -χ(t -τ ) dτ. ( 62
)
It is evident that R(a 0 , t, κ) = 0 and R r (a 0 , t, κ) = 0. Thus, the matching conditions (49) are valid in the main terms for the asymptotic expansion. It is also important to note that substituting (17), ( 18) into ( 14) and collecting the leading order terms, after some algebra, one can recover Eq. ( 100). The following two identities are of great help:

χ a(t) a 0 ξ n e -χ(t-b(ξ)) b ′ (ξ)dξ = a n (t) -K -1 a n (t) (63) n a(t) a 0 ξ n-1 e -χ(t-b(ξ)) dξ = K -1 a n (t) -a n 0 e -χt . ( 64
)
6 Numerical results and discussions 6.1 General remarks.

In order to illustrate the constructed solution of the refined biomechanical contact problem and compare it with that from (Wu et al., 1997a), we will adopt for the computations the same material properties and geometrical parameters of a typical human cartilage reported in [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a). Namely we assume that H A = 0.5 MPa, µ s = 0.25 MPa,

ν = 0.0, H = H A /µ s = 2(1 -ν)/(1 -2ν) = 2.0, k = 2 • 10 -3 mm 4 N -1 s -1 , h = 1 mm, R 0 = 400 mm.
In the sequel we also consider some variation of this set to analyze the sensitivity of the solution to the variation of material characteristics. In all the numerical tests, a constant contact load F (t) = F 0 = 100 N is applied.

Before to present the numerical results themselves, we would like to shortly discuss the computational models adapted here. First, for the reason of comparison, we compute (a) the classic solution reported in [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a) represented here by substituting κ = 0 in Eqs. ( 13), (25) (to compute the radius of the contact zone, a(t), and the punch indentation, δ 0 (t)), alongside with Eqs. ( 17) and ( 18) and the boundary condition (10) to reconstruct the pressure distribution, p(r, t). Note here that the last formulae were presented in (Wu et al., 1997a) in an implicit form.

Then we compute the complete numerical solution in general case (κ > 0) according to formulas (13), ( 14), (41), and (50) together with the boundary condition (10). We utilize here an iterative procedure. All the values of interest: the radius of the contact zone, a(t), the punch indentation, δ 0 (t), and the contact pressure, p(r, t); are computed then in the same common scheme. The procedure appears to be highly unstable and sensitive to variations of the parameters under consideration. We explain this phenomena later. Moreover, since all the parameters in question should be computed at each iteration step, the procedure large amount of memory in comparison with other computational procedures.

Finally, the solution we employ here is the approximation of the previous one with an accuracy of O(κ 2 ), and the solution that is constructed by means of the exact formula (13) and the approximate ones ( 16), ( 17) and ( 18) together with the boundary condition (10). This solution has all advantages of the classic solution described above (Wu et al., 1997a) and simultaneously involves non zero values of the parameter κ, or, in other words, accounts for the influence of the tangential displacements in the model.

Two crucial questions should be answered now. First, whether one should take into account the tangential displacements in the cartilage layer model. And, second, if the answer on the first question is positive, whether one can use for this reason the approximate model instead of the exact one.

To answer the first question, we observe that the small parameter ε enters the right-hand sides of formulas ( 3) and ( 4) as ε and ε 2 , respectively. Thus, the radial tangential displacement becomes more asymptotically important with respect to the vertical normal displacement as ε → 0. However, real biomechanical contact problems should be investigated for finite values of their parameters.

Now, let us rewrite the governing integral equation ( 8) in the dimensionless form. Denoting a characteristic measure of the contact radius by ā0 and employing formulas (5), we obtain

1 r ′ ∂ ∂r ′ r ′ ∂p ′ (r ′ , t ′ ) ∂r ′ + 1 H t ′ 0 1 r ′ ∂ ∂r ′ r ′ ∂p ′ (r ′ , τ ) ∂r ′ dτ + ǫr ′ ∂p ′ (r ′ , t ′ ) ∂r ′ = Cā 0 r ′2 -δ ′ 0 (t ′ ). ( 65
)
We introduced the following standard normalised notation [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a)

p ′ (r ′ , t ′ ) = 1 mā 3 0 p(r, t), δ ′ 0 (t ′ ) = 1 ā0 δ 0 (t), ǫ = ā2 0 κ = 3ā 2 0 2hR 0 . (66)
Here we have taken into account the second formula (9). Finally, taking into account the notation ε = h/ā 0 , we obtain

ǫ = 3 2ε ā0 R 0 . ( 67
)
Thus, formula (67) shows that the new dimensionless parameter ǫ turns out to be inversely proportional to the parameter ε. At the same time, the value of ǫ is proportional to the ratio ā0 /R 0 , which is usually small in real biomechanical contacts. In particular, in the considered example, one can put ā0 = 10 mm [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a), obtaining ε = 0.1 and ǫ = 0.375, what cannot be considered as negligible value. Moreover, the ratio ā0 /R 0 can range in sufficiently wide limits in the experiments on evaluation of material properties by means of indentation. On the other hand, the fact that ǫ = 0.375 does not necessary mean that this should influence greatly all the parameters of the biomechanical contact under consideration.

It is clear that if ǫ tends to zero, the obtained exact solution (as well as the approximate one) approaches the solution given by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF], Wu et al. (1997). It is interesting to note that we cannot calculate exactly the value of ǫ unless we solve the refined contact problem, since this parameter depends on the characteristic measure of the contact radius, ā0 , which is not given before hand. However, one can estimate the value of the parameter ǫ using the solution of the original contact problem. Using formulas (10), ( 27), and (110), we obtain for the case of the external load under consideration the following relation:

3a 2 (t) 2hR 0 ≃ 3.25 F 0 µ s R 2 0 1/3 1 + 3µ s k h 2 t 1/3 . ( 68
)
For the material properties and geometrical parameters of a typical human cartilage reported in [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997), the second factor on the right-hand side of formula ( 68) is less than 1.1 at the moment of time t = 200 s. From formula (68) it is obvious what physical circumstances correspond to diminishing the influence of the tangential displacements effect.

Note also that comparing ( 65) with ( 8) one can rewrite all the results obtained above in the dimensionless form that is more convenient for numerical computations.

Numerical comparison with the classic solution.

Let us first compute and compare the initial contact zone a 0 = a(0) using the aforementioned models. This gives us a

(cl) 0 ≃ 10.8481mm, a (ex) 0 ≃ 10.8494mm, a (ap1) 0 ≃ 10.8481mm, a (ap2) 0 ≃ 10.8496mm.
Here the first two values were obtained from the first two aforementioned models, while the next two correspond to the same approximate model where a (ap1) 0

is the consequence of Eq. ( 16) which is an approximation of the exact solution with accuracy of O(ǫ 2 ) (in the normilised form), whereas a (ap2) 0 was computed from Eq. ( 25) (or equivalently ( 79)) together with (57) that gives the next order of approximation O(ǫ 3 ) as has been mentioned in Section 5.1.

From glance these computations may look surprising and need some explanations. First, according to the equation ( 59) we should expect that the value a (ap1) 0 computed from the first order approximate equation ( 16) gives the same result as that in (Wu et al. (1997a)). This means that, apart of the fact that for arbitrary time the solution differs only on term O(ǫ), at the initial moment the estimate becomes O(ǫ 2 ) instead. This, in fact, has been accurately proven in (52). However, for the value ǫ = 0.375 even this argument could not completely explain the extremely small difference between the obtained numerical results. And only if one looks at all equations related to the asymptotic solution with respect to the parameter κ (see, for example, ( 51) -( 52), (17), and ( 18)), it becomes clear that the value κa 2 0 /4 = ǫ/4 appears there instead of ǫ in the governing equation ( 65).

Other useful observation concerns the contact zone a (cl) (t) computed according to the classic solution [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a), which always gives the lower estimate for the accurate solution a (ex) (t) while that computed according to the second order approximation, a (ap2) (t), always gives the upper estimate. Moreover, the difference between the latter solution and the exact one is negligible at least within reasonable range of parameters appearing in biomechanical contact problem. For this reason we have used this specific approximate model (a (ap2) (t) ∼ = a (ex) (t)) for further computations. Let us remind that the exact solution requires disproportional numerical efforts. Possible errors caused by such replacement are of the order O(ǫ 2 ) according to our estimate or even less as our computations indicated.

In Figs 2345, we present a comparison of the results obtained for the radius of the contact zone, a(t), punch indentation, δ 0 (t), contact pressure, p(r, t), and its derivative with respect to the radial coordinate r, for three stages of deformation corresponding to the moments t = 0, 100, 200 [s]. Computations have been performed for different sets of material and geometrical parameters reported earlier (see Table 1). All curves computed according to the classic model are presented by solid lines whereas the corrected results are given by dashed lines.

Please include Figure 2 somewhere here.

It is clear from the results presented in Fig. 2 that the difference between the obtained solution to the refined problem, which takes into account the radial tangential displacements, and the classic one [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a) increases with time. However, this difference is very small. For example, for the radius of the contact zone, a(t), the discrepancy between the solutions leads to errors of 0.014% for t = 0 s and 0.074% for t = 200 s. However, accounting for the radial tangential displacements leads to the higher impact upon the values of the two most important for application parameters: the punch indentation, δ 0 (t), and the contact pressure, p(r, t).

Please include Figure 3 somewhere here.

In fact, the punch indentation, δ 0 (t), differs for these two models under investigation at the beginning and the end of the indentation process like 3.57% and 3.13%, respectively. Note here that we depicted on the figures the doubled punch indentation, 2δ 0 (t), for reason of the direct comparison with results reported in (Wu et al., 1997a). As concerns the pressure, the difference takes its maximum value at the central point r = 0. For the set of the parameters under consideration we computed the maximum difference of contact pressures as 0.85% and 3.57% for the moments t = 0 s and t = 200 s, respectively. As it has been mentioned above, the refined contact model always yields higher contact pressure maximum.

Please include Figure 4 somewhere here.

Summarizing, for the set of parameters discussed above, the difference between the classic solution and the refined one is not essential, but for the punch indentation and the contact pressure should be important in some cases though. These cases may include inverse problems of determining biomechanical and the geometrical parameters of the cartilage layer from experiments. This situation may change even drastically for other set of biomechanical parameters as they essentially vary with progress of osteoarthritis [START_REF] Wu | Joint contact mechanics in the early stages of osteoarthitis[END_REF]. To illustrate this circumstance, we choose two other sets of parameters corresponding to a "diseased" cartilage. The second set differs from that discussed above only by one parameter -the cartilage permeability, k. Namely, following [START_REF] Wu | Joint contact mechanics in the early stages of osteoarthitis[END_REF], we adopt now the bigger value k ′ = 4k. The corresponding numerical results are presented in Fig. 3. All other notation used remains the same.

As in the previous case, the computed results for the radius of the contact zone, a(t), are close enough with the difference estimated as 0.014% for t = 0 s and 0.595% for t = 200 s. Although at the end of the indentation process the difference even for this parameter becomes pronounced.

Please include Figure 5 somewhere here.

The punch indentation deviates from the classic solution with practically the same accuracy 3.57% and 3.24% at the beginning and the end of the computations. Finally, the contact pressure, p(0, t), differs rather essential at the end of the indentation process. The respective values are 0.82% and 4.73%.

Then we consider the third set of biomechanical and geometrical parameters from the table 1. Now three parameters changed their values [START_REF] Wu | Joint contact mechanics in the early stages of osteoarthitis[END_REF]: the thickness of the cartilage layer h ′ = 2h, the permeability k ′ = 2k and the shear modulus µ ′ s = 0.5µ s , while all other parameters are the same as in [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a).

The corresponding results presented in Fig. 4 are in the same convention as above. In this particular case, the trend is different. Indeed, in contrast to the previous two cases, the difference between the radii of the contact zone stays practically constant with time (varies from 0.0215% to 0.0213%).

The punch indentation, δ 0 (t), exhibits similar tendency slightly changing with time while the difference between the classic solution and that obtained in the frame of the refined model slowly decreases (from 4.465% to 4.266%). But the most interesting conclusion concerns the contact pressure, p(r, t). Namely the difference between the models stays practically without changing (0.954% for t = 0 s and 0.962% for t = 200 s). Moreover, contrary to the previous cases, the variation of the pressure itself remains practically the same during the indentation process, whatever the model is used.

Finally, we consider one set of biomechanical and geometrical parameters which is also slightly different with respect to the first one. Namely, now the thickness of the cartilage layer assumed to be smaller, h ′ = h/2, while all other parameters have the same values as in the case of healthy cartilage layer from [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]Wu et al., 1997a).

Corresponding results presented in Fig. 5. The computed results for the radius of the contact zone, a(t), show the maximal error at level 0.6% for t = 200 s that is similar to the second set of parameters discussed above. The same conclusion is valid for the contact pressure with the discrepancy of 5%. However, in this case, the difference between the values of the indentation parameter, 2δ 0 (t), computed according to the classical and refined asymptotic models is much higher. Namely, it now equal to 9%. The last example shows that it is impossible to estimate a priori the maximum difference as we have discussed above from the analytical considerations. On the other hand, we would like to underline this again that the exact solution for the refined contact problem lies somewhere between those two computed by the asymptotical models. In other words, the difference between the solutions gives the aforementioned estimate.

Main conclusions

The present study results in the analytical solution to the refined unilateral axisymmetric contact problem for a biphasic cartilage layer. It has been shown that the approximate solution (represented by formulas ( 13), ( 58), (17), and ( 18)) is of the same complexity as the classic one constructed by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]; Wu et al. (1997a). Moreover, the asymptotical models not only allow to avoid computations based on the accurate exact result complicated though, but also to determine an accuracy of the approximate computations.

The results of numerical simulations have been summarized in Table 1. The refined solution always leads to some increase in the radius of the contact zone, the punch indentation as well as in the contact pressure at the center of the contact zone. This, in turns, according to the constructed solution, means that the gradient of the contact pressure is also higher than that in the classic model. Utilization of this refined solution in analysis of contact problems for a biphasic cartilage layer is desirable if one investigates different stages of osteoarthritis when biomechanical and geometrical parameters vary rather significantly. However, it may be even more important to take account of the refined contact problem formulation in solving an inverse problem of determining biomechanical and geometrical parameters of cartilage layers from experiments.

It should be noticed that the limitations of the refined formulation of contact problem to biphasic isotropic materials, axisymmetric configurations, and monotonic loadings is something that can be dealt with using numerical models. The developed approach can be generalized for nonaxisymmetric configurations as well as for non-monotonic loadings. However, the introduction of a new material model for articular cartilage (see the rather extensive review by [START_REF] Wilson | The role of computational models in the search for the mechanical behaviour and damage mechanisms of articular cartilage[END_REF]) essentially changes the governing integral equation, and to solve it analytically is a real challenge.

Finally, we stress that that from the geometrical point of view, the formulation of the contact problem is too simplified. However, the simplified geometry does not impose an obstacle for studying the effect of tangential displacements. On the other hand, to estimate this effect by means of FEM simulations is more difficult since we do not know a priori the governing nondimensional parameter controlling this effect. The employed simplifying hypotheses allow us to obtain relatively simple analytical results that can be used for estimating the effect of tangential displacements in solving the contact problems with geometries closer to the reality by means of a finite element method.

Note that this asymptotics is uniformly valid for any positive t.

To estimate behaviour of the kernel L(t, r) near the point t = 0, we rewrite Eq. ( 37) in a slightly different form than Eq. ( 69):

L(t, r) = 1 2πi i∞ -i∞ e G(s,r) -1 -G(s, r) - 1 2 G 2 (s, r) e st -1 ds + 1 8 κχr 2 4 + κχr 2 t e -χt 1 2πi i∞ -i∞ e G(s,r) -1 -G(s, r) - 1 2 G 2 (s, r) ds.
As it follows from the residuum theorem, the last integral is identically equal to zero, while the first integral gives the worst estimated term

L(t, r) = 1 2 κχr 2 + O(t), t → 0, (71)
which is consistent with (70).

Another possibility to effectively compute the kernel function L(t, r) is to deform the integration line in (37) into a unit circle with a center located at the singular point. This gives:

L(t, r) = 1 2π e -χt π -π e κχr 2
2 e -iφ -1 e e iφ t+iφ dφ.

(72)

A2. Kernel function M(t, r). Evaluation and properties

Apart of the fact that there is the simple relationship (40) between the kernel functions L(t, r) and M (t, r) the last one should be computed from (39) or with the help of following equivalent formula:

M (t, r) = 1 2πi i∞-η -i∞-η e G(s,r) -1 -G(s, r) - 1 2 G 2 (s, r) e st ds s - 1 8 κr 2 4 + κr 2 (1 + χt) e -χt . (73)
Here, 0 < η < χ. As κ → 0, the latter formula immediately gives the following asymptotic expansion:

M (t, r) = - κr 2 2 e -χt 1 + κ r 2 4 (1 + χt) + O(κ 3 r 6 ). ( 74
)
The asymptotic formula (74) is valid for any time separated from zero (i. e., 0 < t 0 ≤ t). Moreover, this asymptotics holds true also when r → 0. To evaluate the asymptotic bahaviour of the function M (t, r) as t → 0, we rewrite (73) in the form

M (t, r) = 1 2πi i∞-η -i∞-η e G(s,r) -1 -G(s, r) e st -1 ds s - 1 2 κr 2 e -χt + 1 2πi i∞-η -i∞-η e G(s,r) -1 -G(s, r) ds s ,
where 0 < η < χ. The first integral can be easily estimated to give O(t), while the second one can be evaluated in the closed form using the residuum theorem at point s = 0:

M (t, r) = 1 -e κr 2 2 + O(t), t → 0, (75)
Another representation of the kernel function M (t, r) useful for practical computation is

M (t, r) = e -χt 2π π -π
e κr 2 e -iφ -1 e χte iφ +2iφ 2 dφ e iφ -2 .

(76)

Finally, the following identity for the operators involving the kernel functions L(t, r) and M (t, r) can be established:

q(r, t) + t 0 L(t -τ, r)q(r, τ )dτ = e -κ r 2 2 q(r, t) (77)
+ M (t, r)q(r, 0) + t 0 M (t -τ, r)q(r, τ )dτ.

A3. Recursive relations for the integral characteristics of the contact pressure density

Let us introduce the generalized polar moments of inertia of the contact pressure by means of the formula

F n (t) = 2π a(t) 0 p(ρ, t)ρ 2n+1 dρ. (78)
Thus, Eq. ( 25) can be written in the form

πm 48 Ca 6 (t) = KF 0 (t) + κ a 2 (t) 4 F 0 (t) -F 1 (t) . (79)
Now, taking into account formula (23), we represent Eq. ( 22) in the following form:

Kp(r, t) + κ      - a(t) r p(ρ, t)ρ dρ + 2 a(t) 0 p(ρ, t)ρ ln a(t) ρ dρ -2 r 0 p(ρ, t)ρ ln r ρ dρ      = mC 16 a(t) 2 -r 2 2 + a(t) 2 -r 2 κF 0 (t) 2πa 2 (t) . (80)
The following formulas can be verified by direct calculations using an obvious change of the integration variable: 2π

Now, multiplying the both sides of Eq. ( 80) by 2πr 2n+1 and integrating with respect to r, we obtain the following equation in view of formulas ( 81) -( 83) and the notation (15):

F n (t) = πmC 8(n + 1)(n + 2)(n + 3) K -1 a 2n+6 (t) (84) +κK -1 - a 2n+2 (t) 2(n + 1) 2 (n + 2) F 0 (t) + (n + 2) 2(n + 1) 2 F n+1 (t) .

A4. Short-time asymptotics

Now we can analyze two important limit cases. First, suppose that t → 0. Then we can immediately conclude from the behaviour of the kernel of the integral operator (37) that

κ m ∂p -(r, 0) ∂r = Cr + exp -κr 2 2 -1 r 2C κ + δ 0 (0) . (85)
It is important to check whether this solution satisfy also condition (11). Substituting r = a 0 ≡ a(0), we have the following identity to be valid:

Ca 2 0 = 1 -exp - κa 2 0 2 2C κ + δ 0 (0) . (86)
Next, we determine the distribution of the pressure along the contact area at t = 0 integrating Eq. ( 85). Taking into account the boundary condition (10), we find

κ m p -(r, 0) = C 2 r 2 -a 2 0 + 2C κ + δ 0 (0) g(r, a 0 , κ), (87)
where the following new notation has been introduced:

g(r, b, κ) = r b exp -κξ 2 2 -1 ξ dξ. (88)
Multiplying the both parts of the latter equation by r and integrating within the interval (0, a 0 ), we obtain

κ πm F 0 = - C 4 a 4 0 + 2C κ + δ 0 (0) exp - κa 2 0 2 -1 κ + a 2 0 2 . (89)
Note that Eqs. ( 86) and (89) constitute a system of nonlinear equations allowing to determine two parameters a 0 = a(0) and δ 0 (0).

Below we will show that this approach is equivalent to that discussed in Sections 2 and 3. In fact, eliminating the exponential term from Eqs. ( 86) and ( 89) we obtain

δ 0 (0) = C a 2 0 2 + 2κ m F 0 πa 2 0 . (90)
Notice that Eq. ( 90) coincides with (13) for t = 0. Substituting now this equation into any of Eqs. ( 86) or ( 89), we get

Φ(a 0 ) ≡ Ca 2 0 exp - κa 2 0 2 -1 + 2C κ + C a 2 0 2 + 2κ m F 0 πa 2 0 = 0. ( 91
)
We have obtained here the transcendental equation to determine unknown value a 0 = a(0). One can directly show that Φ(ξ) → +∞ as ξ → +0, Φ(ξ) → -∞ as ξ → +∞ and Φ ′ (ξ) < 0 for any ξ > 0. This means that there exists the unique solution a 0 > 0 to equation ( 91).

On the other hand, one can determine the value of a 0 using a formally different routine. Namely, from Eq. ( 25) at t = 0 it follows that

πm 48 Ca 6 0 = F 0 + κ a 2 0 4 F 0 -2π a 0 0 p(ρ, 0)ρ 3 dρ , (92)
In view of (87) after computation of the integrals, we get

κ m a 0 0 p -(ρ, 0)ρ 3 dρ = -C a 6 0 24 + 2C κ + δ 0 (0 + ) a 4 0 16 + a 2 0 4κ e -κ a 2 0 2 - 1 2κ 2 1 -e -κ a 2 0 2 . (93)
Finally, after some algebra, with taking into account Eqs. ( 86) and ( 89), one can check that ( 92) and ( 93) lead to Eq. ( 91).

If we consider the limit κ → 0, all three conditions will degenerate to the only one

δ 0 (0) = C 2 a 2 0 ( 94
)
coinciding with that obtained in [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]; Wu et al. (1997a). Passing to the same limit as κ → 0 in formula (85), we obtain

p -(r, 0) = m δ 0 (0) a 2 0 -r 2 4 -C a 4 0 -r 4 16 .
Now, taking into account (94), we rewrite the latter equation in the form

p -(r, 0) = mδ 0 (0) 8a 2 0 (a 2 0 -r 2 ) 2 ≡ m 16 C(a 2 0 -r 2 ) 2 (95)
that again corresponds to the result obtained previously by [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF].

A5. Long-time asymptotics

Finally, in order to prove that we have constructed the uniformly valid asymptotic solution for the contact radius, let us verify that the solution of the approximate equation ( 16) satisfies the asymptotic condition (53). In order to simplify the considerations, we assume that F (t) = F 0 , where F 0 = const. Then, formula (53) takes the form

a(t) ∼ c t 1/6 , t → ∞. (96)
Here and afterwards, c denotes the constant which does not depend on t.

Consequently, the following asymptotic conditions must hold

KF (t) ∼ c t, a 2 (t)F (t) ∼ c t 1/3 , t → ∞. (97)
and one can immediately conclude that the first term in the right-hand side of Eq. ( 16) is the dominant one as t → ∞ under the asymptotic condition (96).

Further, integrating by parts, we obtain

K -1 g(t) ∼ t 0 g ′ (τ )e -χ(t-τ ) dτ, t → ∞. (98)
Hence, taking account of the asymptotic condition (96), we find

K -1 a 8 (t) ∼ c t 0 τ 1/3 e -χ(t-τ ) dτ, t → ∞. (99)
Integrating by parts in the integral in (99), we get t 0 τ 1/3 e -χ(t-τ ) dτ = t 1/3 χ -e -χt 3χ t 0 e χτ τ 2/3 dτ.

The last integral can be estimated as follows: Finally, taking into account both relations, we obtain K -1 a 8 (t) ∼ ct 1/3 , t → ∞.

(100)

The asymptotic formula (100) implies that the first term in the round brackets ( 14) is a dominant one as t → ∞ under the asymptotic condition (96). This proves that the approximate equation ( 16) presents the uniformly valid asymptotic approximation to the equation for the radius of the contact area (25) in the whole time range interval t ∈ [0, +∞). By an argument similar to the one above we can consider the general case of monotonic loading when the function F (t) increases with time.

The same line of the reasoning can be provided for the equation ( 58). The behavior at infinity gives the same agreement whereas condition at the initial point, t = 0, reads now πmC 48F (0) a 6 0 = 240 + 60κa 2 0 + 10κ 2 a 4 0 240 + 60κa 2 0 + 9κ 2 a 4 0 .

(101)

Together with the conclusions from Section 4.3, this an additional explanation why the first order approximation (59) has not contained a term of the order O(κ). Notice that the rational approximation (101) coincides with that reported in (52) with accuracy to O(κ 3 ). Figure 2. Prediction of the radius of the contact zone, a(t), indentation 2δ 0 (t) as functions of time and the distribution of the contact pressure p(r, t) and its derivative with respect to the radial coordinate r for three different moments of time t = 0, 100, 200 [s] according to the model reported by Wu et al. (1997a) and the present one for the set of parameters according to [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the first set of parameters in Table 1). , 200 [s] according to the model reported by Wu et al. (1997a) and the present one for the set of parameters: k ′ = 4k and others according to [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the second set of parameters in Table 1). . Prediction of the radius of the contact zone, a(t), indentation 2δ 0 (t) as functions of time and the distribution of the contact pressure p(r, t) and its derivative with respect to the radial coordinate r for three different moments of time t = 0, 100, 200 [s] according to the model reported by Wu et al. (1997a) and and the present one for the set of parameters: h ′ = 2h, k ′ = 2k, µ ′ s = 0.5µ s , while others according to [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the third set of parameters in Table 1). , 200 [s] according to the model reported by Wu et al. (1997a) and and the present one for the set of parameters: h ′ = h/2 and others according to [START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the fourth set of parameters in Table 1).

  = 2e + e χt .

Figure 3 .

 3 Figure3. Prediction of the radius of the contact zone, a(t), indentation 2δ 0 (t) as functions of time and the distribution of the contact pressure p(r, t) and its derivative with respect to the radial coordinate r for three different moments of time t = 0, 100, 200 [s] according to the model reported byWu et al. (1997a) and the present one for the set of parameters: k ′ = 4k and others according to[START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the second set of parameters in Table1).

  Figure4. Prediction of the radius of the contact zone, a(t), indentation 2δ 0 (t) as functions of time and the distribution of the contact pressure p(r, t) and its derivative with respect to the radial coordinate r for three different moments of time t = 0, 100, 200 [s] according to the model reported byWu et al. (1997a) and and the present one for the set of parameters: h ′ = 2h, k ′ = 2k, µ ′ s = 0.5µ s , while others according to[START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the third set of parameters in Table1).

Figure 5 .

 5 Figure5. Prediction of the radius of the contact zone, a(t), indentation 2δ 0 (t) as functions of time and the distribution of the contact pressure p(r, t) and its derivative with respect to the radial coordinate r for three different moments of time t = 0, 100, 200 [s] according to the model reported byWu et al. (1997a) and and the present one for the set of parameters: h ′ = h/2 and others according to[START_REF] Ateshian | An asymptotic solution for the contact of two biphasic cartilage layers[END_REF]. (See the fourth set of parameters in Table1).

Table 1

 1 Maximum difference between the two models on the time interval t ∈ [0,200 s]

	µ s	k	h	Error in Error in Error in
	MPa mm 4 N -1 s -1 mm	a(t)	p(r, t)	δ 0 (t)
	0.25	2 • 10 -3	1	0.07%%	3.57%	3.57%
	0.25	8 • 10 -3	1	0.59%	4.73%	3.57%
	0.125	4 • 10 -3	2	0.02%	0.96%	4.46%
	0.25	2 • 10 -3	0.5	0.6%	5%	9%

Mishuris, G. and Argatov, I. Exact solution to a refined contact problem for biphasic cartilage layers, Proc. 1st International Conference on Mathematical and Computational Biomedical Engineering -CMBE2009

June 29 -July 1, 2009, Swansea, UK. P. Nithiarasu and R. Löhner (Eds.). P. 1-4

Acknowledgment

The authors would like to express their gratitude to the Referees for their helpful comments and discussions.

Appendix

A1. Kernel function L(t, r). Evaluation and properties

At first glance the function L(t, r) can be easily computed from (37) as the integrand tends to zero at infinity under properly chosen value of the parameter η. Moreover, it has only one singular point s = -χ in the entire complex domain. However, this point is not a pole but an essential singularity, and an attempt to compute the corresponding residuum leads to an infinite series that is not the most convenient way for computations.

On the other hand, the integral (37) can be easily computed numerically (if η = 0 is chosen, it gives, in fact, the standard Fourier transform). To improve the quality of the computation one can rewrite the formula in the equivalent form

Moreover, instead of the simple choice η = 0 in (37) and ( 69), it is always better to take in computations a negative value of the parameter η. For example, we take η = -χ/2.

For some limiting values of the problem parameters we can work out the asymptotics of the kernel L(t, r). Thus, using the reasoning like in (69) we can derive the asymptotic expansion as r → 0 for any fixed value t separated from zero (i. e., t ≥ t 0 > 0):

Complete set of figures (1 -5)

Figure 1. Support of the solution of Eq. ( 27).

t r a 0 = a(0) r = a(t) Ω + Ω -