Analysis of the Unilateral Contact Problem for Biphasic Cartilage Layers with an Elliptic Contact Zone and Accounting for Tangential Displacements

Analysis of the Unilateral Contact Problem for Biphasic Cartilage Layers with an Elliptic Contact


Introduction
Biomechanical contact problems involving transmission of forces across biological joints are of considerable practical interest (see, e.g.[2,3,11,14]).Many analytical solutions to the problem of contact interaction of articular cartilage surfaces in joints are available.In particular, Ateshian et al. [8] obtained an asymptotic solution for the axisymmetric contact problem for two identical biphasic cartilage layers consisting of a solid phase and a fluid phase and attached to two rigid impermeable spherical bones of equal radii.Later, Wu et al. [16] extended this solution to a more general axisymmetric model by combining the assumption of the kinetic relationship from classical contact mechanics [12] with the joint contact model [8] for the contact of two biphasic cartilage layers.An improved solution for the contact of two biphasic cartilage layers in the axisymmetric setting, which can be used for dynamic loading, was obtained by Wu et al. [15].
An asymptotic modeling approach to study the contact problem for biphasic cartilage layers has been performed by Argatov and Mishuris in a series of articles (see [4,6,7]).In particular, it was shown [6] that accounting for the tangential displacements is important in the case of diseased cartilage where the measurement of indentation depth may differ even as much as 10% in comparison with the healthy case.In [4], the unilateral contact problem for articular cartilages bonded to subchondral bones with a contact zone in the shape of an arbitrary ellipse has been considered, and a closed form analytic solution was found.Exploiting this exact result, Argatov and Mishuris [7] have performed perturbation analysis of the contact problem with approximate geometry of the contact surfaces.Other analytic solutions for the contact problem were found using the viscoelastic cartilage model for elliptic contact zone in [5].A new methodology for modeling articular tibio-femoral contact based on the developed asymptotic model of frictionless elliptical contact interaction between thin biphasic cartilage layers was presented in [2].The mathematical model of articular contact was extended to the case of contact between arbitrary viscoelastic incompressible coating layers.
The constitutive model for biphasic cartilage layers has been extensively discussed in the literature.Our formulation most closely resembles the model proposed by Ateshyan et al [8].We omit a detailed description of the modelling due to a lack of space.Instead, we restrict the discussion, by appropriate citation, to the basic model, with clear identification of the origins of the asymptotic model.
The principal originality of this work, with contrast to papers [6] and [4], is in the accounting for tangential displacements in the contact problem for cartilage layers while using a contact zone of elliptical shape, based on the biphasic model.Although the load is normal, the displacements of the material points on the contact zone have both normal and tangential components, since the surface of the bone is not flat.Despite an absence of friction, the tangential displacement is small but present, and perhaps essential, as has been shown in contact mechanics (with reference to the book by Johnson [12]).Comparing our results with those of other authors we come to the conclusion that accounting for the tangential displacements is important in determining a more accurate approximation of the real behaviour of the complex "bone-cartilage".
Note that the perturbation method proposed in [7] could be one of the options for the analysis, however, the procedure is too complex to perform even a few asymptotic steps.Here, employing some technique and ideas from [6] and [4], we propose another way to construct the asymptotics which utilizes the assumption that the shape of the contact zone is an ellipse at the initial stage of deformation and can be regarded as a small perturbation of the ellipse at any other stage of deformation.
The paper is organized as follows.The unilateral contact problem formulation and its linearization are presented in Section 2, where a special case of the contact configuration with one cartilage layer being plane and rigid is also considered in detail.In Section 3, we derive exact general relationships between the contact approach and some integral characteristics of the contact pressure, including the contact force.In Section 3.3, we obtain asymptotic representations for the contact pressure integral characteristics in terms of the contact approach and some integral characteristics of the contact zone.The zero-order and first-order asymptotic approximations for the solution to the contact problem are obtained in Subsections 4.1 and 4.2, respectively.Detailed calculations which led to the corresponding sets of equations are presented in [13].The first-order approximation problem constitutes the main result of the present study.Section 5 presents a numerical analysis of the model.On the basis of this discussion of the obtained numerical results we make some conclusions concerning the model.

