Numerical analysis of the loosened total hip replacements (THR)1

Mathematics and Computers in Simulation 50 (1999) 285±304

Numerical analysis of the loosened total hip replacements (THR)1 J. Nedomaa,*, M. BartosÏb, H. HornaÂtovaÂc, Z. KestrÏaÂnek, Jr.a,c, J. StehlõÂkd a

Institute of Computer Science AS CR, Pod VodaÂrenskou veÏzÏÂõ 2, 182 07 Prague 8, Czech Republic b Central Military Hospital, U Vojenske nemocnice 1200, 160 00 Prague 6, Czech Republic c JaÂchymov Spa, LeÂcÏebne laÂzneÏ a.s., 362 51 JaÂchymov, Czech Republic d Orthopaedic Clinic, 3rd Medical Faculty of the Charles University, SÏrobaÂrova 50, 100 00 Prague 10, Vinohrady, Czech Republic

Abstract The paper deals with mathematical simulation of a loosened hip joint, simulation of mechanical processes taking place during static and dynamic burdening and their mathematical description. The main result of the paper consists of the existence theorem and the analysis of the loosened total hip replacement (THR). # 1999 IMACS/Elsevier Science B.V. All rights reserved. Keywords: Quasi-static contact problems in elasticity; Biomechanics; Loosened THR

1. Introduction In the past 30 years, a lot of experience has been assembled with a change in the surgical technique of osteoarthritis of the human joints, namely the hip joint. In former times, readjustment of the bone was the only possible therapy, but as a consequence of the extensive development in the area of medical techniques, a total joint replacement has become more and more widespread. Although many positive results have been obtained, the length of life of the total hip replacement (THR) is limited. Its service life (vitality) depends on several factors. Firstly on the quality of the applied surgical technique and on the quality of the THR used and secondly on the high values of static and dynamic loading. As orthopaedic observations show, the THR starts to loosen after a certain period. Disengagement starts to loosen in the area of the upper part of the stem of the artificial endoprothesis (THR) and continues to its lower parts in dependence on the static and dynamic loading of the artificial hip joint. This paper will discuss the problem of simulation of the mechanical processes taking place during the static and dynamic loading and the mathematical analysis of the model problem [3,4,6±8,10]. ÐÐÐÐ * Corresponding author. E-mail: [email protected] 1 Supported by the GA AS CR grant No. 308-95-0304. 0378-4754/99/$20.00 # 1999 IMACS/Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 9 9 ) 0 0 0 8 2 - 8

286

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

