Analysis of the unilateral contact problem for biphasic cartilage layers with an elliptic contact zone and accounting for the tangential displacements
aa r X i v : . [ phy s i c s . b i o - ph ] M a r Analysis of the unilateral contact problem for biphasiccartilage layers with an elliptic contact zone andaccounting for the tangential displacements
A.A. Koroleva ∗ , S.V. Rogosin † , G.S. Mishuris ‡ March 27, 2018
Abstract:
A three-dimensional unilateral contact problem for articular cartilage layers at-tached to subchondral bones shaped as elliptic paraboloids is considered in the frameworkof the biphasic cartilage model. The main novelty of the study is in accounting not onlyfor the normal (vertical), but also for tangential vertical (horisontal) displacements of thecontacting surfaces. Exact general relationships have been established between the contactapproach and some integral characteristics of the contact pressure, including the contactforce. Asymptotic representations for the contact pressure integral characteristics are ob-tained in terms of the contact approach and some integral characteristics of the contact zone.The main result is represented by the first-order approximation problem.
Biomechanical contact problems involving transmission of forces across biological joints areof considerable practical interest (see, e.g. [2, 3, 11, 13]). Many analytical solutions tothe problem of contact interaction of articular cartilage surfaces in joints are available. Inparticular, Ateshian et al. [8] obtained an asymptotic solution for the axisymmetric contactproblem for two identical biphasic cartilage layers consisting of a solid phase and a fluid ∗ Belarusian State University, Belarus. † Aberystwyth University, UK, and Belarusian State University, Belarus. ‡ Aberystwyth University, UK.The research leading to these results has received funding from the People Programme (Marie CurieActions) of the European Union’s Seventh Framework Programme FP7/2007-2013/ under REA grant agree-ment PIRSES-GA-2013-610547 - TAMER and by FP7-PEOPLE-2012-IAPP through the grant PIAP-GA-2012-284544-PARM2. The authors are thankful to Dr. I. Argatov for important suggestions improving theproposed model.
Formulation of the contact problem
We consider a frictionless contact between two thin linear biphasic cartilage layers firmlyattached to rigid bones shaped like elliptic paraboloids. In the Cartesian co-ordinates( x , x , z ) = ( x , z ) the equations for the two cartilage surfaces can be written in the form z = ( − n Φ ( n ) ( x ), n = 1 ,
2, whereΦ ( n ) ( x ) = x R ( n )1 + x R ( n )2 (2.1)with R ( n )1 , R ( n )2 being the curvature radii of the n -th bone surface at its apex.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 theorigin of the co-ordinate system.We denote by w ( x , t ), w ( x , t ) the local vertical displacements of the correspondingcartilage surfaces. Let also u ( x , t ), u ( x , t ) be the local horizontal (tangential) displacementsof the corresponding surface of the cartilages. Finally, we denote by P ( x , t ) the contactpressure density. In this notation the equations for the cartilage surfaces can be written inthe following form: z = δ ( t ) − Φ (1) ( x + u ( x , t )) + w ( x , t ) ,z = − δ ( t ) + Φ (2) ( x + u ( x , t )) − w ( x , t ) . (2.2)Here, δ , δ are some (positive) vertical displacements of the rigid bones. Note also thatthe vertical displacements w , w are positive, while the tangential displacements u , u aredirected outside of the contact zone. Denoting by δ ∗ ( t ) = δ ( t ) + δ ( t ) the contact approachof the bones, we get from (2.2) the following inequality: δ ∗ ( t ) + w ( x , t ) + w ( x , t ) ≤ Φ (1) ( x + u ( x , t )) + Φ (2) ( x + u ( x , t )) . (2.3)It was shown in [8] (see also [4]) that vertical and the tangential displacements of eachbone can be represented in the form w n ( x , t ) = h n ǫ n µ s,n ∆ P ( x , t ) + 3 H n t Z ∆ P ( x , τ ) dτ , n = 1 , , (2.4) u n ( x , t ) = − h n ǫ n µ s,n ∇ P ( x , t ) , n = 1 , . (2.5)Here ǫ n = h n /a are dimensionless small parameters, h , h mean the thicknesses of thecartilage layers, and a denotes a characteristic measure of the contact zone (see the detaileddescription of the role of this parameter in [4]), H n = ( λ s,n + 2 µ s,n ) /µ s,n are material pa-rameters of cartilages, where λ s,n and µ s,n represent the first Lame coefficient and the shear3odulus of the solid phase of the n -th cartilage tissue. Note that u and u in (2.5) do notnecessarily coincide, they depend on both spatial variables x , x , and on the time variable t . Following [8], we introduce new spatial variables and time variable via formulas x ′ j = x j a , j = 1 , , t ′ = χtµ , where χ = 3 µ s, k h + 3 µ s, k h , µ = µ s, λ s, + 2 µ s, + µ s, λ s, + 2 µ s, ,a is a characteristic measure of the contact zone, and k , k are the cartilage’s permeabilities.In these variables we have the following relations on the contact area ω ( t ) encircled by thecurve Γ( t ) = ∂ω ( t ): w ( x ′ , t ′ ) + w ( x ′ , t ′ ) = (cid:18) h µ s, + h µ s, (cid:19) ∆ P ( x ′ , t ′ ) + χ t ′ Z ∆ P ( x ′ , τ ′ ) dτ ′ , (2.6)Φ ( n ) ( x ′ + u n ( x ′ , t ′ )) ≃ Φ ( n ) ( x ′ ) − h n a µ s,n ∇ Φ ( n ) ( x ′ ) · ∇ P ( x ′ , t ′ ) , n = 1 , . (2.7)Further the equality in (2.3), i.e., δ ∗ ( t ′ ) + w ( x ′ , t ′ ) + w ( x ′ , t ′ ) = Φ (1) ( x + u ( x , t )) + Φ (2) ( x + u ( x , t )) , (2.8)determines the contact area ω ( t ).Now we substitute (2.6), (2.7) into (2.8) and obtain the governing equation relating thecontact pressure with the vertical approach of the bones δ ∗ ( t ) in the following form (fromnow on we keep the names of new unknown functions, e.g. Φ( x ) := Φ( x ′ a ) etc.):∆ P ( x , t ) + χ t Z ∆ P ( x , τ ) dτ = m (cid:16) Φ( x ) − δ ∗ ( t ) − ∇ e Φ( x ) · ∇ P ( x , t ) (cid:17) . (2.9)Here we have introduced the notation m = (cid:18) h µ s, + h µ s, (cid:19) − , (2.10)Φ( x ′ ) = Φ (1) ( x ′ ) + Φ (2) ( x ′ ) . (2.11)Thus, it follows from (2.1) and (2.11) that the functions Φ and e Φ are given byΦ( x ) = Φ( x , x ) = Ax + Bx (2.12)4ith A = 12 R (1)1 + 12 R (2)1 , B = 12 R (1)2 + 12 R (2)2 and e Φ( x ) = e Φ( x , x ) = e Ax + e Bx . (2.13)Note that the coefficients in e A and e B are positive dimensionless numbers, which are lessthan unit.Without loss of generality, one can assume that A > B . Then, Eq. (2.9) can berewritten in an equivalent form, using all dimensionless parameters:∆ P ε ( x , t ) + χ t Z ∆ P ε ( x , τ ) dτ = µ (cid:0) Ψ ( x ) − δ ε ( t ) − ε ∇ Ψ ( x ) · ∇ P ε ( x , t ) (cid:1) , (2.14)where the following notation has been introduced:Ψ j ( x ) = x + e j x , j = 1 , , δ ε ( t ) = 1 A δ ∗ ( t ) , (2.15) µ = Am, e = p B/A, e = q e B/ e A, ε = e AA .
It is important to note that χ = O (1) , µε ≪ χ. (2.16)Discussion of the characteristic values of the introduced parameters is presented in Section5 (see also [8, 4]).Since the solution of (2.14) depends on the parameter ε , it is customer to denote anunknown contact pressure by P = P ε in what follows. Note that the problem for ε = 0coincides with that considered in [5], where an exact solution to this problem was found.Equation (2.14) is the equation for determination of the contact pressure P ε ( x , t ) ≥ x ∈ ω ε ( t ). In particular, in the case when the contact domain is represented by an ellipse ω ε ( t ) = (cid:26) x ∈ R : x b ( t, ε ) + β ( t, ε ) x b ( t, ε ) ≤ (cid:27) . (2.17)We supply Eq. (2.14) with the following boundary conditions: P ε ( x , t ) = 0 , x ∈ Γ( t ) , (2.18) ∂P ε ∂n ( x , t ) = 0 , x ∈ Γ( t ) . (2.19) Note that in the axisymmetric case formula (2.14) coinsides with formula [4, (8)]. ZZ ω ε ( t ) P ε ( x , t ) d x = F ( t ) (2.20)connects the external load F ( t ), unknown contact pressure P ε ( x , t ), and unknown contactdomain ω ε ( t ). In order to check the content of formula (2.9) we consider here a special case, namely, wesuppose that the lower part cartilage layer is plane and rigid (the same assumption wasemployed in [14]), it means that µ s, = ∞ and R (1)1 = R (1)2 = ∞ , i.e.,Φ (1) ≡ , Φ ≡ Φ (2) . In this case we have got the following equation for determination of the contact domain ω ( t )in the form similar to (2.9):∆ P ( x , t ) + χ t Z ∆ P ( x , τ ) dτ = m (cid:16) Φ( x ) − δ ∗ ( t ) − ∇ e Φ( x ) · ∇ P ( x , t ) (cid:17) . (2.21)Here we will have m = 3 µ s, h , χ = 3 µ s, k h . (2.22)At the same time, small changes have to be made in the right-hand side of Eq. (2.21) asfollows: Φ( x ) = x R (2)1 + x R (2)2 , ˜Φ( x ) = h a x µ s, R (2)1 + h a x µ s, R (2)1 . Thus Eq. (2.21) can be rewritten as∆ P ( x , t ) + 3 µ s, k h t Z ∆ P ( x , τ ) dτ = 3 µ s, h x R (2)1 + x R (2)2 − δ ∗ ( t ) ! (2.23) − a h " x R (2)1 ∂ x P ( x , t ) + x R (2)2 ∂ x P ( x , t ) . (2.24)It can be easily checked that in the axisymmetric case Eq. (2.23) reduces to the governingdifferential equation obtained in [4]. 6 General relationships between the solution compo-nents
In our model we assume that the external load is non-decreasing. Thus, the contact domainis monotonically expanded, i.e. ω ε ( t ) ⊆ ω ε ( t ) , ∀ t ≤ t . (3.1)It is convenient to suppose also that the contact pressure is defined on the whole plane. Forthis we simply extend the density P ε ( x , t ) by assuming that P ε ( x , t ) = 0 , ∀ x ω ε ( t ) . (3.2)Integrating (2.14) over contact domain ω ( t ), we get ZZ ω ( t ) ∆ P ε ( x , t ) d x + χ ZZ ω ( t ) t Z ∆ P ε ( x , τ ) dτ d x == µ ZZ ω ( t ) (Ψ ( x ) − δ ε ( t )) d x − εµ ZZ ω ( t ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x . (3.3)For simplicity of notation, we omit here (and everywhere in the next two sections) thesubindex ε in ω ε .From the monotonicity of the contact domain (3.1) and assumption (3.2), it follows thatthe second integral on the left-hand side can be written in the form ZZ ω ( t ) t Z ∆ P ε ( x , τ ) dτ d x = t Z ZZ ω ( t ) ∆ P ε ( x , τ ) d x dτ. (3.4)Using the second Green’s formula ZZ ω ( t ) ( u ( x )∆ v ( x ) − v ( x )∆ u ( x )) d x = Z Γ( t ) (cid:18) u ( x ) ∂v∂n ( x ) − v ( x ) ∂u∂n ( x ) (cid:19) ds (3.5)with u ≡ v = P ε ( x , t ) we get the following relation in view of the boundary condition(2.19): ZZ ω ( t ) ∆ P ε ( x , τ ) d x = Z Γ( t ) ∂P ε ∂n ( x , s ) ds = 0 , ∀ τ ≤ t. (3.6)Therefore, the both integrals on the left-hand side of (3.3) vanish.7urther, we use the first Green’s formula ZZ ω ( t ) ( ϕ ∆ ψ + ∇ ϕ · ∇ ψ ) d x = Z Γ( t ) ϕ ∂ψ∂n ds (3.7)with ψ ( x ) = Ψ ( x ) and ϕ ( x ) = P ε ( x , t ). In this case the integral on the right-hand sidevanishes in view of (2.18), and we obtain the relation ZZ ω ( t ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x = − ZZ ω ( t ) P ε ( x , t )∆Ψ ( x ) d x = − e ) F ( t ) , (3.8)where we used the equilibrium equation (2.20) and the identity∆Ψ ( x ) = 2(1 + e ) (3.9)with e being defined in (2.15).In what follows, it is convenient to have the following notation for the integrals of theproduct of k -th power of the function Ψ and l -th power of the function Ψ : A k,l ( ω ) = ZZ ω Ψ k ( x )Ψ l ( x ) d x > , k, l = 0 , , , . . . (3.10)In particular, A , ( ω ) 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 avoidcumbersome expressions. Computations of A k,l ( ω ) for the elliptic domain (2.17) we includedinto Appendix (see Section 6.1).Taking into account Eqs. (3.6) and (3.8), we get δ ε ( t ) = A , ( ω ε ( t )) A , ( ω ε ( t )) + 2(1 + e ) εA , ( ω ε ( t )) F ( t ) . (3.11)This formula allows us to compute the contact approach δ ε ( t ) as a function of the totalexternal force F ( t ) and the main axes of the ellipse describing the shape of the contact zone,which in fact depends on time too. In order to write out a more informative equation for the contact load, we use the followingtrick. We multiply both sides of (2.14) by the function v ( x ) = Ψ ( x ) and integrate the8btained equation over the contact domain ω ( t ) ZZ ω ( t ) Ψ ( x )∆ P ε ( x , t ) d x + χ ZZ ω ( t ) t Z Ψ ( x )∆ P ε ( x , τ ) dτ d x == µ ZZ ω ( t ) Ψ ( x )Ψ ( x ) d x − µδ ε ( t ) ZZ ω ( t ) Ψ ( x ) d x − µε ZZ ω ( t ) Ψ ( x ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x . (3.12)Let us calculate the integrals in this relation by using Green’s formulas. For the first integralon the left-hand side we use formula (3.5) with u = Ψ , v = P ε and the boundary conditions(2.18), (2.19). Hence, we obtain ZZ ω ( t ) Ψ ( x )∆ P ε ( x , t ) d x = ZZ ω ( t ) ∆Ψ ( x ) P ε ( x , t ) d x . Now taking into account (3.9), we get ZZ ω ( t ) Ψ ( x )∆ P ε ( x , t ) d x = 2(1 + e ) F ( t ) . (3.13)For the second integral on the left-hand side, we apply the same approach, but interchangefirst the integrals over ω ε ( t ) and over τ ∈ (0 , t ) exploiting the load monotonicity. Therefore,we arrive at the equation ZZ ω ( t ) t Z Ψ ( x )∆ P ε ( x , τ ) dτ d x = t Z ZZ ω ( t ) Ψ ( x )∆ P ε ( x , τ ) dτ d x = 2(1 + e ) t Z F ( τ ) dτ. (3.14)For the first and second integrals on the right-hand side, we simply use the notation (3.10),which gives ZZ ω ( t ) Ψ ( x )Ψ ( x ) d x = A , ( b ; β ) , ZZ ω ( t ) Ψ ( x ) d x = A , ( b ; β ) . (3.15)Finally, for the third integral on the right-hand side, we make use of the following simpleformula which follows immediately from the definition of Ψ :Ψ ∇ Ψ = 12 ∇ Ψ . ZZ ω ( t ) Ψ ( x ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x = − ZZ ω ( t ) ∆Ψ ( x ) P ε ( x , t ) d x . By applying the second Green’s formula (3.5) with u = P ε , v = Ψ , and the boundaryconditions (2.18), (2.19), we represent this integral in the form ZZ ω ( t ) Ψ ( x ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x = − ZZ ω ( t ) Ψ ( x )∆ P ε ( x , t ) d x . (3.16)This integral still contains the unknown density of contact pressure P ε ( x , t ). Let us define M ( j ) P ε ( t ) ≡ ZZ ω ( t ) Ψ j ( x )∆ P ε ( x , t ) d x . (3.17)Now we rewrite the relation (3.12) by using the results for all integrals (3.13)–(3.16) inthe following form:2(1 + e ) K F ( t ) = µA , ( ω ε ( t )) − µδ ε ( t ) A , ( ω ε ( t )) + µε M (2) P ε ( t ) . (3.18)Here, we have introduced the Volterra operator K as follows: K F ( t ) = F ( t ) + χ t Z F ( τ ) dτ. (3.19)Note that the integral in the right-hand side of the equation (3.18) allows to continuethe same procedure to deliver an asymptotic estimate for this equation.We continue to proceed with Eq. (3.18) on the next steps. M ( j ) P ε ( t ) Now we proceed to calculate the last integral in (3.18). For this we multiply the governingintegral equation (2.14) by Ψ j ( x ) ( j ≥
2) and integrate over contact domain ω ( t ): ZZ ω ( t ) Ψ j ( x )∆ P ε ( x , t ) d x + χ ZZ ω ( t ) t Z Ψ j ( x )∆ P ε ( x , τ ) dτ d x == µ ZZ ω ( t ) Ψ j ( x )Ψ ( x ) d x − µδ ε ( t ) ZZ ω ( t ) Ψ j ( x ) d x − µε ZZ ω ( t ) Ψ j ( x ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x . (3.20)10y using the same argument as on the previous step, we get KM ( j ) P ε ( t ) = µA ,j − µδ ε ( t ) A ,j ( a ; β ) − µε ZZ ω ( t ) Ψ j ( x ) ∇ Ψ ( x ) · ∇ P ε ( x , t ) d x . (3.21)For the last integral we use the relationsΨ j ( x ) ∇ Ψ ( x ) = 1 j + 1 ∇ Ψ j +12 ( x )and ZZ ω ( t ) ∇ Ψ j +12 ( x ) · ∇ P ε ( x , t ) d x = − ZZ ω ( t ) ∆Ψ j +12 ( x ) P ε ( x , t ) d x . Therefore, the integral M ( j ) P ε ( t ) = µ K − (cid:26) A ,j ( ω ε ( t )) − δ ε ( t ) A ,j ( ω ε ( t )) + εj + 1 KM ( j +1) P ε ( t ) (cid:27) (3.22)has been obtained as a solution of the integral equation (3.21). Here the inverse operator K − is defined by the formula K − Y ( t ) = Y ( t ) − χ t Z Y ( τ ) e − χ ( t − τ ) dτ. (3.23)Performing the same computation, we obtain the following representation for the integralin the right-hand side of (3.18): M (2) P ε ( t ) = N X j =1 ε j − ( j + 1)! µ j K − j { A ,j +1 ( ω ε ( t )) − δ ε ( t ) A ,j +1 ( ω ε ( t )) } + 2 ε N ( N + 2)! µ N K − N M ( N +2) P ε ( t ) . (3.24)Substituting this representation into Eq. (3.18), we finally get2(1 + e ) K F ( t ) = N X j =0 ε j ( j + 1)! µ j +1 K − j { A ,j +1 ( ω ε ( t )) − δ ε ( t ) A ,j +1 ( ω ε ( t )) } + ε N +1 ( N + 2)! µ N +1 K − N M ( N +2) P ε ( t ) , (3.25)or equivalently2(1 + e ) K N +1 F ( t ) = N X j =0 ε j ( j + 1)! µ j +1 K N − j { A ,j +1 ( ω ε ( t )) − δ ε ( t ) A ,j +1 ( ω ε ( t )) } ε N +1 ( N + 2)! µ N +1 M ( N +2) P ε ( t ) . (3.26)The latter relation allows us to determine the problem parameters asymptotically with anyprescribed accuracy.Note that apart from the fact that the shapes of the contacting bones are ellipticalparaboloids, 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 approximatelyrepresented by an ellipse. This will be done later. Remark 1.
For every t for which the contact pressure P ε ( t ) is bounded and the contactregion ω ( t ) belongs to a bounded domain, the remainder ε N +1 ( N +2)! µ N +1 M ( N +2) P ε ( t ) in formula(3.26) tends to zero as N → ∞ . Thus, the series corresponding to the sum on the righthand-side of (3.26) is converging. First, we get solution of the problem for ε = 0. In this case Eq. (2.14) has the form∆ P (0) ( x , t ) + χ t Z ∆ P (0) ( x , τ ) dτ = µ (cid:0) Ψ ( x ) − δ (0) ( t ) (cid:1) , (4.1)where Ψ ( x ) is defined in (2.15). Since we know from [5] that the contact zone is an ellipseat this stage of approximation we will have δ ε = δ (0) ( t ) = δ ε ( b ( t ); β ( t )) = A , ( ω ( t )) A , ( ω ( t )) . (4.2)Using formula (4.2) and calculations presented in Section 6.1 (see formula (6.6)), onecan find that A , ( ω ( t )) = πb β , A , ( ω ( t )) = πb β (cid:0) β + e (cid:1) , (4.3)and therefore δ (0) ( t ) = b ( β + e )4 β . (4.4)Note that formulas (4.3) and (4.4) contain two known constants e and e defined in (2.15)and two still unknown functions b ( t ) and β ( t ), which are the main semi-axis and theeccentricity of the ellipse ω ( t ) = (cid:26) x ∈ R : x b ( t ) + β ( t ) x b ( t ) ≤ (cid:27) . (4.5)12he leading terms in (3.26) imply (for N = 0) the following equation:2(1 + e ) K F ( t ) = µA , ( ω ( t )) − µδ (0) ( t ) A , ( ω ( t )) . (4.6)Here, K is the Volterra integral operator defined in (3.19).Analogously, using some results from Section 6.1 (see, in particular, formula (6.6)), weobtain A , ( ω ( t )) = πb β (cid:0) β + e (cid:1) (4.7)and A , ( ω ( t )) = πb β (cid:8) β + ( e + e ) β + 3 e e (cid:9) , (4.8)and thus 2(1 + e ) K F ( t ) = µ πb β (cid:8) β − ( e + e ) β + 3 e e (cid:9) . (4.9)To find the functions b ( t ) and β ( t ) together with the pressure distribution over thecontact zone, P (0) ( x , t ), we follow [5] and introduce a new unknown function p (0) ( x , t ) = P (0) ( x , t ) + χ t Z P (0) ( x , τ ) dτ = K P (0) ( x , t ) . (4.10)In the case of monotone external load, this function should satisfy the Poisson equation(following from (2.9)) ∆ p (0) ( x , t ) = µ (cid:0) Ψ ( x ) − δ (0) ( t ) (cid:1) , x ∈ ω ( t ) , (4.11)with the boundary conditions (2.18), (2.19).It is customary to rewrite this relation in the form G ( x , t ) = 0 , (4.12)where G ( x , t ) = G ( b , β , δ ) (4.13) ≡ ∆ p (0) ( x , t ) − µ (cid:0) Ψ ( x ) − δ (0) ( t ) (cid:1) , x ∈ ω ( t ) . (4.14)Bearing in mind that the function Ψ ( x ) is a quadratic polynomial (compare with (2.15)),it is natural to look for the solution of such problem in the form of a polynomial in x , x ofthe fourth degree, that is p (0) ( b , β , η , x , t ) = η ( t ) (cid:18) − x b − β x b (cid:19) Q ( x , x ) . (4.15)13ote that the term in the brackets vanishes on the boundary ω , and thus the condition(2.18) is satisfied automatically.In Section 6.2, it has been shown that Q is a polynomial of the second order having theform Q ( x , x ) = (cid:18) − x b − β x b (cid:19) , (4.16)so that p (0) ( x , x ; t ) = η ( t ) (cid:18) − x b − β x b (cid:19) . (4.17)Taken into account this representation we arrive at the following relations (see Section 6.3): η ( t ) = µδ (0) ( t )4(1 + β ) b , (4.18) η ( t ) = µb β ) = µb β ) , (4.19) η ( t ) = µb e β + 6 β ) = µb e β (1 + 3 β ) . (4.20)This system allows us to determine the unknown functions b ( t ) and β ( t ). Indeed,eliminating η from the last two equations, we get a bi-quadratic equation defining the valueof the parameter β , i.e., 3 β + (1 − e ) β − e = 0 . (4.21)By definition, β is a positive parameter, thus the unique positive solution of (4.21) has theform β = s ( e −
1) + p e + 34 e + 16 . (4.22)Note that at the zero-approximation the parameter β does not depend on time. The otherparameter, η ( t ), can be computed directly from (4.19) or (4.20), if one knows the remainingconstant b ( t ). Moreover, taking into account (4.18) and (4.4), one can use an equivalentformula η ( t ) = µb ( β + e )16 β (1 + β ) . (4.23)In the same way, one can offer, in addition to (4.4), two equivalent representations forthe indentation parameter δ (0) ( t ) = 1 + β β b ( t ) = (1 + β ) e β (1 + 3 β ) b ( t ) . (4.24)Finally, the major semi-axis b of the ellipse ω is determined as follows: b ( t ) = F ( t ) + χ t Z F ( τ ) dτ (cid:18) β (1 + e ) µπ (3 β − β ( e + e ) + 3 e e ) (cid:19) / . (4.25)14ote that the parameters b , η as well as the indentation, δ , depend on time t in contrastto the ellipse eccentricity β .Now, it remains only to find the pressure over the contact area. Using (4.10) and (4.17),we get P (0) ( b , β , η , x , x , t ) = K − (cid:0) η ( t ) Q ( x , x ) (cid:1) . (4.26)If ( x , x ) belongs to the initial contact zone, i.e. 1 − x b ( t ) − β x b ( t ) >
0, then P (0) ( x , x , t ) = η ( t ) (cid:18) − x b ( t ) − β x b ( t ) (cid:19) − χ t Z η ( τ ) (cid:18) − x b ( τ ) − β x b ( τ ) (cid:19) e − χ ( t − τ ) dτ. (4.27)If ( x , x ) lies outside of the initial contact zone, i.e. 1 − x b ( t ) − β x b ( t ) <
0, then P (0) ( x , x , t ) = η ( t ) (cid:18) − x b ( t ) − β x b ( t ) (cid:19) − χ t Z t ∗ ( x ,x ) η ( τ ) (cid:18) − x b ( τ ) − β x b ( τ ) (cid:19) e − χ ( t − τ ) dτ. (4.28)The critical moment of time t ∗ is determined by the formula b ( t ∗ ) = x + β x . Using (4.25), we get F ( t ∗ ) + χ t ∗ Z F ( τ ) dτ = µπ β (cid:18) β − β ( e + e ) + 3 e e e (cid:19) ( x + β x ) . (4.29)If the load is stepwise, we have F ( t ) = F . Hence, we find that t ∗ = µπ β χF (cid:20) (3 β − β ( e + e ) + 3 e e )1 + e ( x + β x ) (cid:21) − χ . (4.30)Note that in this case b ( t ∗ ) = 96 β (1 + e )(1 + χt ∗ ) µπ (3 β − β ( e + e ) + 3 e e ) F . (4.31)This finishes the zero iteration step. Note that the results of this Section after changingthe notation coincide with those obtained in [4].15 .2 First-order approximation problem For the next steps we consider an appropriately deformed contact domain ω (1) ε , defined as aperturbation of the zero-order one ω . Namely, we assume that it can be written in the form ω (1) ε = ω (1) ε ( t ) = n ( x , x ) : Q ( x , t ) + εQ ( x , t ) ≥ o , (4.32)where unknown polynomials are taken in the forms Q ( x , t ) = Q ( x , β , b ) , (4.33) Q ( x , t ) = a ( t ) x + a ( t ) x x + a ( t ) x . (4.34)Note that for ε = 0 the solution form coincides with (4.5), if one take b ≡ b , β ≡ β .The idea behind such choice of the asymptotic anzatz is to satisfy the boundary condi-tions (2.18) and (2.19) automatically. This will be archived by putting P (1) ε = K − (cid:16) η (1) ( t ) (cid:0) Q ( x , x , β ( t ) , b ( t )) + εQ ( x , t ) (cid:1) (cid:17) . (4.35)Now, when the boundary conditions are valid, we will satisfy the governing equation(2.9). Note that P (1) ε = P + εP + O ( ε ) , (4.36)where p j = K ( P j ), j = 0 ,
1, and p = η (1) ( t ) (cid:18) − x b ( t ) − β ( t ) x b ( t ) (cid:19) , (4.37) p = 2 η (1) ( t ) (cid:18) − x b ( t ) − β ( t ) x b ( t ) (cid:19) Q ( x , t ) . (4.38)Substituting this representation into Eq. (2.9), we obtain K (cid:0) ∆( P (0) + εP + O ( ε )) (cid:1) = µ (cid:0) Ψ − δ (1) ε − ε ∇ Ψ · ( ∇ P (0) + ε ∇ P (1) + O ( ε )) (cid:1) , (4.39)where the parameter δ (1) ε is represented in the same form as P (1) ε , i.e., δ (1) ε = δ + εδ + O ( ε ) = δ (1) + O ( ε ) . (4.40)We can write Eq. (4.39) with the accuracy to the terms of O ( ε ) as follows:∆ p (0) + ε ∆ p = µ (cid:0) Ψ − δ (1) − ε ∇ Ψ · ∇ P (0) (cid:1) . (4.41)An extended variant of this equation can be written by using the definition of all com-ponents of the equation and by comparing coefficients at different powers of x , x , so that16 η (1) b (1 + β ) = − µδ (1) , (4.42)4 η (1) (cid:20) β b + ε (6 a + a ) (cid:21) = µ (1 − εθ , ) , (4.43)4 η (1) (cid:20) β (1 + 3 β ) b + ε ( a + 6 a ) (cid:21) = µ ( e − εe θ , ) , (4.44) − ε η (1) b ( a β + a (1 + β ) + a ) = 8 εµ (1 + e ) θ , , (4.45) − ε η (1) b ( a (15 + β ) + a ) = 8 εµθ , , (4.46) − ε η (1) b ( a (15 β + 1) + a β ) = 8 εµe θ , , (4.47)where θ k, l ( t ) = K − (cid:0) η (1) b − k β l (cid:1) , k, l = 0 , , . (4.48)In the system (4.42)–(4.47) we have 6 equations and 7 unknowns: η (1) ( t ) , δ (1) ε , b ( t ) , β ( t ),and a , a , a (coefficients of the polynomial Q ). Therefore, we have to add an extraequation to the above system, namely δ (1) ( t ) = A , ( ω ε ( t )) A , ( ω ε ( t )) + 2(1 + e ) εA , ( ω ε ( t )) F ( t ) , (4.49)where F ( t ) can be represented in the form F ( t ) = Z Z ω (1) ε P (1) ε ( x , t ) d x . We also make use of Eq. (3.26) written for this approximation step with the accuracy of O ( ε ) in the form2(1 + e ) K F ( t ) = X j =0 ε j ( j + 1)! µ j +1 K − j (cid:8) A ,j +1 ( ω ε ( t )) − δ (1) ( t ) A ,j +1 ( ω ε ( t )) (cid:9) . (4.50) Remark 2.
Note that putting ε = 0, the system (4.42)–(4.47), (4.49) transforms to theprevious case evaluated in the previous section. Remark 3.
In the case when ε >
0, the system (4.42)–(4.47), (4.49) has to be solvednumerically. Note that the parameter ε in the last three equations (4.45) – (4.47) can becanceled. We left these multipliers here to explain the limiting case ( ε = 0).17 Discussion and conclusion
First of all, observe that at t = 0, the contact problem for biphasic layers reduces to that forelastic incompressible layers. The contact problem in the latter case were studied in a numberof papers [1, 9, 10, 16], however, without taking into account the tangential displacements.To solve the resulting problem (4.42)–(4.47) and (4.49), we suggest the following iterativealgorithm: • Taking ε = 0, we have computed all values η, b, β, δ = η , b , β , δ from the zero-orderapproximation. • Having them we can compute the quantity θ k, l ( t ) from (4.48), • Then, from the system of three equations (4.45)–(4.47) we compute the constants a , a , a assuming the values of η, b, β as above. • Finally from the system of four equations (4.42)–(4.44) and (4.49) considering the right-hand side known (computed by the values know from the previous computations), wefound new values η, b, β, δ and compare them with the previous computations. If therequired accuracy has achieved we stop the computation, if not we are going to thesecond step of this iterative procedure.We note that formulas (2.4) and (2.5) for the vertical and tangential displacementscontain different powers of parameters ǫ , namely, ǫ and ǫ , respectively. Note also that ouranalysis (with the values of another parameters taken into account) shows, that the roleof these magnitudes (vertical and tangential displacements) is quite opposite. In the finalequation (see (2.14)) the leading terms, corresponding to the vertical displacement, containthe zero power of the new small parameter ε , but the leading terms, corresponding to thetangential displacements, contain the first power of ε . References [1] Aleksandrov, V.M. Asymptotic solution of the axisymmetric contact problem for anelastic layer of incompressible material.
J. Appl. Math. Mech. , 589–593 (2003)[2] Argatov, I. Development of an asymptotic modeling methodology for tibio-femoral con-tact in multibody dynamic simulations of the human knee joint. Multibody Syst. Dyn. , 3–20 (2012).[3] Argatov, I. Contact problem for a thin elastic layer with variable thickness: Applicationto sensitivity analysis of articular contact mechanics. Appl. Math. Model. , 8383–8393(2013). 184] Argatov, I., Mishuris, G. Axisymmetric contact problem for a biphasic cartilage layerwith allowance for tangential displacements on the contact surface. Eur. J. Mech.A/Solids , 1051–1064 (2010).[5] Argatov, I., Mishuris, G. Elliptical contact of thin biphasic cartilage layer: Exact solu-tion for monotonic loading. J. Biomech. , 759–761 (2011).[6] Argatov, I., Mishuris, G. Frictionless elliptical contact of thin viscoelastic layers bondedto rigid substrates. Applied Math. Model. , 3201–3212 (2011).[7] Argatov, I., Mishuris, G. Contact problem for thin biphasic cartilage layers: perturba-tion solution. Quart. J. Mech. Appl. Math. , 297–318 (2011).[8] Ateshian, G.A., Lai, W.M., Zhu, W.B., Mow, V.C. Anasymptotic solution for the con-tact of two biphasic cartilage layers. J. Biomech. , 1347–1360 (1994).[9] Barber, J.R.: Contact problems for the thin elastic layer. Int. J. Mech. Sci. , 129–132(1990).[10] Chadwick, R.S. Axisymmetric indentation of a thin incompressible elastic layer. SIAMJ. Appl. Math. , 1520–1530 (2002).[11] Hunziker, E. B. Articular cartilage repair: basic science and clinical progress. A reviewof the current status and prospects. Osteoarthritis and Cartilage . , 432–463 (2001).[12] Johnson, K. L. Contact Mechanics . Cambridge, Cambridge University Press (1985).[13] Owen, J. R., Wayne, J. S. Contact models of repaired articular surfaces: influence ofloading conditions and the superficial tangential zone.
Biomech Model Mechanobiol. ,461–471 (2011).[14] Wu, J. Z., Herzog, W., and Ronsky, J. Modeling axi-symmetrical joint contact withbiphasic cartilage layers–An asymptotic solution. J. Biomech. , 1263–1281 (1996).[15] Wu, J. Z., Herzog, W., and Epstein, M. An improved solution for the contact of twobiphasic cartilage layers. J. Biomech. , 371–375 (1997).[16] Yang, F. Indentation of an incompressible elastic film. Mech. Mater. , 275–286 (1998).19 Appendix A kl Here, we compute the values of the constants A k,l ( b ; β ) = ZZ ω ( t ) Ψ k ( x )Ψ l ( x ) d x > , k, l = 0 , , , . . . First of all, we note that an unknown contact domain ω ( t ) is of the same type as sections ofthe initial gap elliptical paraboloid, i.e., it is an ellipse coaxial to the ellipse ω ( t ) = ω ε ( t ) = (cid:26) x ∈ R : x b ( t ; ε ) + x β ( t ; ε ) b ( t ; ε ) ≤ (cid:27) . In order to avoid long formulas, we use the short notation for ω ( t ), writing all parameterswithout variables they depend on, i.e., ω ( t ) = (cid:26) x ∈ R : x b + x β b ≤ (cid:27) . Performing the standard change of variables x = br cos θ, x = bβ sin θ, we represent the integral for A k,l ( t ) in the form A k,l ( b ; β ) = Z π Z (cid:18) b r cos θ + b e β r sin θ (cid:19) k (cid:18) b r cos θ + b e β r sin θ (cid:19) l b β rdrdθ = b k +2 l +2 β Z r k +2 l +1 dr π Z k X i =0 k ! i !( k − i )! e i β i sin i θ cos k − i θ (6.1) × l X j =0 l ! j !( l − j )! e j β j sin j θ cos l − j θdθ = b k +2 l +2 (2 k + 2 l + 2) β π Z k X i =0 k ! i !( k − i )! e i β i sin i θ cos k − i θ × l X j =0 l ! j !( l − j )! e j β j sin j θ cos l − j θdθ. (6.2)20ince the trigonometric functions are presented here only in even powers, then the lastintegration can be performed over the interval [0 , π/
2] as follows: A k,l ( b ; β ) = 4 b k +2 l +2 (2 k + 2 l + 2) β π/ Z k X i =0 k ! i !( k − i )! e i β i sin i θ cos k − i θ × l X j =0 l ! j !( l − j )! e j β j sin j θ cos l − j θdθ = 4 b k +2 l +2 (2 k + 2 l + 2) β k X i =0 k ! i !( k − i )! e i β i × l X j =0 l ! j !( l − j )! e j β j π/ Z sin i +2 j θ cos k − i +2 l − j θdθ. (6.3)The integrals in (6.3) are calculated by using formulas π/ Z sin p θ cos q θdθ = 12 Γ( p + 1 / q + 1 / p + q + 1) , p, q > , (6.4)and Legendre’s duplication formula for the Gamma-functionΓ( n + 1 /
2) = √ π Γ(2 n )2 n − / Γ( n ) , n ∈ N , (6.5)as well as the relation Γ( n + 1) = n !. Finally, we arrive at the following representation of A k,l = A k,l ( b ; β ) valid for all k, l ∈ N = N ∪ { } : A k,l = 2 πb ( b/ k +2 l β (2 k + 2 l + 2)( k + l )! k X i =0 k ! i !( k − i )! e i β i (6.6) × l X j =0 l ! j !( l − j )! e j β j (2 i + 2 j )!(2 k − i + 2 l − j )!( i + j )!( k − i + l − j )! . Q In order to determine the coefficients of the polynomial Q ( x , x ) = 1 + q , x + q , x + q , x + q , x x + q , x , we need to compute the normal derivative of the unknown functions p (0) (4.15) along theelliptic boundary Γ: ∂p (0) ∂n | Γ = ∇ p (0) · −→ n | Γ = η ( t ) (cid:18) − x b − β x b (cid:19) Q | Γ = 0 . (6.7)21ere we take into account the fact that, since the contact domain is an ellipse (4.5), thetangential and normal vectors to the boundary Γ = ∂ Ω are given by −→ r = (cid:0) − β x , x (cid:1) , −→ n = (cid:0) x , β x (cid:1) . (6.8)Then, to satisfy the boundary condition (2.19) the following equation should be valid: Q | Γ = 0 . (6.9)This, in turn, is equivalent to the representation Q ( x , x ) = (cid:18) − x b − β x b (cid:19) . (6.10) Since p (0) ( x , t ) = p (0) ( x , x , t ) = η ( t ) (cid:18) − x b − β x b (cid:19) , (6.11)we have ∂p (0) ∂x = 2 η ( t ) (cid:18) − x b − β x b (cid:19) · (cid:18) − x b (cid:19) , (6.12) ∂ p (0) ∂x = 2 η (cid:20) − b (cid:18) − x b − β x b (cid:19) + 2 x b x b (cid:21) . Therefore, by straightforward computations, we find that ∂ p (0) ∂x = 2 η (cid:20) − b + 6 x b + 2 β x b (cid:21) . (6.13) ∂p (0) ∂x = 2 η (cid:18) − x b − β x b (cid:19) · (cid:18) − β x b (cid:19) , (6.14) ∂ p (0) ∂x = 2 η (cid:20) − β b (cid:18) − x b − β x b (cid:19) + 2 β x b β x b (cid:21) . Thus, we obtain ∂ p (0) ∂x = 2 η ( t ) (cid:20) − β b + 2 β x b + 6 β x b (cid:21) . (6.15)Substituting (6.13) and (6.15) into the main equation G ( b , β , δ ) ≡ ∆ p (0) ( x , t ) − µ (cid:0) Ψ ( x ) − δ (0) ( t ) (cid:1) = 0 , (6.16)where G = 2 η ( t ) (cid:20) ( −
2) 1 + β b + (cid:18) β b (cid:19) x + (cid:18) β + 2 β b (cid:19) x (cid:21) − µ (cid:0) Ψ ( x , x ) − δ (0) ( t ) (cid:1) , ( x ) = Ψ ( x , x ) = x + e x , one concludes that the expression for G is represented by a second order polynomial withrespect to the independent variables x and x in the following form: G ( b , β , η , δ (0) ) = q ( b , β , η , δ (0) ) + q ( b , β , η ) x + q ( b , β , η ) x . (6.17)Here the coefficients are defined as follows: q ( b , β , η , δ (0) ) = 4 η µb (1 + β ) − δ (0) , (6.18) q ( b , β , η ) = 4 η b (3 + β ) − µ, (6.19) q ( b , β , η ) = 4 η β b (1 + 3 β ) − µe . (6.20) Taking into account (4.37), we can represent p ( x , t ) in the form p ( x , t ) = η (1) ( t ) (cid:18) − x b − β x b + 2 β x x b + x b + β x b (cid:19) . (6.21)Hence, applying the Laplace equation, we get∆ p ( x , t ) = η (1) ( t ) (cid:18) − b (1 + β ) + x b (3 + β ) + x β b (1 + 3 β ) (cid:19) . (6.22)Next, by using representation (4.38), we can write p ( x , t ) in the form p ( x , t ) = 2 η (1) ( t ) (cid:18) a x + a x x + a x − a x b − a x x b − a x x b (6.23) − a β x x b − a β x x b − a β x b (cid:19) . Therefore, we obtain∆ p ( x , t ) = 2 η (1) ( t ) (cid:18) (12 a + 2 a ) x + (2 a + 12 a ) x (6.24) − β a + 12 a (1 + β ) + 12 a b x x − a (30 + 2 β ) + 2 a b x − a β + a (2 + 30 β ) b x (cid:19) . (6.25)23e also use the following representations:Ψ j ( x ) = x + e j x , j = 1 , . Thus, applying the gradient operator, we simply get ∇ Ψ ( x ) = (cid:0) x , e x (cid:1) and ∇ P ( x , t ) = (cid:0) K − ∇ p ( x , · ) (cid:1) ( t ) . It yields the following representation: ∇ Ψ ( x ) · ∇ P ( x , t ) = − (cid:18) K − (cid:20) η (1) (cid:18) − x b − β x b (cid:19) (cid:18) x b + e β x b (cid:19)(cid:21)(cid:19) ( t ) (6.26)= − x (cid:18) K − (cid:18) η (1) b (cid:19)(cid:19) ( t ) − e x (cid:18) K − (cid:18) η (1) β b (cid:19)(cid:19) ( t ) + 8 x (cid:18) K − (cid:18) η (1) b (cid:19)(cid:19) ( t )+8(1 + e ) x x (cid:18) K − (cid:18) η (1) β b (cid:19)(cid:19) ( t ) + 8 e x (cid:18) K − (cid:18) η (1) β b (cid:19)(cid:19) ( t )=: − x θ , ( t ) − e x θ , ( t ) + 8 x θ , ( t ) + 8(1 + e ) x x θ , ( t ) + 8 e x θ , ( t ) . Here we have introduced the notation θ k, l = (cid:0) K − (cid:0) η (1) b − k β l (cid:1)(cid:1) ( t ) ..