Formulation of the contact problem
We consider a frictionless contact between two thin linear biphasic cartilage layers firmly attached to rigid bones shaped like elliptic paraboloids (see Fig. 1).It is a common assumption in most of papers devoted to the study of the bone-cartilage system to consider the bone as a rigid elliptic paraboloid.Since the stiffness of the bone is obviously much greater than that of the cartilage, this assumption seems physically consistent.The geometrical assumptions are a common simplification in the literature, allowing us to: a) analytically identify basic features of the contact problem, b) consider the solution as a benchmark for FEM simulations.
In the Cartesian co-ordinates (x 1 , x 2 , z) = (x, z) the equations for the two cartilage surfaces can be written in the form z = (−1) n Φ (n) (x), n = 1, 2, where being the curvature radii of the n-th bone surface at its apex.
We note that assuming the bone to have the form of an elliptical paraboloid is practically reasonable for approximation of the real shape of human bones (see [4] and references therein).
In the undeformed state, the cartilage-bone systems occupy convex domains z ≤ −Φ (1) (x) and z ≥ Φ (2) (x), respectively.They are in the initial contact with the plane z = 0 at the origin of the co-ordinate system.
We denote by w 1 (x, t), w 2 (x, t) the local vertical displacements of the corresponding cartilage surfaces.Let also u 1 (x, t), u 2 (x, t) be the local horizontal (tangential) displacements of the corresponding surface of the cartilages.Finally, we denote by P (x, t) the contact pressure density.Following [12] the equations for the cartilage surfaces can be written in the following form: Here, δ 1 , δ 2 are some (positive) vertical displacements of the rigid bones.Note also that the vertical displacements w 1 , w 2 are positive, while the tangential displacements u 1 , u 2 are directed outside of the contact zone.More detailed modelling of the vertical and tangential displacements can be found in [12].
Denoting by δ * (t) = δ 1 (t) + δ 2 (t) the contact approach of the bones, we get from (2.2) the following inequality: 3) It was shown in [8] (see also [6]) that the vertical and the tangential displacements of each bone (taking the asymptotic model of the cartilage layer into account) can be represented in the form (2.5) Here n = h n /a 0 are dimensionless small parameters, h 1 , h 2 mean the thicknesses of the cartilage layers, and a 0 denotes a characteristic measure of the contact zone (see the detailed description of the role of this parameter in [6]; the values taken for numerical analysis of the model are given latter in this section), H n = (λ s,n + 2µ s,n )/µ s,n are material parameters of cartilages, where λ s,n and µ s,n represent the first Lame coefficient and the shear modulus of the solid phase of the n-th cartilage tissue.Note that u 1 and u 2 in (2.5) do not necessarily coincide, they depend on both spatial variables x 1 , x 2 , and on the time variable t.
Following [8], we introduce new spatial variables and time variable via formulas where a 0 is a characteristic measure of the contact zone, and k 1 , k 2 are the cartilage's permeabilities.In these variables we have the following relations on the contact area ω(t) encircled by the curve Γ (t) = ∂ω(t) (here and in the following, we retain the same notation for displacements w n , u n and for the contact pressure P ): Further the equality in (2.3), i.e., determines the contact area ω(t).
Now we substitute (2.6), (2.7) into (2.8) and obtain the governing equation relating the contact pressure with the vertical approach of the bones δ * (t) in the following form: Here we have introduced the notation (2.10) Thus, it follows from (2.1) and (2.10) that the functions Φ and Φ are given by Note that the coefficients in A and B are positive dimensionless numbers, which are less than unit.
Without loss of generality, one can assume that A > B. Then, Equation (2.9) can be rewritten in an equivalent form, using all dimensionless parameters: where the following notation has been introduced: ) It is important to note that χ = O(1), µε χ.
Discussion of the characteristic values of the introduced parameters is presented, e.g., in [6,8].We note that in numerical analysis of the model we can take a 0 = b(0) 1 + e 2 1 as the initial value of the characteristic measure of the contact zone.
Since the solution of (2.11) depends on the parameter ε, it is customer to denote an unknown contact pressure by P = P ε in what follows.Note that the problem for ε = 0 coincides with that considered in [4].
Equation (2.11) is the equation for determination of the contact pressure In particular, in the case when the contact domain is represented by an ellipse We supply Equation (2.11) with the following boundary conditions: The equilibrium equation connects the external load F (t), unknown contact pressure P ε (x, t), and unknown contact domain ω ε (t).
Remark 1.The problem (2.11), (2.15), (2.16), (2.17) has the following form with unknown boundaries for the contact domain ∂Ω, an unknown indentation parameter δ and an unknown contact pressure P (where ε is a small parameter, g and f are given functions in Ω, and K is the Volterra operator).In [4] an exact solution was found, corresponding to the case ε = 0 (in our notation), for elliptical contact.In [6], the existence of a solution was proven, under the assumption of an axisymmetric initial configuration of the contact zone (i.e. when g(x, y) = g(r), f (x, y) = f (r)).Thus, existence of the solution in a more general case, for small values of the parameterε = 0 or small eccentricity, follows from the standard results of perturbation analysis of nonlinear boundary value problems for the Laplace equation.