2. The model Let the investigated part of the human skeleton being in the initial stress±strain state S occupy a Ndimensional (N ˆ 2; 3) domain with a Lipschitz boundary @ . Let ˆ 4ˆ1  , where

 ;  ˆ 1; . . . ; 4 denote stage by stage the pelvis, the artificial acetabulum, the artificial head of the femur (i.e. the THR) and the femur. Let the boundary @ consists of three parts ÿ , ÿu, ÿc such as @ ˆ ÿ [ ÿu [ ÿc , where ÿ denotes the part of the boundary, where the loading is prescribed, ÿu denotes the part of the boundary where the skeleton is fixed (i.e. one part of the femur and of the pelvis) and ÿc denotes the contact boundary between the artificial acetabulum and the artificial head of the femur, and/or between the stem of the THR and the bone and/or between the artificial acetabulum and the pelvis. We shall assume that the problem studied is described, with sufficient accuracy, by a dynamic or by a quasi-static contact problem with friction. Let t 2 I  h0; Ti; T > 0:, Let n ˆ …ni †; t ˆ …ti † be the outward unit normal and unit tangential vectors, respectively, to the boundary @ . Let: ij …x; t† ˆ cijkl …x†ekl …u…x; t††;

i; j; k; l ˆ 1; . . . ; N;

be the stress tensor and: eij …u…x; t†† ˆ 1=2…@ui …x; t†=@xj ‡ @uj …x; t†=@xi †;

i; j ˆ 1; . . . ; N;

the strain tensor and i ˆ ij nj the stress vector. Let its normal and tangential components be defined by t ˆ ~  ÿ n n and the normal and tangential components of the displacement vector n ˆ i ni ˆ ij nj ni ; ~ by un ˆ ui ni ; ut ˆ u ÿ un n. Let us denote by uk,ul (the indices k,l denote the neighbouring components of the artificial hip joint studied) the displacements in the neighbouring components of the hip joint studied. All these quantities are functions of spatial coordinates x and time t. Then we have to solve the following problem: Problem (P). Find u(x,t), t  0 satisfying:   [ @ 2 ui …x; t† @ @uk …x; t†  ‡  …x†fi …x; t† ˆ 0; i; j; k; l ˆ 1; . . . ; N; on ÿ c …x†

  I; …x† ijkl 2 @t @xj @xl  (1) u…x; t† ˆ 0

on ÿu  I;

ij …x; t†nj ˆ Pi …x; t† on ÿ  I;

(2) (3)

ukn …x; t†ÿuln …x; t†  0; nkl …x; t†  nk …x; t† ˆ ÿnl …x; t†  0; …ukn …x; t†ÿuln …x; t††nkl …x; t† ˆ 0; on ÿkl c I; and kl if ukn …x; t† ÿ uln …x; t† ˆ 0 then j~  kl t …x; t†j  gc …x†; kl k l if j~  kl t …x; t†j < gc …x† then ut …x; t† ÿ ut …x; t† ˆ 0;

kl k l  kl if j~  kl t …x; t†j ˆ gc …x† then 9  0 such that ut …x; t† ÿ ut …x; t† ˆ ÿ…x†~ t …x; t†;

(4)

where gkl c …x† denotes the given frictional forces, (x) the density, cijkl(x) the elastic coefficients, f(x,t), P(x,t) the body and surface forces and …x†@ 2 ui …x; t†=@t2 represents the inertia forces. If the inertia

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

287

forces are neglected we have the quasi-static contact problem. Moreover, the elastic coefficients cijkl(x) satisfy: cijkl …x† ˆ cjikl …x† ˆ cklij …x† ˆ cijlk …x† on  ; 0 < c0  cijkl …x†ij kl jjÿ2  c1 < ‡1 8x 2  ;  ~  ˆ … † 2 RNN ; i; j; k; l ˆ 1; . . . ; N:

(5)

ij

3. Variational formulation and the weak solution of the problem Let us denote by (.,.) the scalar product in ‰L2 … †ŠN ; by k
kk  k
k the norm in Hk( i), where H k …  †; k 2 R1 ;  ˆ 1; . . . ; 4 denotes the Sobolev space in the usual sense. In order to formulate the variational formulation of the problem discussed we introduce the space of virtual displacements by: V ˆ fvjv 2 W… † ˆ ‰H 1 … 1 †ŠN 


 ‰H 1 … 4 †ŠN ; v ˆ 0

on ÿu g

and the set of admissible displacements by: [ ÿkl K ˆ fvjv 2 V; vkn ÿ vln  0 on c g: Let us look for a function u,u 2 [H1( )]N, in case the inertia forces are neglected, i.e. the quasistatic case of the following form: Problem (Pqs). Find u(x,t), t  0 satisfying:   [ @ @uk …x; t†  cijkl …x†

  I; (6) ‡  …x†f  …x; t† ˆ 0; i; j; k; l ˆ 1; . . . ; N; on ÿ @xj @xl  u…x; t† ˆ 0

on ÿu  I;

ij …x; t†nj ˆ Pi …x; t† on ÿ  I;

(7) (8)

ukn …x; t†ÿuln …x; t†  0; nkl …x; t†  nk …x; t† ˆ ÿnl …x; t†  0; …ukn …x; t†ÿuln …x; t††nkl …x; t† ˆ 0; on ÿc I; and kl if ukn …x; t† ÿ uln …x; t† ˆ 0; then j~  kl t …x; t†j  gc …x†; kl kl k l if j~  t …x; t†j < gc …x† then ut …x; t† ÿ ut …x; t† ˆ 0; kl if j~  kl then 9  0 such that ukt …x; t† ÿ ult …x; t† ˆ ÿ…x†~  kl t …x; t†j ˆ gc …x†; t …x; t†; (9) kl 1 where gkl c …x† 2 L …ÿc † denotes the 0given frictional forces. By multiplying Eq. (6) by vi ÿ ui , integrating over   I;  ˆ 1; . . . ; 4, using the Green theorem and boundary conditions, we obtain the variational formulation of the problem (Pqs).

Definition 1. Function u(x,t) is called a weak solution of Eqs. (6)±(9), if for any T > 0, u 2 L2(I;K), …@=@xj …cijkl …x†@uk …x; t†=@xl †† 2 L2 …I; L2 … †† …i ˆ 1; . . . ; N†; @u…x; t†=@t 2 L2 …I; H 1 …  ††; v 2 L1 …I; K†,

288

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

the following inequality: ZT

a…u; v ÿ u0 † dt ‡

ZT

0

…j…v† ÿ j…u0 †† dt 

0

ZT 

ZT

2 3 Z Z 6 7 4 fi …vi ÿ u0i †dx ‡ Pi …vi ÿ u0i † ds5 dt

0

ÿ



S…v ÿ u0 † dt;

(10)

0

S kl kl 1 kl 2 2 holds, where gkl Pi 2 L2 …I; L2 …ÿ ††, a…u; v† ˆ c …x† 2 L …ÿc †, R R gc  0 a.e.Ron ÿc , fi 2 L R…I;SL … ††, kl k l

cijkl eij …u†ekl …v† dx, S…v† ˆ fi vi dx ‡ ÿ Pi vi ds, j…v† ˆ ÿkl gc jvt ÿ vt jdx. c

In order to prove the existence theorem of the weak solution of Eq. (10) the penalty and regularization techniques will be used. The penalized problem to be solved is defined as follows: Problem (Pqsp). Let "m > 0. Find u"i(x,t), t  0 satisfying: [   ÿ@=@xj …cijkl …x†@u" …x; t†=@x † ‡  …x†f …x; t† ˆ 0; i; j; k; l ˆ 1; . . . ; N; on

  I; (11) l k u" …x; t† ˆ 0 on ÿu  I; ij …x; t†nj ˆ Pi …x; t† on ÿ  I;  "k  S "l ‡ on ÿkl nkl …u" …x; t†† ˆ ÿ"ÿ1 c  I; m un ÿ un and "l kl if u"k then j~  kl n …x; t† ÿ un …x; t†† ˆ 0; t …x; t†j  gc …x†; if if

j~  kl t …x; t†j j~  kl t …x; t†j

< ˆ

gkl c …x†; gkl c …x†;

then

u"k t …x; t†

ÿ

u"l t …x; t†



(12) (13)

(14)

ˆ 0;

"l then 9  0 such that u"k  kl t …x; t† ÿ ut …x; t† ˆ ÿ…x†~ t …x; t†:

Then we will prove a priori estimates and finally the limitation process "m ! 0 (i.e. u" ! u weakly). Aweak solution of the penalized problem for "m > 0 fixed can be found on the basis of Theorem 1.2 of [5], (p. 298). The problem (11)±(14) leads to seeking u 2 [H1( )]N such that the following condition is satisfied for any v 2 V where the index " will be omitted, i.e.: Z Z ÿ1 cijkl eij …u†ekl …v† dx ‡ "m ‰ukn ÿ uln Š‡ …vkn ÿ vln † ds ‡ …j0 …u†; v† ˆ S…v†; (15) S kl

ÿc

where t is regarded to be a parameter. The problem (15) leads to minimizing the functional Jm(u) defined by Z Z 1 ÿ1 Jm …u† ˆ cijkl eij …u†ekl …u† dx ‡ …2"m † …‰ukn ÿ uln Š‡ †2 ds ‡ j…u† ÿ S…u†  J0 …u† 2 S

‡ …2"m †ÿ1

ÿkl c

Z S

…‰ukn ÿ uln Š‡ †2 ds ˆ J0 …u† ‡ …2"m †ÿ1 P…u† ÿkl c

R on the space V. The penalty function P…u† ˆ 1=2 S ÿkl …‰ukn ÿ uln Š‡ †2 ds is: c

(16)

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

289

(i) Lipschitz continuous, i.e.: jDP…u; w† ÿ DP…v; w†j  cku ÿ vk1 kwk1 : (ii) Monotone, i.e.: DP…u ‡ w; w† ÿ DP…u; w†  0; DP…u; w† ˆ 0 8w , u 2 K; P…u† ˆ 0 , u 2 K: It is evident that the functional J0 …u† ‡ "ÿ1 m P…u† is coercive and weakly lower semi-continuous. Then the functional Jm(u) for fixed "m attains its minimum over V at the point um. For u 2 K, according to the proof of Theorem 1.2 of [5], (p. 298), we have: ÿ1 DJ0 …um ; um ÿ u† ÿ DJ0 …u; um ÿ u† ‡ "ÿ1 m DP…um ; um ÿ u† ÿ "m DP…u; um ÿ u†

ˆ ÿDJ0 …u; um ÿ u†

(17)

and then: c0 kum ÿ uk2  c1 kum ÿ uk…1 ‡ kuk†: Hence kumk  c2. Then there exists a subsequence {um} such that um ! u0 weakly. Since: DJ0 …um ; w† ‡ "ÿ1 m DP…um ; w† ˆ 0 8wV;

(18)

jDP…um ; w†j  "m c3 kwk;

(19)

then:

then parallel to the proof of Theorem 1.2 of [5], (p. 299) and using the so-called Minty's trick, we will prove that u0 2 K. From (17), putting u ˆ u0 , we find: ckum ÿ u0 k2  ÿDJ0 …u0 ; um ÿ u0 †: Hence um ! u0. Let w 2 K, then: DJ0 …um ; um ÿ w† ‡ "ÿ1 m …DP…um ; um ÿ w† ÿ DP…w; um ÿ w†† ˆ 0: Hence: DJ0 …um ; um ÿ w†  0: For m ! 1 we have: DJ0 …u0 ; u0 ÿ w†  0

8w 2 K;

DJ0 …u0 ; w ÿ u0 †  0

8w 2 K:

or Since (20) has one and only one solution, then um ! u0 is valid for the whole sequence {um}.

(20)

290

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

Thus we have proved that there exists a unique solution u 2 K of the equivalent variational problem DJ…u; v ÿ u†  0 8v 2 K;

(21)

where u ˆ limm!1 um , where um is a solution of (16), t 2 I being fixed. Now we shall prove that u 2 L2(I;K), @u/@t 2 L2(I; [H1( )]N). Let us define the space Ws, s > 0, as:    @ @vk cijkl 2 L2 …I; ‰L2 … †ŠN †; i ˆ 1; . . . ; N; Ws ˆ v 2 L2 …I; ‰H 1 … †ŠN †; @xj @xl )







@v @ @v @v k

s 2 L2 …I; ‰H 1 … †ŠN †;

@x cijkl @x  s; i ˆ 1; . . . ; N; @t 2 @t j l L …I;‰H 1 … †ŠN † Space Ws is a Banach space with the norm kvkWs ˆ kvkL2 …I;‰H 1 … †ŠN † , which follows from the completeness of L2(I;[H1( )]N) with respect to the norm k
kWs and the weak compactness of balls in L2(I;[L2( )]N) or L2(I;[H1( )]N), respectively (see [9]). In order to prove that u 2 Ws we will estimate u and @u/@t. We shall use the following regularized formulation [2]: "

a…u ; v† ‡

…j0" …u" †; v†

Z ÿ

S

cijkl ÿkl c

@u"k vi nj ds ˆ S…v† @xl

8v 2 V;

(22)

R ÿ1 …1‡"† " k l " where we put j" …v† ˆ S ÿkl gkl ; " > 0. Let us put c …x†' …vt ÿ vt † ds; ' …† ˆ …1 ‡ "† jj c 0 0 v ˆ u" …t† (where we denoted u" …t†  @t u" …t††=@t, integrate over t and since … j0 …v†; v† > 0, then: 1 a…u" …t†; u" …t††  2 Zt ‡

0

Zt 0

Z S

cijkl ÿkl c

@u"k "0 1 ui nj ds dt ‡ a…u" …0†; u" …0†† 2 @xl

0

Z ÿ

S

ÿkl c

Z

0

‰…f; u" …t†† ‡ …P; u" …t††Š dt 

@u" …0† 0 cijkl k u"i …0†nj ds ‡ @xl

S

Zt

cijkl ÿkl c 0

@u"k "0 u nj ds @xl i 0

‰…f; u" …t†† ‡ …P; u" …t††Š dt

0

Let us estimate the first term on the right-hand side. We put vh …:; t† ˆ hÿ1 …v…:; t ‡ h† ÿ v…:; St††; ÿ‡ ˆ fx 2 ÿkl u …:; t† ˆ hÿ1 …u…:; t ‡ h† ÿ u…:; t††; h > 0; uh …:; t ‡ h†  u1 …:; t†; uh …:; t†  u2 …:; t†, c juk1n ÿ ul1n > 0; uk2n ÿ ul2n > 0g. Then, according to [9]) we have: Z h h a…u ; u † ÿ cijkl @uhk =@xl uhi nj ds ˆ S…uh †: S kl h

ÿc

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

Thus: Z S

Z @uhk h cijkl u nj ds ˆ @xl i S

ÿkl c

Z

ÿ

S

…n1 …uk2n

Z ˆ

Then: S

S

T

ÿkl c

ÿc

ÿ‡

T

Z S

‡ ÿkl c jÿ

Z S

n2 …uk1n

ÿ

ul1n †† ds

Z 1 ˆÿ "m S



 j‰uk1n ÿ ul1n Š‡ j2 ‡ j‰uk2n ÿ ul2n Š‡ j2 ds

ÿkl c

Z . . . ds ‡

ÿ‡

1 . . . ds ˆ ÿ "m

. . . ds ˆ ÿ ÿ‡

S

. . . ds ‡ ÿkl c nÿ

Z

 k 2 …u1n ÿ ul1n † ÿ …uk2n ÿ ul2n † ds  0

ÿ‡

1 "m S

Z ‡ ÿkl c jÿ



 1 j‰uk1n ÿ ul1n Š‡ j2 ‡ j‰uk2n ÿ ul2n Š‡ j2 ds ‡ "m

ÿ k
1 ‰u1n ÿ ul1n Š‡ …uk2n ÿ ul2n † ‡ ‰…uk2n ÿ ul2n Š‡ …uk1n ÿ ul1n †† ds  ÿ "m 

 ‰uk1n ÿ ul1n Š‡ j2 ‡ j‰uk2n ÿ ul2n Š‡ j2 ds  0

‡ ÿkl c jÿ

Hence: Z S

‡

ÿkl c

Z S

ÿkl c

…‰uk1n ÿ ul1n Š‡ …uk2n ÿ ul2n † ‡ ‰uk2n ÿ ul2n Š‡ …uk1n ÿ ul1n †† ds

T kl

Z

ÿkl c

…n1 …uk1n ÿ ul1n † ‡ n2 …uk2n ÿ ul2n †† ds

ÿkl c

Z 1 "m S

‡

ÿ

ul2n †

291

cijkl ÿkl c

@uhk h u nj ds  0: @xl i

(24)

The second term ˆ 0 as u" …0† ˆ 0. From Eqs. (23) and (24), assuming the time differentiability of f, P and using the integration by parts, we obtain: Zt 2 " ku …t†k  c ku" …t†k2 dt ‡ ju…0†j2 0

292

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

and thus using the Gronwall lemma we obtain: ku" …t†k2  c;

(25)

0

i.e. u" is bounded in L2(I;H1( )]N) and …@=@xj …cijkl @u"k =@xl †† …i ˆ 1; . . . ; N† is bounded in L2(I;[L2( )]N). 00 Now we shall differentiate Eq. (22) with respect to t and let v ˆ u" …t†. Let us put X…t† ˆ 0 …j0" …u" …t†; u" …t†††. Thus:   Z d 0 d 0 k l k l j" …w†; v ˆ gkl c '" …wt …t† ÿ wt …t††…vt ÿ vt † ds dt dt S kl Z

ˆ

d dt S Z

ˆ

S

ÿc

k l k l gkl c " …wt …t† ÿ wt …t††…vt ÿ vt † ds ÿkl c ÿ1 k l k l k l gkl c lim h ‰ " …wt …t ‡ h† ÿ wt …t ‡ h†† ÿ " …wt …t† ÿ wt …t††Š…vt ÿ vt † ds h!0

ÿkl c

Therefore:   Z d 0 j" …w…t††; w…t† ˆ dt S

ÿ1 k l k gkl c lim h ‰ " …wt …t ‡ h† ÿ wt …t ‡ h†† ÿ " …wt …t† h!0

ÿkl c

ÿ wlt …t††Šhÿ1 …wkt …t ‡ h† ÿ wlt † ds: With respect to the monotony property of the functional   d 0 j …w…t††; w…t†  0: dt " Then we have: "0

Z

"00

a…u …t†; u …t†† ÿ

S

"

the last integral 0, and therefore,

0

@u" …t† 00 00 00 cijkl k u"i …t†nj ds  …f 0 ; u" …t†† ‡ …P0 ; u" …t††: @xl

ÿkl c

Integrating with respect to t, we obtain: "0

2



2

"0

ku …t†k  c ku …0†k



Zt   0 ‡ ku" …t†k2 ‡ ju" 0…t†j2 dt; 0

where the second integral on the left-hand side can be estimated as in the previous case. Using the Gronwall lemma we obtain: 0

ku" …t†k2  c; "0

(26) 2

1

N

i.e. u (t) is bounded in L (I;[H ( )] ).

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

293

With respect to estimates (25) and (26), we can select such subsequence, we also denote it as u", such that: u" ! u weakly-star in L1 …I; V†; u" ! u" weakly in L1 …I; ‰H 1 … †ŠN †; 0

…cijkl @u"k …t†=@xl † ! …cijkl @uk …t†=@xl † …i ˆ 1; . . . ; N† weakly in L2 …I; ‰L2 … †ŠN †: With respect to Eq. (22) we have "0

"

a…u ; v ÿ u † ‡

…j0" …u" †; v

Z

"0

ÿu †ÿ

S

@u"k …t† 0 …vi ÿ u"i †nj ds ˆ S…v ÿ u" † @xl

cijkl ÿkl c

8vV;

hence: 0

Z

0

a…u" ; v ÿ u" † ‡ j" …v† ÿ j" …u" † ÿ

S

0

cijkl ÿkl c

@u"k …t† 0 0 …vi ÿ u"i †nj ds ÿ S…v ÿ u" † @xl

0

ˆ j" …v† ÿ j" …u" † ÿ …j0" …u" †; v ÿ u" †  0:

S Thus, for v ˆ v…t†; v 2 L2 …I; V†; vkn ÿ vln  0 on ÿkl c we obtain 3 2 Z ZT 7 6 " @u"k …t† 6a…u ; v† ‡ j" …v† ÿ S…v† ÿ c vi nj ds7 ijkl 5 dt 4 @xl S kl 0 ZT 

ÿc

2

6 " "0 6a…u ; u † ‡ j" …u"0 † ÿ 4

0

"

ZT

"

‡ a…u …0†; u …0††Š ‡

S

cijkl ÿkl c

"0

j" …u † dt ÿ 0

ZT

"0

ZT

j" …u …t†† dt ˆ

‡ 0

3

Z

7 @u"k …t† "0 1 ui nj ds7 dt ˆ ‰a…u" …T†; u" …T†† 5 @xl 2

Z S

cijkl ÿkl c

@u"k …T† " 1 ui …T†nj ds  a…u" …T†; u" …T†† @xl 2

a…u; u0 † ‡ j" …u0 †† dt

(27)

0

Then from (27) we obtain: ZT

‰a…u; v ÿ u0 † ‡ j…v† ÿ j…u0 † ÿ S…v ÿ u0 †Š dt  0

0

which conclude the validity of Eq. (10).

8v 2 L2 …I; V†; vkn ÿ vln  0 on

[

ÿkl c;

294

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

We have further proved the existence of a solution to the quasi-static contact problem. The uniqueness of the solution can be proved in the usual way (see [1]). We have proved the following result: Theorem 1. Let N ˆ 2; 3. Let us assume that: f; f 0 2 L2 …I; ‰L2 … †ŠN ; P; P0 2 L2 …I; ‰L2 …ÿ †ŠN †: Then there exists a solution u to the problem (Pqs) in the sense of Definition 1. 4. Analysis of numerical results Classical biomechanics is based on the vector analysis of the loaded human skeleton and modern geometry methods. Such approaches cannot analyse the loosening of the THRs. The real situation of the loosening of the THRs is very complicated and therefore must be simulated by a model problem discussed in Section 2 and based on the variational formulation of Section 3 and its finite element approximation. We must say that the loosening of the THR is a model problem at the time when in the first approximation the inertia forces can be omitted. In our model problem discussed a real patient was analysed. Fig. 1(a)±(e) shows and analyses the deformation of a part of his skeleton and the distribution of stresses in this part of the skeleton when the THR is not loosened. In Figs. 2±4 the loosened THR is analysed numerically; the case in the Fig. 2

Fig. 1. The unloosened total hip replacement (THR): (a) the deformation of THR; (b)±(e) the vertical stress component and the principal stresses and their zooming.

Fig. 1. (Continued )

296

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

Fig. 1. (Continued )

corresponds to a real situation of the patient investigated, while the last cases correspond to a future situation of loosening the THR. Fig. 1(a) shows in detail that the skeleton with the THR is only deformed, while the detailed Fig. 3(a), Fig. 4(a) show that the stem is loosened and bulged out

Fig. 2. The loosened THR ± the real stage: (a)±(d) the vertical stress component and the principal stresses and their zooming.

298

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

Fig. 2. (Continued )

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

299

Fig. 3. The loosened THR ± the future stage: (a) the deformation of THR; (b)±(e) the vertical stress component and the principal stresses and their zooming.

300

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

Fig. 3. (Continued )

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

301

Fig. 3. (Continued )

Fig. 4. The loosened THR ± the future stage: (a) The normal and tangential components of displacement; (b) The normal and tangential components of stresses.

302

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

Fig. 4. (Continued )

Fig. 5. Convergence analysis of preconditioned conjugate gradient method: (a) the Jacobi preconditioner; (b) the SOR preconditioner; (c) the ILL preconditioner.

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

303

Fig. 5. (Continued )

approximately by about 0.01±0.06 mm. Fig. 1(b)±(e),2(a)±(d),3(b)±(e) show the vertical stress component and the principal stresses and their zooming. In Fig. 2(a),(b) we can see that the lines of constant stresses ± the isostresses ± between values ``G'' and ``H'' denote the lines between

304

J. Nedoma et al. / Mathematics and Computers in Simulation 50 (1999) 285±304

compression and extension (and similarly in other cases). It is evident that the left parts of the THR are only compressed while in the right parts of the stem extensions occur. This last character of the stress distribution is evident from the principal stresses (see Fig. 1 (d),(e),2(d),3(d),(e)). In Fig. 1(d),(e) we can find relatively great extensions in the trochanteric part of the femur near the collar of the THR, which become weak with enlarged loosened domain, which evident from Fig. 3(d),(e). Greater values of extension change their values to the places below the loosened parts of the THRs stem. Although the principal values of compression and/or extension are dominant in one direction, the extension and/or compression in the second direction of the principal stresses has a relatively great value. When stresses in the area of the trochanteric part of the femur bone vanish, the bone there starts to lose robustness and the bone structure gets thinner (more porous). The numerical results indicate that this osteoporosis is the result of the loosened THR. In Fig. 5(a)±(c) the convergence of preconditioned conjugate gradient method (by Jacobi, SOR and ILL preconditioners) are given. References [1] K.A. Anes, L.E. Payne, Uniqueness and continuous dependence of solutions to a multi-dimensional thermo-elastic contact problem, J. Elasticity 34 (1994) 139±148. [2] G. Duvaut, J.L. Lions, Les ineÂquations en MeÂcanique et en Physique, Dunod, Paris, 1972. [3] Eck Ch., Existenz und RegularitaÈt der LoÈsungen fuÈr Kontaktprobleme mit Reibung, Thesis, Math. Institut A der UniversitaÈt Stuttgart, 1996, 168 pp. [4] J. Haslinger, I. HlavaÂcÏek, Contact between two elastic bodies. I. Continuous problems, Aplikace matematiky 25 (1978) 324±347. [5] J. NecÏas, I. HlavaÂcÏek, Mathematical Theory of Elastic and Elasto-Plastic Bodies: An Introduction, Elsevier, Amsterdam, 1981. [6] J. Nedoma, On the Signorini problem with friction in linear thermo-elasticity, the quasi-coupled 2D-case, Appl. Math. 32 (1987) 186±199. [7] J. Nedoma, Mathematical modelling in biomechanics. Bone- and vascular-implant systems, Habilitation Thesis, Institute of Computer Science ASCR, Prague, 1993, 220 pp. [8] J. Nedoma, Contact Problem in Biomechanics and Iterative Solution Methods for Constrained Optimization Theory, Technical Report no. V-756, ICS AS CR, Prague, 1998. [9] J.E. Rivera, R. Racke, Multi-dimensional Contact Problems in Thermo-elasticity, Technical Report no. 26/97, Lab. Nacional de Comput. Cientifca-LNCC, PetroÂpolis, Brazil, 1997. [10] J. StehlõÂk, J. Nedoma, Mathematical Simulation of Function of Great Human Joints and Optimal Design of Their Artificial Substitutes, Parts I,II, Research Reports V-406,V-407, Institute of Computer Science ASCR, Prague (in Czech), 1989, 1990.