- Email: [email protected]
Boundary element analysis of stress distribution around a crack in plane micropolar elasticity
International Journal of Engineering Science 45 (2007) 199–209 www.elsevier.com/locate/ijengsci
Boundary element analysis of stress distribution around a crack in plane micropolar elasticity E. Shmoylova, S. Potapenko *, L. Rothenburg Department of Civil Engineering, University of Waterloo, 200 University Avenue West, Waterloo, Ontario, Canada N2L 3G1 Received 20 October 2006; received in revised form 22 November 2006; accepted 30 November 2006 Available online 7 June 2007
Abstract In this paper, we use the boundary element method to find a semi-analytical solution to the problem of stress concentration around a crack in plane micropolar elasticity. We provide an example demonstrating the effect of material microstructure. 2007 Elsevier Ltd. All rights reserved. Keywords: Micropolar elasticity; Weak solutions; Boundary element method; Cracks in micropolar medium
1. Introduction The theory of micropolar elasticity (also known as Cosserat or asymmetric theory of elasticity) was introduced by Eringen [1] (see [2] for a review of works in this area and an extensive bibliography) to eliminate discrepancies between classical theory of elasticity and experiments in cases when effects of material microstructure were known to contribute significantly to the body’s overall deformation, for example, materials with granular microstructure such as polymers or human bones (see [3–6]). These cases are becoming increasingly important in the design and manufacture of modern day advanced materials as small-scale effects become very important in the prediction of the overall mechanical behavior of these materials. For the last 30 years numerous investigations have been performed in this direction. For example, threedimensional problems of Cosserat elasticity have been formulated in a rigorous setting and solved by means of potential theory methods by Kupradze [7]. The corresponding boundary value problems for plane and antiplane shear deformations of a micropolar homogeneous, linearly elastic solid were shown to be well-posed and subsequently solved using the boundary integral equation method by Iesan [8], Schiavone [9], Potapenko et al. [10,11] and Potapenko [12,13].
*
Corresponding author. Tel.: +1 519 888 4567x37156; fax: +1 519 888 6197. E-mail addresses: [email protected] (E. Shmoylova), [email protected] (S. Potapenko), [email protected] (L. Rothenburg). URL: www.civil.uwaterloo.ca/spotapen (E. Shmoylova). 0020-7225/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2007.04.006
200
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Several studies relating to investigations of stress distributions around a crack have been undertaken under assumptions of a simplified theory of plane Cosserat elasticity by Mu¨hlhaus and Pasternak [14], Atkinson and Leppington [15], experimentally by Lakes and Nakamura [5], and using the finite element method (see [6]). There also has been some activity lately in the area of crack analysis in the three-dimensional Cosserat elasticity (see [16–19]). Recently, Chudinovich and Constanda [20] used the boundary integral equation method in a weak (Sobolev) space setting to obtain the solution for fundamental boundary value problems in a theory of bending of classical elastic plates. This approach has wide practical applicability because it also covers domains with reduced boundary smoothness. In addition, it allows to give an answer for the fundamental question relating to the existence and uniqueness of the solution and, second, gives an opportunity to employ an effective numerical procedure which allows the construction of a numerical solution which can be very useful for practical purposes. Further, Chudinovich and Constanda [20] extended their method to accommodate several problems relating to the investigation of stress concentrations around a crack in classical plates. In a number of very recent works by Shmoylova et al. [21–23], the authors performed the rigorous analysis of interior and exterior Dirichlet and Neumann boundary value problems in plane Cosserat elasticity and also obtained the solution to the crack problem arising in this theory in the form of modified integral potentials with unknown distributional densities. Unfortunately, it is very difficult, if not impossible, to find these densities analytically, consequently, we have to find a numerical technique which will allow us to obtain a numerical approximation of the solution. One of the most effective approaches to achieve this goal is to use the boundary element method. This method has been developed by Brebbia [24] and has become very popular among researchers in different areas including fracture mechanics (see, for example, [25] for references on applications of the boundary element method in science and engineering). The boundary element method has originated from works on classical integral equations and finite elements and incorporates advantages of both techniques. On the one hand, it allows us to reduce the dimension of a problem by one and defines domains extending to infinity with a high degree of accuracy exactly as the boundary integral equation method. One the other hand, the boundary element method does not require the differentiation of shape functions, which is the major requirement of the finite element method when we have to find stresses, but allows us to differentiate the matrix of fundamental solutions instead, which makes the calculation of stresses easier and more accurate. In this paper, we use the boundary element method to find the solution for an infinite domain weakened by a crack in plane Cosserat elasticity, when stresses and couple stresses are prescribed along both sides of the crack (Neumann boundary value problem), and discuss its convergence. To illustrate the effectiveness of the method for applications we consider a crack in a human bone which is modelled under assumptions of plane micropolar elasticity. We find the numerical solution for stresses around the crack and show that the solution may be reduced to the classical one if we set all micropolar elastic constants equal to zero. We come to the conclusion that there could be up to 26% difference in quantitative characteristics of the stress around a crack in the micropolar case in comparison with the model when microstructure is ignored (classical case, see for example, [26]).
2. Preliminaries In what follows Greek and Latin indices take the values 1, 2 and 1, 2, 3, respectively, the convention of summation over repeated indices is understood, Mmn is the space of ðm nÞ-matrices, En is the identity element in Mnn , the columns of a (3 · 3)-matrix P are denoted by P ðiÞ , a superscript T indicates matrix transposition, the generic symbol c denotes various strictly positive constants, and ð Þ;a oð Þ=oxa . Let S be a domain in R2 occupied by a homogeneous and isotropic linearly elastic micropolar material with elastic constants k, l, a, c and e. We use the notations k k0;S and h; i0;S for the norm and inner product in L2 ðSÞ \ Mm1 for any m 2 N. When S ¼ R2 , we write k k0 and h; i0 . T The state of plane micropolar strain is characterized by a displacement field uðx0 Þ ¼ ðu1 ðx0 Þ; u2 ðx0 Þ; u3 ðx0 ÞÞ T and a microrotation field /ðx0 Þ ¼ ð/1 ðx0 Þ; /2 ðx0 Þ; /3 ðx0 ÞÞ of the form
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
ua ðx0 Þ ¼ ua ðxÞ; u3 ðx0 Þ ¼ 0; /a ðx0 Þ ¼ 0; /3 ðx0 Þ ¼ /3 ðxÞ;
201
ð1Þ
where x0 ¼ ðx1 ; x2 ; x3 Þ and x ¼ ðx1 ; x2 Þ are generic points in R3 and R2 , respectively. The equilibrium equations of plane micropolar strain written in terms of displacements and microrotations are given by [8,9] Lðox ÞuðxÞ ¼ qðxÞ;
x 2 S;
ð2Þ T
in which now, denoting /3 by u3 , we have uðxÞ ¼ ðu1 ; u2 ; u3 Þ , the matrix partial differential operator Lðox Þ ¼ Lðo=oxa Þ is defined by 0 1 ðk þ l aÞn1 n2 2an2 ðl þ aÞD þ ðk þ l aÞn21 B C LðnÞ ¼ Lðna Þ ¼ @ A; ðk þ l aÞn1 n2 ðl þ aÞD þ ðk þ l aÞn22 2an1 2an2
2an1
ðc þ eÞD 4a
T
where D ¼ na na , and the vector q ¼ ðq1 ; q2 ; q3 Þ represents body forces and body couples. Together with L, we consider the boundary stress operator T ðox Þ ¼ T ðo=oxa Þ defined by 0 1 ðl aÞn1 n2 þ kn2 n1 2an2 ðk þ 2lÞn1 n1 þ ðl þ aÞn2 n2 B C T ðnÞ ¼ T ðna Þ ¼ @ ðl aÞn2 n1 þ kn1 n2 ðk þ 2lÞn2 n2 þ ðl þ aÞn1 n1 2an1 A; 0 0 ðc þ eÞna na T
where n ¼ ðn1 ; n2 Þ is the unit outward normal to oS. To guarantee the ellipticity of system (2), in what follows we assume that k þ l > 0;
l > 0;
c þ e > 0;
a > 0:
The internal energy density is given by 2Eðu; vÞ ¼ 2E0 ðu; vÞ þ lðu1;2 þ u2;1 Þðv1;2 þ v2;1 Þ þ aðu1;2 u2;1 þ 2u3 Þðv1;2 v2;1 þ 2v3 Þ þ ðc þ eÞðu3;1 v3;1 þ u3;2 v3;2 Þ; 2E0 ðu; vÞ ¼ ðk þ 2lÞðu1;1 v1;1 þ u2;2 v2;2 Þ þ kðu1;1 v2;2 þ u2;2 v1;1 Þ: Clearly, Eðu; uÞ is a positive quadratic form. The space of rigid displacements and microrotations F is spanned by columns of the matrix 0 1 1 0 x2 B C F ¼ @ 0 1 x1 A 0 0 1 from which it can be seen that LF ¼ 0 in R2 , T F ¼ 0 on oS and a general rigid displacement can be written as Fk, where k 2 M31 is constant and arbitrary. A Galerkin representation for the solution of (2) when qðxÞ ¼ dðjx yjÞ, where d is the Dirac delta distribution, yields the matrix of fundamental solutions (see [9]) Dðx; yÞ ¼ L ðoxÞtðx; yÞ;
ð3Þ
where L is the adjoint of L a 2 tðx; yÞ ¼ f½k 2 jx yj þ 4 ln jx yj þ 4K 0 ðkjx yjÞg: 4 8pk K 0 is the modified Bessel function of order zero and the constants a; k 2 are defined by a1 ¼ ðc þ eÞðk þ 2lÞðl þ aÞ; In view of (3) and (4) Dðx; yÞ ¼ DT ðx; yÞ ¼ Dðy; xÞ:
k2 ¼
4la : ðc þ eÞðl þ aÞ
ð4Þ
202
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Along with matrix Dðx; yÞ we consider the matrix of singular solutions T
P ðx; yÞ ¼ ðT ðoyÞDðy; xÞÞ : ðiÞ
ð5Þ ðiÞ
2
It is easy to verify that D ðx; yÞ and P ðx; yÞ satisfy (2) with qðxÞ ¼ 0 at all x 2 R ; x 6¼ y. We consider an infinite domain with a crack modelled by an open arc C0 and assume that C0 is a part of a simple closed C 2 -curve C that divides R2 into interior and exterior domains Xþ and X . In what follows we denote by the superscripts + and the limiting values of functions as x ! C from within Xþ or X . Further, we define X ¼ R2 n C0 and C1 ¼ C n C0 . Let H 1;x ðXÞ be the space of all u ¼ fuþ ; u g such that uþ 2 H 1 ðXþ Þ, ± u 2 H 1;x ðX Þ and cþ 1 uþ ¼ c1 u , where c0 are the trace operators on C1 from within X , respectively. The þ space H 1 ðX Þ is a standard Sobolev space and H 1;x ðX Þ is a weighted Sobolev space defined in [21]. Further, we introduce the corresponding single layer and double layer potentials given respectively by Z ðV uÞðxÞ ¼ Dðx; yÞuðyÞdsðyÞ; C Z0 P ðx; yÞuðyÞdsðyÞ; ðW uÞðxÞ ¼ C0
where u 2 M31 is an unknown density matrix. 3. Boundary value problem Let us consider Neumann boundary value problem. We seek u 2 C 2 ðXÞ \ C 1 ðXÞ, u 2 A such that LuðxÞ ¼ 0; þ
x 2 X; þ
ðTuÞ ðxÞ ¼ g ðxÞ;
ðTuÞ ðxÞ ¼ g ðxÞ;
x 2 C0 ;
ðNÞ
where gþ and g are prescribed on C0 . Asymptotic class A was introduced in [9]. The variational formulation of problem (N) is as follows. We seek u 2 H 1;x ðXÞ such that 8v 2 H 1;x ðXÞ; bðu; vÞ ¼ hdg; cþ 0 vþ i0;C0 þ hg ; dvi0;C0 R R where bðu; vÞ ¼ 2 X Eðu ; v Þdx þ 2 Xþ Eðuþ ; vþ Þdx. Note that (6) is solvable only if
hz; dgi0;C0 ¼ 0
8z 2 F:
ð6Þ
ð7Þ
± We introduce the restriction operators p± to X± and the trace operators c 0 on C0 from within X , respectively and define the modified single layer potential V of density u by ðVuÞðxÞ ¼ ðV uÞðxÞ ðV uÞ0 ; ~zðiÞ 0;C ~zðiÞ ðxÞ; x 2 R2 ; 0
where V u is the single layer potential, and V0 is the boundary operator defined by ðV uÞ0 ¼ c 0 p V u, and 3 f~zðiÞ gi¼1 is an L2 ðC0 Þ-orthonormal basis for F. Also we introduce the modified double layer potential W of density w ðWwÞðxÞ ¼ ðW wÞðxÞ p0 W þ w; ~zðiÞ 0;C ~zðiÞ ðxÞ; x 2 X; 0
where p0 is the restriction operator to C0. Solution of problem (6) may be represented in the form u ¼ ðVuÞX þ Ww þ z;
ð8Þ
where u and w are unknown densities, and z 2 F is arbitrary. The detailed procedure of obtaining solution (8) has been developed by Shmoylova et al. in [23]. 4. Boundary element method Consider problem (6). As shown in [23], the solution to this problem may be represented in the form (8) and the corresponding boundary integral equations are uniquely solvable with respect to distributional densities u
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
203
and w. As we stated above, these densities cannot be found analytically. To approximate them numerically we use the boundary element method [27], which makes use of the following classical result. Lemma 1 (Somigliana formula). Using classical techniques, it can be proved that if u 2 H 1;x ðXÞ is a solution of Lu = 0 in X, then Z
1 ½Dðx; yÞdðT ðoy ÞuðyÞÞ P ðx; yÞduðyÞdsðyÞ ¼ duðxÞ; 2 C0
x 2 C0 ;
ð9Þ
where dðT ðoy ÞuðyÞÞ denotes the jump of T ðoy ÞuðyÞ on the crack. It has been shown in [23], that the density of the modified single layer potential may be found in the form u ¼ dðT ðoy ÞuðyÞÞ ¼ dg. Now we need to find the density of the modified double layer potential w ¼ du. To ðkÞ achieve this goal we divide C0 into n elements C0 , each of which possesses one node nðkÞ located in the middle of the element. The values of dg and du are constant throughout the element and correspond to the value at the node dgðnðkÞ Þ and duðnðkÞ Þ. Then (9) becomes n Z X 1 ½Dðx; yÞdgðnðkÞ Þ P ðx; yÞduðnðkÞ ÞdsðyÞ ¼ duðxÞ; x 2 C0 : ðkÞ 2 C0 k¼1 Placing x sequentially at all nodes, we obtain the linear algebraic system of equations ! ! Z Z n n X X 1 ðiÞ ðkÞ ðiÞ Dðn ; yÞdsðyÞ dgðn Þ P ðn ; yÞdsðyÞ duðnðkÞ Þ ¼ duðnðiÞ Þ; i; k ¼ 1; n ðkÞ ðkÞ 2 C0 C0 k¼1 k¼1
ð10Þ
with respect to duðn R i Þ. We note that CðkÞ DðnðiÞ ; yÞdsðyÞ are defined for any i and k, as in [9]. 0 Solving (10) we construct the approximation to w. If we introduce the shape function Uk ðxÞ by ( ðkÞ 1; x 2 C0 ; Uk ðxÞ ¼ ðkÞ 0; x 2 C0 n C0 ; Pn Pn then the approximated densities u and w are uðnÞ ðxÞ ¼ k¼1 Uk ðxÞdgðnðkÞ Þ and wðnÞ ðxÞ ¼ k¼1 Uk ðxÞduðnðkÞ Þ and the approximate solution is becoming of the form uðnÞ ¼ ðVuðnÞ ÞX þ WwðnÞ þ z, where z is arbitrary. We now have to prove that the approximate solution uðnÞ will converge to the exact analytical solution u when n ! 1. Let us formulate the following theorem. Theorem 2. uðnÞ ! u as n ! 1. Proof. Since we consider the Neumann problem, rigid displacement terms are not determined. Consequently, it is enough to show that V uðnÞ ! V u and W wðnÞ ! W w as n ! 1. Consider V uðnÞ . For x 2 X ! Z 3 Z n X X ðkÞ ðnÞ ðiÞ ðiÞ jV uðxÞ V u ðxÞj 6 D ðx; yÞuðyÞdsðyÞ D ðx; yÞdsðyÞ uðn Þ ðkÞ C0 C0 i¼1 k¼1 3 X n Z X ¼ ½DðiÞ ðx; yÞuðyÞ DðiÞ ðx; yÞuðnðkÞ ÞdsðyÞ ðkÞ C0 i¼1 k¼1 3 X n Z X ¼ DðiÞ ðx; yÞ½uðyÞ uðnðkÞ ÞdsðyÞ ðkÞ C i¼1 k¼1 0
6
3 X
n X
i¼1
k¼1
kDðiÞ ðx; Þk1;x;X juðyÞ uðnðkÞ Þjhk ;
where hk is the length of the kth element assuming that all hk ¼ h ¼ Ln, L is the length of C0 , and P3 P3 P ðkÞ ðkÞ juðyÞ uðn Þj 6 i¼1 jui ðyÞ ui ðn Þj ¼ i¼1 2a¼1 joa ui ðnðkÞ Þjh þ Oðh2 Þ.
204
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Denote by M 1 ¼ maxa¼1;2;i¼1;3;k¼1;n joa ui ðnðkÞ Þj. Since kDðiÞ ðx; Þk1;x;X is uniformly bounded [8,9], i.e. there exists M 2 > 0 such that kDðiÞ ðx; Þk1;x;X 6 M 2 for any x 2 X, we obtain that jV uðxÞ V uðnÞ ðxÞj 6
18M 1 M 2 L2 ! 0 as n ! 1: n
Repeating the proof for W wðnÞ , we conclude that uðnÞ ! u as n ! 1. h 5. Example As an example, we consider a longitudinal crack inside a human bone in the case when constant normal stretching pressure of magnitude p is applied on both sides of the crack. If we consider a typical transversal cross-section of the bone and assume that this cross-section is small enough then the deformation of each cross-section under the prescribed load will be the same throughout the length of the bone and will develop in the plane of the cross-section. Consequently, such deformations may be considered under assumptions of plane micropolar elasticity. Such a model is not an idealization that lies far from reality, as it may seem at first, but, as shown, for example, in [28–30], can describe actual cracks in bones very closely, since orthopedic biomechanics usually deals with cracks of a very small size. We model a crack as an open arc of the circle given by equations: x1 ¼ a cos h, x2 ¼ a sin h, h 2 ð0; p=6Þ (Fig. 1). By changing the radius a of the circle we will change the length of the crack. We are interested in how the normal traction distributes at a distance from the crack tip along the line: x1 ¼ a, x2 < 0. Clearly, this problem can be considered as the Neumann problem described above. Elastic constants for a human bone have been measured in [5] and take the following values: a = 4000 MPa, c = 193.6 N, e = 3047 N, k = 5332 GPa, l = 4000 MPa. In our example we construct solutions for cracks of lengths equal to 0.26 mm, 0.52 mm, 0.75 mm, and 10 mm to show good agreement of our results with those presented in the experimental study by Nakamura and Lakes [5] performed on human bone cracks of same lengths. We also assume that the normal stretching pressure p takes a value of 2 MPa. Let the distance from the tip of the crack be q ¼ jx2 j. The numerical solution for boundary tractions and moments has been found to approximate the exact solution to five decimal places for n = 52 elements of C0 (see Table 1). Let us now compare the results for the normal traction in the micropolar case with the results of the classical theory. The classical case may be obtained from the solution for micropolar elasticity setting the micropolar elastic constants equal to zero. In Figs. 2–5 there is a graphical representation for the distribution of the normal traction at a distance from the lower crack tip for crack lengths equal to 0.26 mm, 0.52 mm, 0.75 mm
Fig. 1. Crack in x1 x2 -plane. Length: 0.5 mm – solid line (a = 1 mm).
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
205
Table 1 Approximate solution for a crack 0.1
0.5
0.7
1
n=4
Distance q (mm) Tn (MPa) Ts (MPa) M 3 (N/m)
0.983256 1.175489 162.5682
0.108643 0.538926 91.24553
0.082351 0.213505 79.45362
0.057282 0.756437 70.64523
n = 10
Tn (MPa) Ts (MPa) M 3 (N/m)
0.569361 0.634936 84.24634
0.068735 0.264282 49.86301
0.052678 0.091475 44.26856
0.032536 0.045343 36.09357
n = 30
Tn (MPa) Ts (MPa) M 3 (N/m)
0.427549 0.506874 69.24764
0.061862 0.184756 34.25447
0.039754 0.075982 30.25879
0.021830 0.033547 24.62579
n = 50
Tn (MPa) Ts (MPa) M 3 (N/m)
0.398462 0.454906 61.11549
0.053918 0.108063 29.69127
0.030028 0.061039 22.59611
0.015244 0.026688 15.80458
n = 52
Tn (MPa) Ts (MPa) M 3 (N/m)
0.398456 0.454899 61.11542
0.053910 0.108054 29.69119
0.030021 0.061033 22.59604
0.015239 0.026678 15.80451
Length: 0.52 mm. Tn – normal traction, Ts – tangential traction, M 3 – moment about x3-axis.
Fig. 2. Normal traction on the edge of the crack. Length: 0.26 mm, (a) micropolar, (b) classical.
and 10 mm correspondingly. The traction is divided by the applied load p to represent the data in non-dimensional values. The bold curve characterizes the stress distribution in the micropolar case while the classical case is plotted by the ordinary curve. The distance between the first point, in which we compute the normal traction, and the tip of the crack is equal to one fifth of the length of the crack. It may be seen that the normal traction is significantly higher in the vicinity of the crack tip in the micropolar case in comparison with the case when microstructure is ignored (classical theory), particularly, for cracks of 0.26 mm, 0.52 mm and 0.75 mm long. In the case when the length of the crack is equal to 0.75 mm we can observe that the normal traction in the vicinity of the crack tip is 26.8% higher under
206
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Fig. 3. Normal traction on the edge of the crack. Length: 0.52, (a) micropolar, (b) classical.
Fig. 4. Normal traction on the edge of the crack. Length: 0.75, (a) micropolar, (b) classical.
assumptions of Cosserat elasticity in comparison with the classical case. At the same time, for the crack of the linear size equal to 10 mm, the traction in the vicinity of the crack tip differs only by 13%. The fact that this difference is still present and significantly lower than in the case of the shorter cracks may be explained by the size effect of the crack opening, which is of the order of 105 for this crack in comparison with that of 106 in the other cases. When it comes to the consideration of stresses at a distance from the crack tip, we can conclude that the traction in the micropolar case decays faster than in the classical case and approximately at a distance of one
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
207
Fig. 5. Normal traction on the edge of the crack. Length: 10 mm, (a) micropolar, (b) classical.
crack length the values of the normal traction in both cases become equal to each other. Farther from the crack tip, the traction in the micropolar case becomes lower than in the classical case, especially when we consider the crack of 0.26 mm long for which, as it may be seen from Fig. 2, the difference is drastic and may be up to 19.8%. Further, it may be observed that at a distance of approximately three crack lengths from the tip of the crack the effect of the crack on stresses is negligible in agreement with Saint-Venant’s principle. The results presented in Figs. 2–5 may be compared with earlier investigations undertaken in the area of Cosserat solids. Lakes and Nakamura [5] performed an experimental study on stress concentrations around cracks in a human bone. They considered the same crack lengths as in the present article. Later, the same authors validated their results presented in [5] using finite element method [6]. However, it should be noted, that direct comparison with the results presented in [5,6] does not seem feasible since our formulation of the crack problem is different from that adopted in [5,6], where a crack is considered as a ‘blunt’ notch similarly to the approach of classical fracture mechanics [26,31,32], in which a crack is usually modelled as a ‘squashed’ ellipse of small eccentricity. The tip of the crack under consideration in [5,6], therefore, is smooth, meanwhile in our study, the crack is represented by a piece of a plane curve whose edges have sharp corners. It has been found in [6] that the difference in stress concentrations near the crack edge between the classical and micropolar case for a crack with a ‘blunt’ tip may be up to 30% in the case when the crack length is equal to 0.26 mm and that for longer cracks, for example 10 mm, this difference is almost negligible. The order of difference is in agreement with our results but in [6] the stress concentration near the crack tip in the micropolar case is lower than in the classical case. At a distance from the crack tip the results obtained in [6] are almost identical to the results of the present investigation. The explanation of discrepancies between our results and those presented in [6] in the vicinity of a crack tip lies in the field of crack geometry. If we consider a smooth contour such as an ellipse of small eccentricity as in [6], then material particles at every point of the contour can rotate under the applied load and generate couple stresses. Consequently, the applied load is distributed between stresses and couple stresses, therefore, the resultant traction is reduced in comparison with the classical case. If we consider a crack whose tip is sharp, then material particles located at the crack tip, i.e. in the corners, can no longer rotate, consequently, couple stresses cannot absorb any part of the applied load but, on the contrary, contribute into the growth of the resultant traction. Hence, the resultant stress in the vicinity of a sharp edge increases. When we move away from the
208
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
corner the particles gain back their ability to rotate and the load is redistributed between the stresses and the couple stresses so the resultant traction starts decreasing. Indirect confirmation of our explanation may be found in [33,34]. In [33] Park and Lakes performed an experimental investigation of torsion of a rectangular micropolar beam. The boundary of any typical crosssection of such beam by a plane perpendicular to the generators contains sharp corners. It has been found that the stress distribution on the boundary of the beam cross-section is significantly higher than in the case when microstructure is ignored. At the same time, the study by Potapenko et al. [34] performed for an elliptic micropolar bar (in this case any typical cross-section of the bar is bounded by a smooth curve) shows that stress concentrations on the boundary of the bar cross-section may be up 15% lower than in the classical case. Similar conclusions may be made when we compare cracks with sharp and ‘blunt’ tips in micropolar medium. 6. Summary In this paper, we have shown that the method introduced by Shmoylova et al. in [23] may be applied to the investigation of stress distribution around a crack with a sharp tip in micropolar medium. We came to the conclusion that material microstructure does have a significant effect on the stress distribution around a crack and demonstrated it using the example of a crack in a human bone. The effect of material microstructure depends on the crack length, crack geometry and has the strongest influence in the vicinity of the crack tip. References [1] A.C. Eringen, Linear theory of micropolar elasticity, J. Math. Mech. 15 (1966) 909–923. [2] W. Nowacki, Theory of Asymmetric Elasticity, Polish Scientific Publishers, Warsaw, 1986. [3] R. Lakes, Experimental methods for study of Cosserat elastic solids and other generalized elastic continua, in: H.B. Muhlhaus (Ed.), Continuum Models for Materials with Microstructure, John Wiley and Sons, New York, 1995, pp. 1–22. [4] R. Lakes, Dynamical study of couple stress effects in human compact bone, J. Biomed. Eng. 104 (1982) 6–11. [5] S. Nakamura, R. Lakes, Finite element analysis of stress concentration around a blunt crack in a Cosserat elastic solid, Comp. Methods Appl. Mech. Eng. 66 (1988) 257–266. [6] R. Lakes, S. Nakamura, J. Behiri, W. Bonfield, Fracture mechanics of bone with short cracks, J. Biomech. 23 (1990) 967–975. [7] V.D. Kupradze et al., Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity, North-Holland, Amsterdam, 1979. [8] D. Iesan, Existence theorems in the theory of micropolar elasticity, Int. J. Eng. Sci. 8 (1970) 777–791. [9] P. Schiavone, Integral equation methods in plane asymmetric elasticity, J. Elasticity 43 (1996) 31–43. [10] S. Potapenko, P. Schiavone, A. Mioduchowski, Antiplane shear deformations in a linear theory of elasticity with microstructure, ZAMP 56 (2005) 516–528. [11] S. Potapenko, P. Schiavone, A. Mioduchowski, On the solution of mixed problems in antiplane micropolar elasticity, Math. Mech. Solids 8 (2003) 151–160. [12] S. Potapenko, A Generalised Fourier approximation in anti-plane micropolar elasticity, J. Elasticity 81 (2005) 159–177. [13] S. Potapenko, Fundamental sequences of functions in the approximation of solutions to mixed boundary-value problems of Cosserat elasticity, Acta Mech. 177 (2005) 61–71. [14] H-B. Mu¨hlhaus, E. Pasternak, Path independent integrals for Cosserat continua and application to crack problems, Int. J. Fract. 113 (2002) 21–26. [15] C. Atkinson, F.G. Leppington, The effect of couple stresses on the tip of a crack, Int. J. Solids Struct. 13 (1977) 1103–1122. [16] R. De Borst, E. van der Giessen, Material Instabilities in Solids, John Wiley, Chichester, 1998. [17] A. Yavari, S. Sarkani, E. Moyer, On fractal cracks in micropolar elastic solids, J. Appl. Mech. 69 (2002) 45–54. [18] E. Diegele, R. Eisasser, C. Tsamakis, Linear micropolar elastic crack-tip fields under mixed mode loading conditions, Int. J. Fract. 4 (2004) 309–339. [19] M. Garajeu, E. Soos, Cosserat models versus crack propagation, Math. Mech. Solids 8 (2003) 189–218. [20] I. Chudinovich, C. Constanda, Variational and Potential Methods in the Theory of Bending of Plates with Transverse Shear Deformation, Chapman&Hall/CRC, Boca Raton, London, New York, Washington, DC, 2000. [21] E. Shmoylova, S. Potapenko, L. Rothenburg, Weak solutions of the interior boundary value problems of plane Cosserat elasticity, ZAMP 57 (2006) 506–522. [22] E. Shmoylova, S. Potapenko, L. Rothenburg, Weak solutions of the exterior boundary value problems of plane Cosserat elasticity, J. Int. Equat. Appl. 19 (2007) 71–92. [23] E. Shmoylova, S. Potapenko, L. Rothenburg, Stress distribution around a crack in plane micropolar elasticity, J. Elasticity 86 (2007) 19–39. [24] C.A. Brebbia, The Boundary Element Method for Engineers, Pentech Press, London, 1978.
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
209
[25] M.H. Aliabadi, C.A. Brebbia, Advances in Boundary Element Methods for Fracture Mechanics, Computational Mechanics Publications, Elsevier Applied Science, Southampton, Boston, London, New York, 1993. [26] I. Sneddon, Crack Problems in the Classical Theory of Elasticity, Wiley, New York, 1969. [27] L. Gaul, M. Ko¨gl, M. Wagner, Boundary Element Methods for Engineers and Scientists, Springer-Verlag, Berlin, Heidelberg, New York, 2003. [28] F. Bouyge, I. Jasiuk, M. Ostoja-Starzewski, A micromechanically based couple-stress model of an elastic two-phase composite, Int. J. Solids Struct. 38 (2001) 1721–1735. [29] I. Jasiuk, M. Ostoja-Starzewski, Modeling of bone at a single lamella level, Biomech. Model. Mechanobiol. 3 (2004) 67–74. [30] F. Bouyge, I. Jasiuk, M. Ostoja-Starzewski, A micropolar model of trabecular bone, in: Proc. Bioeng. Conf., Big Sky, Montana, USA., June 16–20, 1999. [31] D. Miannay, Fracture Mechanics, Springer, New York, Berlin, 1998. [32] M. Janssen, Fracture Mechanics, Spon Press, London, New York, 2004. [33] H. Park, R. Lakes, Cosserat micromechanics of human bone: strain redistribution by a hydration sensitive constituent, J. Biomech. 19 (1986) 385–397. [34] S. Potapenko, P. Schiavone, A. Mioduchowski, Generalized Fourier series solution of torsion of an elliptic beam with microstructure, Appl. Math. Lett. 17 (2004) 189–192.
Boundary element analysis of stress distribution around a crack in plane micropolar elasticity E. Shmoylova, S. Potapenko *, L. Rothenburg Department of Civil Engineering, University of Waterloo, 200 University Avenue West, Waterloo, Ontario, Canada N2L 3G1 Received 20 October 2006; received in revised form 22 November 2006; accepted 30 November 2006 Available online 7 June 2007
Abstract In this paper, we use the boundary element method to find a semi-analytical solution to the problem of stress concentration around a crack in plane micropolar elasticity. We provide an example demonstrating the effect of material microstructure. 2007 Elsevier Ltd. All rights reserved. Keywords: Micropolar elasticity; Weak solutions; Boundary element method; Cracks in micropolar medium
1. Introduction The theory of micropolar elasticity (also known as Cosserat or asymmetric theory of elasticity) was introduced by Eringen [1] (see [2] for a review of works in this area and an extensive bibliography) to eliminate discrepancies between classical theory of elasticity and experiments in cases when effects of material microstructure were known to contribute significantly to the body’s overall deformation, for example, materials with granular microstructure such as polymers or human bones (see [3–6]). These cases are becoming increasingly important in the design and manufacture of modern day advanced materials as small-scale effects become very important in the prediction of the overall mechanical behavior of these materials. For the last 30 years numerous investigations have been performed in this direction. For example, threedimensional problems of Cosserat elasticity have been formulated in a rigorous setting and solved by means of potential theory methods by Kupradze [7]. The corresponding boundary value problems for plane and antiplane shear deformations of a micropolar homogeneous, linearly elastic solid were shown to be well-posed and subsequently solved using the boundary integral equation method by Iesan [8], Schiavone [9], Potapenko et al. [10,11] and Potapenko [12,13].
*
Corresponding author. Tel.: +1 519 888 4567x37156; fax: +1 519 888 6197. E-mail addresses: [email protected] (E. Shmoylova), [email protected] (S. Potapenko), [email protected] (L. Rothenburg). URL: www.civil.uwaterloo.ca/spotapen (E. Shmoylova). 0020-7225/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2007.04.006
200
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Several studies relating to investigations of stress distributions around a crack have been undertaken under assumptions of a simplified theory of plane Cosserat elasticity by Mu¨hlhaus and Pasternak [14], Atkinson and Leppington [15], experimentally by Lakes and Nakamura [5], and using the finite element method (see [6]). There also has been some activity lately in the area of crack analysis in the three-dimensional Cosserat elasticity (see [16–19]). Recently, Chudinovich and Constanda [20] used the boundary integral equation method in a weak (Sobolev) space setting to obtain the solution for fundamental boundary value problems in a theory of bending of classical elastic plates. This approach has wide practical applicability because it also covers domains with reduced boundary smoothness. In addition, it allows to give an answer for the fundamental question relating to the existence and uniqueness of the solution and, second, gives an opportunity to employ an effective numerical procedure which allows the construction of a numerical solution which can be very useful for practical purposes. Further, Chudinovich and Constanda [20] extended their method to accommodate several problems relating to the investigation of stress concentrations around a crack in classical plates. In a number of very recent works by Shmoylova et al. [21–23], the authors performed the rigorous analysis of interior and exterior Dirichlet and Neumann boundary value problems in plane Cosserat elasticity and also obtained the solution to the crack problem arising in this theory in the form of modified integral potentials with unknown distributional densities. Unfortunately, it is very difficult, if not impossible, to find these densities analytically, consequently, we have to find a numerical technique which will allow us to obtain a numerical approximation of the solution. One of the most effective approaches to achieve this goal is to use the boundary element method. This method has been developed by Brebbia [24] and has become very popular among researchers in different areas including fracture mechanics (see, for example, [25] for references on applications of the boundary element method in science and engineering). The boundary element method has originated from works on classical integral equations and finite elements and incorporates advantages of both techniques. On the one hand, it allows us to reduce the dimension of a problem by one and defines domains extending to infinity with a high degree of accuracy exactly as the boundary integral equation method. One the other hand, the boundary element method does not require the differentiation of shape functions, which is the major requirement of the finite element method when we have to find stresses, but allows us to differentiate the matrix of fundamental solutions instead, which makes the calculation of stresses easier and more accurate. In this paper, we use the boundary element method to find the solution for an infinite domain weakened by a crack in plane Cosserat elasticity, when stresses and couple stresses are prescribed along both sides of the crack (Neumann boundary value problem), and discuss its convergence. To illustrate the effectiveness of the method for applications we consider a crack in a human bone which is modelled under assumptions of plane micropolar elasticity. We find the numerical solution for stresses around the crack and show that the solution may be reduced to the classical one if we set all micropolar elastic constants equal to zero. We come to the conclusion that there could be up to 26% difference in quantitative characteristics of the stress around a crack in the micropolar case in comparison with the model when microstructure is ignored (classical case, see for example, [26]).
2. Preliminaries In what follows Greek and Latin indices take the values 1, 2 and 1, 2, 3, respectively, the convention of summation over repeated indices is understood, Mmn is the space of ðm nÞ-matrices, En is the identity element in Mnn , the columns of a (3 · 3)-matrix P are denoted by P ðiÞ , a superscript T indicates matrix transposition, the generic symbol c denotes various strictly positive constants, and ð Þ;a oð Þ=oxa . Let S be a domain in R2 occupied by a homogeneous and isotropic linearly elastic micropolar material with elastic constants k, l, a, c and e. We use the notations k k0;S and h; i0;S for the norm and inner product in L2 ðSÞ \ Mm1 for any m 2 N. When S ¼ R2 , we write k k0 and h; i0 . T The state of plane micropolar strain is characterized by a displacement field uðx0 Þ ¼ ðu1 ðx0 Þ; u2 ðx0 Þ; u3 ðx0 ÞÞ T and a microrotation field /ðx0 Þ ¼ ð/1 ðx0 Þ; /2 ðx0 Þ; /3 ðx0 ÞÞ of the form
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
ua ðx0 Þ ¼ ua ðxÞ; u3 ðx0 Þ ¼ 0; /a ðx0 Þ ¼ 0; /3 ðx0 Þ ¼ /3 ðxÞ;
201
ð1Þ
where x0 ¼ ðx1 ; x2 ; x3 Þ and x ¼ ðx1 ; x2 Þ are generic points in R3 and R2 , respectively. The equilibrium equations of plane micropolar strain written in terms of displacements and microrotations are given by [8,9] Lðox ÞuðxÞ ¼ qðxÞ;
x 2 S;
ð2Þ T
in which now, denoting /3 by u3 , we have uðxÞ ¼ ðu1 ; u2 ; u3 Þ , the matrix partial differential operator Lðox Þ ¼ Lðo=oxa Þ is defined by 0 1 ðk þ l aÞn1 n2 2an2 ðl þ aÞD þ ðk þ l aÞn21 B C LðnÞ ¼ Lðna Þ ¼ @ A; ðk þ l aÞn1 n2 ðl þ aÞD þ ðk þ l aÞn22 2an1 2an2
2an1
ðc þ eÞD 4a
T
where D ¼ na na , and the vector q ¼ ðq1 ; q2 ; q3 Þ represents body forces and body couples. Together with L, we consider the boundary stress operator T ðox Þ ¼ T ðo=oxa Þ defined by 0 1 ðl aÞn1 n2 þ kn2 n1 2an2 ðk þ 2lÞn1 n1 þ ðl þ aÞn2 n2 B C T ðnÞ ¼ T ðna Þ ¼ @ ðl aÞn2 n1 þ kn1 n2 ðk þ 2lÞn2 n2 þ ðl þ aÞn1 n1 2an1 A; 0 0 ðc þ eÞna na T
where n ¼ ðn1 ; n2 Þ is the unit outward normal to oS. To guarantee the ellipticity of system (2), in what follows we assume that k þ l > 0;
l > 0;
c þ e > 0;
a > 0:
The internal energy density is given by 2Eðu; vÞ ¼ 2E0 ðu; vÞ þ lðu1;2 þ u2;1 Þðv1;2 þ v2;1 Þ þ aðu1;2 u2;1 þ 2u3 Þðv1;2 v2;1 þ 2v3 Þ þ ðc þ eÞðu3;1 v3;1 þ u3;2 v3;2 Þ; 2E0 ðu; vÞ ¼ ðk þ 2lÞðu1;1 v1;1 þ u2;2 v2;2 Þ þ kðu1;1 v2;2 þ u2;2 v1;1 Þ: Clearly, Eðu; uÞ is a positive quadratic form. The space of rigid displacements and microrotations F is spanned by columns of the matrix 0 1 1 0 x2 B C F ¼ @ 0 1 x1 A 0 0 1 from which it can be seen that LF ¼ 0 in R2 , T F ¼ 0 on oS and a general rigid displacement can be written as Fk, where k 2 M31 is constant and arbitrary. A Galerkin representation for the solution of (2) when qðxÞ ¼ dðjx yjÞ, where d is the Dirac delta distribution, yields the matrix of fundamental solutions (see [9]) Dðx; yÞ ¼ L ðoxÞtðx; yÞ;
ð3Þ
where L is the adjoint of L a 2 tðx; yÞ ¼ f½k 2 jx yj þ 4 ln jx yj þ 4K 0 ðkjx yjÞg: 4 8pk K 0 is the modified Bessel function of order zero and the constants a; k 2 are defined by a1 ¼ ðc þ eÞðk þ 2lÞðl þ aÞ; In view of (3) and (4) Dðx; yÞ ¼ DT ðx; yÞ ¼ Dðy; xÞ:
k2 ¼
4la : ðc þ eÞðl þ aÞ
ð4Þ
202
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Along with matrix Dðx; yÞ we consider the matrix of singular solutions T
P ðx; yÞ ¼ ðT ðoyÞDðy; xÞÞ : ðiÞ
ð5Þ ðiÞ
2
It is easy to verify that D ðx; yÞ and P ðx; yÞ satisfy (2) with qðxÞ ¼ 0 at all x 2 R ; x 6¼ y. We consider an infinite domain with a crack modelled by an open arc C0 and assume that C0 is a part of a simple closed C 2 -curve C that divides R2 into interior and exterior domains Xþ and X . In what follows we denote by the superscripts + and the limiting values of functions as x ! C from within Xþ or X . Further, we define X ¼ R2 n C0 and C1 ¼ C n C0 . Let H 1;x ðXÞ be the space of all u ¼ fuþ ; u g such that uþ 2 H 1 ðXþ Þ, ± u 2 H 1;x ðX Þ and cþ 1 uþ ¼ c1 u , where c0 are the trace operators on C1 from within X , respectively. The þ space H 1 ðX Þ is a standard Sobolev space and H 1;x ðX Þ is a weighted Sobolev space defined in [21]. Further, we introduce the corresponding single layer and double layer potentials given respectively by Z ðV uÞðxÞ ¼ Dðx; yÞuðyÞdsðyÞ; C Z0 P ðx; yÞuðyÞdsðyÞ; ðW uÞðxÞ ¼ C0
where u 2 M31 is an unknown density matrix. 3. Boundary value problem Let us consider Neumann boundary value problem. We seek u 2 C 2 ðXÞ \ C 1 ðXÞ, u 2 A such that LuðxÞ ¼ 0; þ
x 2 X; þ
ðTuÞ ðxÞ ¼ g ðxÞ;
ðTuÞ ðxÞ ¼ g ðxÞ;
x 2 C0 ;
ðNÞ
where gþ and g are prescribed on C0 . Asymptotic class A was introduced in [9]. The variational formulation of problem (N) is as follows. We seek u 2 H 1;x ðXÞ such that 8v 2 H 1;x ðXÞ; bðu; vÞ ¼ hdg; cþ 0 vþ i0;C0 þ hg ; dvi0;C0 R R where bðu; vÞ ¼ 2 X Eðu ; v Þdx þ 2 Xþ Eðuþ ; vþ Þdx. Note that (6) is solvable only if
hz; dgi0;C0 ¼ 0
8z 2 F:
ð6Þ
ð7Þ
± We introduce the restriction operators p± to X± and the trace operators c 0 on C0 from within X , respectively and define the modified single layer potential V of density u by ðVuÞðxÞ ¼ ðV uÞðxÞ ðV uÞ0 ; ~zðiÞ 0;C ~zðiÞ ðxÞ; x 2 R2 ; 0
where V u is the single layer potential, and V0 is the boundary operator defined by ðV uÞ0 ¼ c 0 p V u, and 3 f~zðiÞ gi¼1 is an L2 ðC0 Þ-orthonormal basis for F. Also we introduce the modified double layer potential W of density w ðWwÞðxÞ ¼ ðW wÞðxÞ p0 W þ w; ~zðiÞ 0;C ~zðiÞ ðxÞ; x 2 X; 0
where p0 is the restriction operator to C0. Solution of problem (6) may be represented in the form u ¼ ðVuÞX þ Ww þ z;
ð8Þ
where u and w are unknown densities, and z 2 F is arbitrary. The detailed procedure of obtaining solution (8) has been developed by Shmoylova et al. in [23]. 4. Boundary element method Consider problem (6). As shown in [23], the solution to this problem may be represented in the form (8) and the corresponding boundary integral equations are uniquely solvable with respect to distributional densities u
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
203
and w. As we stated above, these densities cannot be found analytically. To approximate them numerically we use the boundary element method [27], which makes use of the following classical result. Lemma 1 (Somigliana formula). Using classical techniques, it can be proved that if u 2 H 1;x ðXÞ is a solution of Lu = 0 in X, then Z
1 ½Dðx; yÞdðT ðoy ÞuðyÞÞ P ðx; yÞduðyÞdsðyÞ ¼ duðxÞ; 2 C0
x 2 C0 ;
ð9Þ
where dðT ðoy ÞuðyÞÞ denotes the jump of T ðoy ÞuðyÞ on the crack. It has been shown in [23], that the density of the modified single layer potential may be found in the form u ¼ dðT ðoy ÞuðyÞÞ ¼ dg. Now we need to find the density of the modified double layer potential w ¼ du. To ðkÞ achieve this goal we divide C0 into n elements C0 , each of which possesses one node nðkÞ located in the middle of the element. The values of dg and du are constant throughout the element and correspond to the value at the node dgðnðkÞ Þ and duðnðkÞ Þ. Then (9) becomes n Z X 1 ½Dðx; yÞdgðnðkÞ Þ P ðx; yÞduðnðkÞ ÞdsðyÞ ¼ duðxÞ; x 2 C0 : ðkÞ 2 C0 k¼1 Placing x sequentially at all nodes, we obtain the linear algebraic system of equations ! ! Z Z n n X X 1 ðiÞ ðkÞ ðiÞ Dðn ; yÞdsðyÞ dgðn Þ P ðn ; yÞdsðyÞ duðnðkÞ Þ ¼ duðnðiÞ Þ; i; k ¼ 1; n ðkÞ ðkÞ 2 C0 C0 k¼1 k¼1
ð10Þ
with respect to duðn R i Þ. We note that CðkÞ DðnðiÞ ; yÞdsðyÞ are defined for any i and k, as in [9]. 0 Solving (10) we construct the approximation to w. If we introduce the shape function Uk ðxÞ by ( ðkÞ 1; x 2 C0 ; Uk ðxÞ ¼ ðkÞ 0; x 2 C0 n C0 ; Pn Pn then the approximated densities u and w are uðnÞ ðxÞ ¼ k¼1 Uk ðxÞdgðnðkÞ Þ and wðnÞ ðxÞ ¼ k¼1 Uk ðxÞduðnðkÞ Þ and the approximate solution is becoming of the form uðnÞ ¼ ðVuðnÞ ÞX þ WwðnÞ þ z, where z is arbitrary. We now have to prove that the approximate solution uðnÞ will converge to the exact analytical solution u when n ! 1. Let us formulate the following theorem. Theorem 2. uðnÞ ! u as n ! 1. Proof. Since we consider the Neumann problem, rigid displacement terms are not determined. Consequently, it is enough to show that V uðnÞ ! V u and W wðnÞ ! W w as n ! 1. Consider V uðnÞ . For x 2 X ! Z 3 Z n X X ðkÞ ðnÞ ðiÞ ðiÞ jV uðxÞ V u ðxÞj 6 D ðx; yÞuðyÞdsðyÞ D ðx; yÞdsðyÞ uðn Þ ðkÞ C0 C0 i¼1 k¼1 3 X n Z X ¼ ½DðiÞ ðx; yÞuðyÞ DðiÞ ðx; yÞuðnðkÞ ÞdsðyÞ ðkÞ C0 i¼1 k¼1 3 X n Z X ¼ DðiÞ ðx; yÞ½uðyÞ uðnðkÞ ÞdsðyÞ ðkÞ C i¼1 k¼1 0
6
3 X
n X
i¼1
k¼1
kDðiÞ ðx; Þk1;x;X juðyÞ uðnðkÞ Þjhk ;
where hk is the length of the kth element assuming that all hk ¼ h ¼ Ln, L is the length of C0 , and P3 P3 P ðkÞ ðkÞ juðyÞ uðn Þj 6 i¼1 jui ðyÞ ui ðn Þj ¼ i¼1 2a¼1 joa ui ðnðkÞ Þjh þ Oðh2 Þ.
204
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Denote by M 1 ¼ maxa¼1;2;i¼1;3;k¼1;n joa ui ðnðkÞ Þj. Since kDðiÞ ðx; Þk1;x;X is uniformly bounded [8,9], i.e. there exists M 2 > 0 such that kDðiÞ ðx; Þk1;x;X 6 M 2 for any x 2 X, we obtain that jV uðxÞ V uðnÞ ðxÞj 6
18M 1 M 2 L2 ! 0 as n ! 1: n
Repeating the proof for W wðnÞ , we conclude that uðnÞ ! u as n ! 1. h 5. Example As an example, we consider a longitudinal crack inside a human bone in the case when constant normal stretching pressure of magnitude p is applied on both sides of the crack. If we consider a typical transversal cross-section of the bone and assume that this cross-section is small enough then the deformation of each cross-section under the prescribed load will be the same throughout the length of the bone and will develop in the plane of the cross-section. Consequently, such deformations may be considered under assumptions of plane micropolar elasticity. Such a model is not an idealization that lies far from reality, as it may seem at first, but, as shown, for example, in [28–30], can describe actual cracks in bones very closely, since orthopedic biomechanics usually deals with cracks of a very small size. We model a crack as an open arc of the circle given by equations: x1 ¼ a cos h, x2 ¼ a sin h, h 2 ð0; p=6Þ (Fig. 1). By changing the radius a of the circle we will change the length of the crack. We are interested in how the normal traction distributes at a distance from the crack tip along the line: x1 ¼ a, x2 < 0. Clearly, this problem can be considered as the Neumann problem described above. Elastic constants for a human bone have been measured in [5] and take the following values: a = 4000 MPa, c = 193.6 N, e = 3047 N, k = 5332 GPa, l = 4000 MPa. In our example we construct solutions for cracks of lengths equal to 0.26 mm, 0.52 mm, 0.75 mm, and 10 mm to show good agreement of our results with those presented in the experimental study by Nakamura and Lakes [5] performed on human bone cracks of same lengths. We also assume that the normal stretching pressure p takes a value of 2 MPa. Let the distance from the tip of the crack be q ¼ jx2 j. The numerical solution for boundary tractions and moments has been found to approximate the exact solution to five decimal places for n = 52 elements of C0 (see Table 1). Let us now compare the results for the normal traction in the micropolar case with the results of the classical theory. The classical case may be obtained from the solution for micropolar elasticity setting the micropolar elastic constants equal to zero. In Figs. 2–5 there is a graphical representation for the distribution of the normal traction at a distance from the lower crack tip for crack lengths equal to 0.26 mm, 0.52 mm, 0.75 mm
Fig. 1. Crack in x1 x2 -plane. Length: 0.5 mm – solid line (a = 1 mm).
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
205
Table 1 Approximate solution for a crack 0.1
0.5
0.7
1
n=4
Distance q (mm) Tn (MPa) Ts (MPa) M 3 (N/m)
0.983256 1.175489 162.5682
0.108643 0.538926 91.24553
0.082351 0.213505 79.45362
0.057282 0.756437 70.64523
n = 10
Tn (MPa) Ts (MPa) M 3 (N/m)
0.569361 0.634936 84.24634
0.068735 0.264282 49.86301
0.052678 0.091475 44.26856
0.032536 0.045343 36.09357
n = 30
Tn (MPa) Ts (MPa) M 3 (N/m)
0.427549 0.506874 69.24764
0.061862 0.184756 34.25447
0.039754 0.075982 30.25879
0.021830 0.033547 24.62579
n = 50
Tn (MPa) Ts (MPa) M 3 (N/m)
0.398462 0.454906 61.11549
0.053918 0.108063 29.69127
0.030028 0.061039 22.59611
0.015244 0.026688 15.80458
n = 52
Tn (MPa) Ts (MPa) M 3 (N/m)
0.398456 0.454899 61.11542
0.053910 0.108054 29.69119
0.030021 0.061033 22.59604
0.015239 0.026678 15.80451
Length: 0.52 mm. Tn – normal traction, Ts – tangential traction, M 3 – moment about x3-axis.
Fig. 2. Normal traction on the edge of the crack. Length: 0.26 mm, (a) micropolar, (b) classical.
and 10 mm correspondingly. The traction is divided by the applied load p to represent the data in non-dimensional values. The bold curve characterizes the stress distribution in the micropolar case while the classical case is plotted by the ordinary curve. The distance between the first point, in which we compute the normal traction, and the tip of the crack is equal to one fifth of the length of the crack. It may be seen that the normal traction is significantly higher in the vicinity of the crack tip in the micropolar case in comparison with the case when microstructure is ignored (classical theory), particularly, for cracks of 0.26 mm, 0.52 mm and 0.75 mm long. In the case when the length of the crack is equal to 0.75 mm we can observe that the normal traction in the vicinity of the crack tip is 26.8% higher under
206
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
Fig. 3. Normal traction on the edge of the crack. Length: 0.52, (a) micropolar, (b) classical.
Fig. 4. Normal traction on the edge of the crack. Length: 0.75, (a) micropolar, (b) classical.
assumptions of Cosserat elasticity in comparison with the classical case. At the same time, for the crack of the linear size equal to 10 mm, the traction in the vicinity of the crack tip differs only by 13%. The fact that this difference is still present and significantly lower than in the case of the shorter cracks may be explained by the size effect of the crack opening, which is of the order of 105 for this crack in comparison with that of 106 in the other cases. When it comes to the consideration of stresses at a distance from the crack tip, we can conclude that the traction in the micropolar case decays faster than in the classical case and approximately at a distance of one
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
207
Fig. 5. Normal traction on the edge of the crack. Length: 10 mm, (a) micropolar, (b) classical.
crack length the values of the normal traction in both cases become equal to each other. Farther from the crack tip, the traction in the micropolar case becomes lower than in the classical case, especially when we consider the crack of 0.26 mm long for which, as it may be seen from Fig. 2, the difference is drastic and may be up to 19.8%. Further, it may be observed that at a distance of approximately three crack lengths from the tip of the crack the effect of the crack on stresses is negligible in agreement with Saint-Venant’s principle. The results presented in Figs. 2–5 may be compared with earlier investigations undertaken in the area of Cosserat solids. Lakes and Nakamura [5] performed an experimental study on stress concentrations around cracks in a human bone. They considered the same crack lengths as in the present article. Later, the same authors validated their results presented in [5] using finite element method [6]. However, it should be noted, that direct comparison with the results presented in [5,6] does not seem feasible since our formulation of the crack problem is different from that adopted in [5,6], where a crack is considered as a ‘blunt’ notch similarly to the approach of classical fracture mechanics [26,31,32], in which a crack is usually modelled as a ‘squashed’ ellipse of small eccentricity. The tip of the crack under consideration in [5,6], therefore, is smooth, meanwhile in our study, the crack is represented by a piece of a plane curve whose edges have sharp corners. It has been found in [6] that the difference in stress concentrations near the crack edge between the classical and micropolar case for a crack with a ‘blunt’ tip may be up to 30% in the case when the crack length is equal to 0.26 mm and that for longer cracks, for example 10 mm, this difference is almost negligible. The order of difference is in agreement with our results but in [6] the stress concentration near the crack tip in the micropolar case is lower than in the classical case. At a distance from the crack tip the results obtained in [6] are almost identical to the results of the present investigation. The explanation of discrepancies between our results and those presented in [6] in the vicinity of a crack tip lies in the field of crack geometry. If we consider a smooth contour such as an ellipse of small eccentricity as in [6], then material particles at every point of the contour can rotate under the applied load and generate couple stresses. Consequently, the applied load is distributed between stresses and couple stresses, therefore, the resultant traction is reduced in comparison with the classical case. If we consider a crack whose tip is sharp, then material particles located at the crack tip, i.e. in the corners, can no longer rotate, consequently, couple stresses cannot absorb any part of the applied load but, on the contrary, contribute into the growth of the resultant traction. Hence, the resultant stress in the vicinity of a sharp edge increases. When we move away from the
208
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
corner the particles gain back their ability to rotate and the load is redistributed between the stresses and the couple stresses so the resultant traction starts decreasing. Indirect confirmation of our explanation may be found in [33,34]. In [33] Park and Lakes performed an experimental investigation of torsion of a rectangular micropolar beam. The boundary of any typical crosssection of such beam by a plane perpendicular to the generators contains sharp corners. It has been found that the stress distribution on the boundary of the beam cross-section is significantly higher than in the case when microstructure is ignored. At the same time, the study by Potapenko et al. [34] performed for an elliptic micropolar bar (in this case any typical cross-section of the bar is bounded by a smooth curve) shows that stress concentrations on the boundary of the bar cross-section may be up 15% lower than in the classical case. Similar conclusions may be made when we compare cracks with sharp and ‘blunt’ tips in micropolar medium. 6. Summary In this paper, we have shown that the method introduced by Shmoylova et al. in [23] may be applied to the investigation of stress distribution around a crack with a sharp tip in micropolar medium. We came to the conclusion that material microstructure does have a significant effect on the stress distribution around a crack and demonstrated it using the example of a crack in a human bone. The effect of material microstructure depends on the crack length, crack geometry and has the strongest influence in the vicinity of the crack tip. References [1] A.C. Eringen, Linear theory of micropolar elasticity, J. Math. Mech. 15 (1966) 909–923. [2] W. Nowacki, Theory of Asymmetric Elasticity, Polish Scientific Publishers, Warsaw, 1986. [3] R. Lakes, Experimental methods for study of Cosserat elastic solids and other generalized elastic continua, in: H.B. Muhlhaus (Ed.), Continuum Models for Materials with Microstructure, John Wiley and Sons, New York, 1995, pp. 1–22. [4] R. Lakes, Dynamical study of couple stress effects in human compact bone, J. Biomed. Eng. 104 (1982) 6–11. [5] S. Nakamura, R. Lakes, Finite element analysis of stress concentration around a blunt crack in a Cosserat elastic solid, Comp. Methods Appl. Mech. Eng. 66 (1988) 257–266. [6] R. Lakes, S. Nakamura, J. Behiri, W. Bonfield, Fracture mechanics of bone with short cracks, J. Biomech. 23 (1990) 967–975. [7] V.D. Kupradze et al., Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity, North-Holland, Amsterdam, 1979. [8] D. Iesan, Existence theorems in the theory of micropolar elasticity, Int. J. Eng. Sci. 8 (1970) 777–791. [9] P. Schiavone, Integral equation methods in plane asymmetric elasticity, J. Elasticity 43 (1996) 31–43. [10] S. Potapenko, P. Schiavone, A. Mioduchowski, Antiplane shear deformations in a linear theory of elasticity with microstructure, ZAMP 56 (2005) 516–528. [11] S. Potapenko, P. Schiavone, A. Mioduchowski, On the solution of mixed problems in antiplane micropolar elasticity, Math. Mech. Solids 8 (2003) 151–160. [12] S. Potapenko, A Generalised Fourier approximation in anti-plane micropolar elasticity, J. Elasticity 81 (2005) 159–177. [13] S. Potapenko, Fundamental sequences of functions in the approximation of solutions to mixed boundary-value problems of Cosserat elasticity, Acta Mech. 177 (2005) 61–71. [14] H-B. Mu¨hlhaus, E. Pasternak, Path independent integrals for Cosserat continua and application to crack problems, Int. J. Fract. 113 (2002) 21–26. [15] C. Atkinson, F.G. Leppington, The effect of couple stresses on the tip of a crack, Int. J. Solids Struct. 13 (1977) 1103–1122. [16] R. De Borst, E. van der Giessen, Material Instabilities in Solids, John Wiley, Chichester, 1998. [17] A. Yavari, S. Sarkani, E. Moyer, On fractal cracks in micropolar elastic solids, J. Appl. Mech. 69 (2002) 45–54. [18] E. Diegele, R. Eisasser, C. Tsamakis, Linear micropolar elastic crack-tip fields under mixed mode loading conditions, Int. J. Fract. 4 (2004) 309–339. [19] M. Garajeu, E. Soos, Cosserat models versus crack propagation, Math. Mech. Solids 8 (2003) 189–218. [20] I. Chudinovich, C. Constanda, Variational and Potential Methods in the Theory of Bending of Plates with Transverse Shear Deformation, Chapman&Hall/CRC, Boca Raton, London, New York, Washington, DC, 2000. [21] E. Shmoylova, S. Potapenko, L. Rothenburg, Weak solutions of the interior boundary value problems of plane Cosserat elasticity, ZAMP 57 (2006) 506–522. [22] E. Shmoylova, S. Potapenko, L. Rothenburg, Weak solutions of the exterior boundary value problems of plane Cosserat elasticity, J. Int. Equat. Appl. 19 (2007) 71–92. [23] E. Shmoylova, S. Potapenko, L. Rothenburg, Stress distribution around a crack in plane micropolar elasticity, J. Elasticity 86 (2007) 19–39. [24] C.A. Brebbia, The Boundary Element Method for Engineers, Pentech Press, London, 1978.
E. Shmoylova et al. / International Journal of Engineering Science 45 (2007) 199–209
209
[25] M.H. Aliabadi, C.A. Brebbia, Advances in Boundary Element Methods for Fracture Mechanics, Computational Mechanics Publications, Elsevier Applied Science, Southampton, Boston, London, New York, 1993. [26] I. Sneddon, Crack Problems in the Classical Theory of Elasticity, Wiley, New York, 1969. [27] L. Gaul, M. Ko¨gl, M. Wagner, Boundary Element Methods for Engineers and Scientists, Springer-Verlag, Berlin, Heidelberg, New York, 2003. [28] F. Bouyge, I. Jasiuk, M. Ostoja-Starzewski, A micromechanically based couple-stress model of an elastic two-phase composite, Int. J. Solids Struct. 38 (2001) 1721–1735. [29] I. Jasiuk, M. Ostoja-Starzewski, Modeling of bone at a single lamella level, Biomech. Model. Mechanobiol. 3 (2004) 67–74. [30] F. Bouyge, I. Jasiuk, M. Ostoja-Starzewski, A micropolar model of trabecular bone, in: Proc. Bioeng. Conf., Big Sky, Montana, USA., June 16–20, 1999. [31] D. Miannay, Fracture Mechanics, Springer, New York, Berlin, 1998. [32] M. Janssen, Fracture Mechanics, Spon Press, London, New York, 2004. [33] H. Park, R. Lakes, Cosserat micromechanics of human bone: strain redistribution by a hydration sensitive constituent, J. Biomech. 19 (1986) 385–397. [34] S. Potapenko, P. Schiavone, A. Mioduchowski, Generalized Fourier series solution of torsion of an elliptic beam with microstructure, Appl. Math. Lett. 17 (2004) 189–192.