Special case of the contact configuration
In order to check the content of formula (2.9) we consider here a special case, namely, we suppose that the lower part cartilage layer is plane and rigid (the same assumption was employed in [16]), it means that µ s,2 = ∞ and R ≡ 0, Φ ≡ Φ (2) .
In this case we have got the following equation for determination of the contact domain ω(t) in the form similar to (2.9): Here we will have At the same time, small changes have to be made in the right-hand side of Equation (2.18) as follows: .
Thus Equation (2.18) can be rewritten as It can be easily checked that in the axisymmetric case Equation (2.19) reduces to the governing differential equation obtained in [6].
3 A priori estimate of the solution

Estimates of the indentation parameter
In our model we assume that the external load is non-decreasing.Thus, the contact domain is monotonically expanded, i.e.
It is convenient to suppose also that the contact pressure is defined on the whole plane.For this we simply extend the density P ε (x, t) by assuming that Integrating (2.11) over contact domain ω(t), we get For simplicity of notation, we omit here (and everywhere in the next two sections) the subindex ε in ω ε .From the monotonicity of the contact domain (3.1) and assumption (3.2), it follows that the second integral on the left-hand side can be written in the form Using the second Green's formula with u ≡ 1 and v = P ε (x, t) we get the following relation in view of the boundary condition (2.16): Therefore, the both integrals on the left-hand side of (3.3) vanish.
Further, we use the first Green's formula with ψ(x) = Ψ 2 (x) and ϕ(x) = P ε (x, t).In this case the integral on the righthand side vanishes in view of (2.15), and we obtain the relation where we used the equilibrium equation (2.17) and the identity with e 2 being defined in (2.12).
In what follows, it is convenient to have the following notation for the integrals of the product of k-th power of the function Ψ 1 and l-th power of the function Ψ 2 : In particular, A 0,0 (ω) is the area of the contact domain.It is to remember that the constants A k,l (ω) depend finally on t, but we omitted this fact in the notation in order to avoid cumbersome expressions.Computations of A k,l (ω) for the elliptic domain (2.14) are given in [13, Appendix, Sec.6.1].
Taking into account Equations (3.5) and (3.7), we get This formula allows us to compute the contact approach δ ε (t) as a function of the total external force F (t) and the main axes of the ellipse describing the shape of the contact zone, which in fact depends on time too.

Integral identities for the contact pressure
In order to write out a more informative equation for the contact load, we use the following trick.We multiply both sides of (2.11) by the function v(x) = Ψ 2 (x) and integrate the obtained equation over the contact domain ω(t) Let us calculate the integrals in this relation by using Green's formulas.For the first integral on the left-hand side we use formula (3.4) with u = Ψ 2 , v = P ε and the boundary conditions (2.15), (2.16).Hence, we obtain Now taking into account (3.8), we get For the second integral on the left-hand side, we apply the same approach, but interchange first the integrals over ω ε (t) and over τ ∈ (0, t) exploiting the load monotonicity.Therefore, we arrive at the equation (3.12) Math.Model.Anal., 2016, iFirst:1-25.
For the first and second integrals on the right-hand side, we simply use the notation (3.9), which gives Finally, for the third integral on the right-hand side, we make use of the following simple formula which follows immediately from the definition of Ψ 2 : Then we can apply Green's formula (3.6) and the boundary conditions (2.15), (2.16) to find By applying the second Green's formula (3.4) with This integral still contains the unknown density of contact pressure P ε (x, t).
Let us define Here, we have introduced the Volterra operator K as follows: Note that the integral in the right-hand side of the equation (3.15) allows to continue the same procedure to deliver an asymptotic estimate for this equation.
We continue to proceed with Equation (3.15) on the next steps.

Posteriori estimates for the contact pressure
Now we proceed to calculate the last integral in (3.15).For this we multiply the governing integral equation (2.11) by Ψ j 2 (x) (j ≥ 2) and integrate over contact domain ω(t): By using the same argument as on the previous step, we get (3.17) For the last integral we use the relations Therefore, the integral has been obtained as a solution of the integral equation (3.17).Here the inverse operator K −1 is defined by the formula Performing the same computation, we obtain the following representation for the integral in the right-hand side of (3.15): Substituting this representation into Equation (3.15), we finally get The latter relation allows us to determine the problem parameters asymptotically with any prescribed accuracy.
Note that apart from the fact that the shapes of the contacting bones are elliptical paraboloids, no additional assumptions on the shape of the contact zone have been made.On the other hand, no proof was offered to show that the contact zone is approximately represented by an ellipse.This will be done later.
Remark 2. For every t for which the contact pressure P ε (t) is bounded and the contact region ω(t) belongs to a bounded domain, the remainder ε N +1 (N +2)! µ N +1 M (N +2) P ε (t) in formula (3.18) tends to zero as N → ∞.Thus, the series corresponding to the sum on the right hand-side of (3.18) is converging.
4 Asymptotic solution to the contact problem

Zero-order approximation
First, we get solution of the problem for ε = 0.In this case Equation (2.11) has the form where Ψ 1 (x) is defined in (2.12).Since we know from [4] that the contact zone is an ellipse at this stage of approximation we will have Using formula (4.1) and calculations presented in [13, Appendix, Sec.6.1], one can find that and therefore Note that formulas (4.2) and (4.3) contain two known constants e 1 and e 2 defined in (2.12) and two still unknown functions b 0 (t) and β 0 (t), which are the main semi-axis and the eccentricity of the ellipse The leading terms in (3.18) imply (for N = 0) the following equation: Here, K is the Volterra integral operator defined in (3.16).
Analogously, using some results from [13, Appendix, Sec.6.1], we obtain and thus To find the functions b 0 (t) and β 0 (t) together with the pressure distribution over the contact zone, P (0) (x, t), we follow [4] and introduce a new unknown function In the case of monotone external load, this function should satisfy the Poisson equation (following from (2.9)) with the boundary conditions (2.15), (2.16).

It is customary to rewrite this relation in the form
where Bearing in mind that the function Ψ 1 (x) is a quadratic polynomial (compare with (2.12)), it is natural to look for the solution of such problem in the form of a polynomial in x 1 , x 2 of the fourth degree, that is Note that the term in the brackets vanishes on the boundary ω 0 , and thus the condition (2.15) is satisfied automatically.
In [13,Appendix,Sec. 6.2], it has been shown that Q 0 is a polynomial of the second order having the form Taken into account this representation we arrive at the following relations (see [13,Appendix,Sec. 6.3]): , (4.9) This system allows us to determine the unknown functions b 0 (t) and β 0 (t).
Indeed, eliminating η 0 from the last two equations, we get a bi-quadratic equation defining the value of the parameter β 0 , i.e., By definition, β 0 is a positive parameter, thus the unique positive solution of (4.11) has the form Note that at the zero-approximation the parameter β 0 does not depend on time.
The other parameter, η 0 (t), can be computed directly from (4.9) or (4.10), if one knows the remaining constant b 0 (t).Moreover, taking into account (4.8) and ( 4.3), one can use an equivalent formula η 0 (t) = µb 4 0 (β 2 0 + e 2 1 )/16β 2 0 (1 + β 2 0 ).In the same way, one can offer, in addition to (4.3), two equivalent representations for the indentation parameter Finally, the major semi-axis b 0 of the ellipse ω 0 is determined as follows: 1/6 . (4.13) Note that the parameters b 0 , η 0 as well as the indentation, δ 0 , depend on time t in contrast to the ellipse eccentricity β 0 .Now, it remains only to find the pressure over the contact area.Using (4.5) and (4.7), we get If (x 1 , x 2 ) belongs to the initial contact zone, i.e. 1 − If (x 1 , x 2 ) lies outside of the initial contact zone, i.e. 1 − The critical moment of time t * is determined by the formula Using (4.13), we get If the load is stepwise, we have F (t) = F 0 .Hence, we find that Note that in this case F 0 .
This finishes the zero iteration step.Note that the results of this section after changing the notation coincide with those obtained in [6].

First-order approximation problem
For the next steps we consider an appropriately deformed contact domain ω ε , defined as a perturbation of the zero-order one ω 0 .Namely, we assume that it can be written in the form where unknown polynomials are taken in the forms Note that for ε = 0 the solution form coincides with (4.4), if one take The idea behind such choice of the asymptotic anzatz is to satisfy the boundary conditions (2.15) and (2.16) automatically.This will be achieved by putting Now, when the boundary conditions are valid, we will satisfy the governing equation (2.9).Note that P , where p j = K(P j ), j = 0, 1, and Math.Model.Anal., 2016, iFirst:1-25.
Substituting this representation into Equation (2.9), we obtain where the parameter δ (1) ε is represented in the same form as P We can write Equation (4.16) with the accuracy to the terms of O(ε 2 ) as follows: An extended variant of this equation can be written by using the definition of all components of the equation and by comparing coefficients at different powers of x 1 , x 2 , so that (1 + β 2 1 ) = −µδ (1) , ( 4η (1) 4η (1) where In the system (4.17)-(4.22)we have 6 equations and 7 unknowns: ε , b 1 (t), β 1 (t), and a 40 , a 22 , a 04 (coefficients of the polynomial Q 1 ).Therefore, we have to add an extra equation to the above system, namely where F 1 (t) can be represented in the form We also make use of Eq. (3.18) written for this approximation step with the accuracy of O(ε 2 ) in the form Remark 3. Note that putting ε = 0, the system (4.17)-(4.22),(4.24) transforms to the previous case evaluated in the previous section.

Discussion of the proposed asymptotic procedure
First of all, observe that at t = 0, the contact problem for biphasic layers reduces to that for elastic incompressible layers.The contact problem in the latter case were studied in a number of papers [1,9,10,17], however, without taking into account the tangential displacements.
• Having them we can compute the quantity θ 2k,2l (t) from (4.23), • Then, from the system of three equations (4.20)-(4.22)we compute the constants a 40 , a 22 , a 04 assuming the values of η, b, β as above.
• Finally from the system of four equations (4.17We note that formulas (2.4) and (2.5) for the vertical and tangential displacements contain different powers of parameters , namely, 2 and , respectively.Note also that our analysis (with the values of another parameters taken into account) shows, that the role of these magnitudes (vertical and tangential displacements) is quite opposite.In the final equation (see (2.11)) the leading terms, corresponding to the vertical displacement, contain the zero power of the new small parameter ε, but the leading terms, corresponding to the tangential displacements, contain the first power of ε.
An extended discussion of the model is presented in the next section.

Numerical results and conclusions
In this section we present a numerical analysis of the algorithm and a discussion of its fundamental peculiarities.We will then address the main question of this analysis, specifically the importance of accounting for the tangential displacement of the contact problem, without an assumption of axisymmetry.
We also compare our approximation to the other available results.
In the axisymmetric case (see [8]), it is commonly assumed that the human bone is approximated by a paraboloid with curvature radius R = 400 mm.We investigate this case (i.e. with R 1 = R 2 = 400 mm) and also a few other possible cases with curvature radii R 1 = 200, R 2 = 300, R 1 = 350, R 2 = 400 and R 1 = 300, R 2 = 600.Our numerical results are provided for two different cartilages.
They are characterized by the constants (n = 1, 2) For these two different cartilages the thicknesses are taken to be h n = 1 mm and h n = 0.5 mm (where the first thickness corresponds to healthy cartilage).
Finally, the average external load is taken to be F = 100 N and for the maximal time of observation we take t = 200 s.These choices for the parameters are in common with many other papers devoted to the cartilage model (cf., [7], [8], [15]).

Numerical results
Here we analyze the convergence of the proposed scheme only the parameters which characterize our solution, namely, β -the eccentricity of the contact zone, b -its smallest semi-axis, δ -the indentation parameter, and η -the maximum of the function η(t) related to the contact pressure P (see (4.5), (4.6)).With this we take into account the application goal of this paper.
We have estimated the convergence rates of the parameters for all analyzed cartilages but present here in Figures 2, 3 only two distinctive cases: large eccentricity in Figure 2, and small eccentricity in Figure 3.We observe the following features of the algorithm: • it converges more rapidly in the case of larger eccentricity, where even 20 iterations are sufficient to reach the "good" rate; • the slowest case is the circular contact, where the same rate is reached after more that 50 iterations; • the level of the convergence rate for all analyzed parameters (β, b, δ, η) is essentially the same; • the convergence is the most accurate when considering eccentricity β (in comparison with three other parameters b, δ, η); • the worst level of convergence is that found when considering η.
The results for successive rates of convergence for those cartilages not discussed in Figures 2, 3 look similarly.As a result, to guarantee the best convergence we choose to make 50 iterations for further computations.

Comparison of the results in the case of the circular contact zone
Here we compare the results of our algorithm, in the case of a circular contact, with those available in the literature, specifically • Wu et al (1997) [15], where the axisymmetric contact problem was analytically solved without accounting for the tangential displacement; • Argatov-Mishuris (A&M (2010)) [6], where the Wu model was extended to take tangential displacement into consideration and to estimate its impact.
In Figures 4, 5 we present for such a comparison the results from [6], [15] alongside ours (red line).
The following immediate conclusions can be made from these figures: one term asymptotic expansions do not guarantee that a approximate result will be very close to the exact numerical solution of Argatov-Mishuris (2010); the  results of our model are close enough to previous results to be of the same order.
Our calculations have been made by taking only the first term of the asymptotic expansion.If greater accuracy is required for the computed parameters, it is necessary to consider at least two terms of asymptotics.In particular, it can be accomplished using the analytic calculations presented in the paper.
For the purposes of this paper, the above accuracy is sufficient, as will be seen in the next subsection.

Comparison of the present approximate solution for the elliptic contact zone
In this subsection we compare the parameters computed on the basis of our approximate solution with the exact solution presented in [4].We note that the exact result in [4] was obtained for an elliptic contact zone but without accounting for tangential displacement.Since the only one term approximation is not particularly accurate we evaluate further on only average characteristics like eccentricity β of the contact zone and the indentation parameter δ.Four types of cartilages with different eccentricity are analyzed in Figures 6, 7 From Figures 6, 7 we can reach the following conclusions: • the general tendency is the same for all parameters, specifically that the deviation grows with the ratio R 2 /R 1 ; • the maximal relative error (20%) for the eccentricity β is found for the radii R 1 = 300, R 2 = 600 and the minimal (less than 1%) for the circular case; • the indentation parameter δ increases by one order of magnitude for the largest values of the ratio R 2 /R 1 .and taking into account that

Figure 1 .
Figure 1.Schematic representation of the the system bone-cartilage.
Now we rewrite the relation (3.10) by using the results for all integrals (3.11)-(3.14) in the following form:
)-(4.19) and (4.24) considering the right-hand side known (computed by the values know from the previous computations), we found new values η, b, β, δ and compare them with the previous computations.If the required accuracy has achieved we stop the computation, if not we are going to the second step of this iterative procedure.

Figure 2 .
Figure 2. Successive rate of the convergence for the parameters in standard and logarithmic scale.R 1 = 300, R 2 = 600, left -for h = 1, right -for h = 0.5.

Figure 3 .
Figure 3. Successive rate of the convergence for the parameters in standard and logarithmic scale.R 1 = 350, R 2 = 400, left -for h = 1, right -for h = 0.5.

Figure 4 .Figure 5 .
Figure 4. Comparison of the values of the parameter b in different models; left -for h = 1, right -for h = 0.5.
. The respective relative deviations (or relative errors) are given.The leftmost figures correspond to the thickness of the cartilage h = 1, and rightmost figures, to the thickness h = 0.5.

Figure 6 .Figure 7 .. 2 −
Figure 6.Relative error for the parameter β.The base of comparison is our approximate solution; left -for h = 1, right -for h = 0.5.