HUANG Cheng，YIN Maoshu，QI Xiao，GUO Hun

1.School of Aviation and Mechanical Engineering，Changzhou Institute of Technology，Changzhou 213032，P.R.China；2.Physical Power Division，Shanghai Institute of Space Power?Sources，Shanghai 200240，P.R.China

Abstract: This paper proposes a straightforward and concise approach to analyze the Saint-Venant’s torsion of a circular shaft containing multiple elliptical inclusions or cracks based on the complex variable method. The complex potentials are first derived for the shaft with N elliptical inclusions by introducing Faber series expansion，and then the shear stresses and torsional rigidity are calculated. When the inclusions degenerate into cracks，the solutions for the intensity factors of stress are obtained. Finally，several numerical examples are carried out to discuss the effects of geometry parameters，different shear modulus ratios and array-types of the elliptical inclusions/cracks on the fields of stresses. The obtained results show that the proposed approach has advantages such as high accuracy and good convergence.

Key words：Saint-Venant’s torsion；complex variable method；Faber series

0 Introduction

Circular shafts under torsion are widely used in engineering. To raise the undergoing-load level of the shafts，they are often holed or made into com?posite shafts reinforced by other materials. Thus，it is of not only theoretical interest but also practical importance to study the torsion of circular shafts containing holes or inclusions.

Ling［1］investigated the Saint-Venant’s torsion problem of a circular bar with a ring of uniformly dis?tributed circular holes of equal radii using a special class of harmonic functions introduced by How?land［2］. Kuo et al.［3-4］studied the torsion of a circular tube with circular holes and a cylinder which is rein?forced with circular inclusions by constructing a real stress function，respectively. Jaswon et al.［5］pro?posed an integral equation solution for the classical torsion problem of Saint-Venant through numerical?ly solving a Neumann-type boundary-value equation.

The solution for torsion problem of a circular shaft can always be solved by using numerical tech?niques［6-11］. Katsikadelis et al.［7］presented the bound?ary element solution for the Saint-Venant torsion problem of composite cylindrical bars of arbitrary cross section. Spountzakis et al.［8］developed the boundary element method for the nonuniform tor?sion of composite bars of arbitrary constant cross section by using domain discretization and an effec?tive Gaussian integration over domains of arbitrary shape. Li et al.［9］studied the Saint-Venant’s torsion problem of the arbitrarily shaped bar made of differ?ent materials based on finite element method. Refs.［10-11］adopted the null-field approach to solve the Saint-Venant’s problem of a circular bar with circu?lar holes or inclusions，respectively.

Recently，the research on torsion of advanced materials has attracted extensive attention. Ecsedi et al.［12］generalized the known elastic solution of Saint-Venant’s torsional problem developed by Prandtl to piezoelectric beams，and later they inves?tigated the Saint-Venant torsion of non-homoge?neous and circular cylinder made of orthotropic piezoelectric material［13］. Wang et al.［14］studied the effects of surface elasticity in the Saint-Venant tor?sion problem. Hassani et al.［15］analyzed the Saint-Venant torsion of an orthotropic bar with multiple curved cracks. However，it should be noted that the above work was all made by real variable method.

Muskhelishvili［16］developed a complex variable method to address the Saint-Venant’s torsion of composite circular shafts and solved the problem of a circular shaft containing an eccentric circular inclu?sion. Based on the Muskhelishvili’s method，Yue et al.［17］dealt with the torsion problem of a compos?ite cylinder with cracks and inclusions by introduc?ing the Mellin transforms and solving a set of mixedtype integral equations. Refs.［18-19］showed that，based on the complex variable method，the interac?tion between porous/inclusions can be effectively solved. However，to the author’s knowledge，no work can be found for the solution to the Saint-Ve?nant’s torsion of a circular shaft containing multiple elliptical inclusions based on the complex variable method.

In this paper，we propose a straightforward and concise approach to analyze the problem of interact?ing elliptical inclusions in a circular shaft of torsion based on complex variable theory. The key step in the present work is to express the complex poten?tials in the matrix with elliptical holes（a multiplyconnected region）in the form of Faber series，and then the continuous conditions between the inclu?sions and the matrix are used to determine the un?known coefficients involved in these complex poten?tials. Thus，the novel feature of this paper is to pres?ent a straightforward and concise method to solve the problem of the Saint-Venant’s torsion of a circu?lar shaft containing multiple elliptical inclusions or cracks effectively with high accuracy.

1 Basic Equations

In a rectangular coordinate system x?y?z，con?sider a circular shaft containing N elliptical inclu?sions which are parallel to each other along the z di?rection. The cross-section of the shaft is shown in Fig.1，where apand bp( p=1，2，…，N ) are the lengths of the elliptical inclusions’semi-axis，and zp0are the center coordinates of inclusions，respec?tively. All the inclusions are assumed to be com?pletely bounded to the matrix. The boundaries of the inclusions and the outer contour of the shaft are denoted by Lk(k=0，1，2，…，N). We now study Saint-Venant torsion problem of the composite shaft loaded by the torque T applied at its two ends.

Fig.1 Torsion of a circular shaft containing multiple ellipti?cal inclusions

In this case，the components of displacement（u，v，w）can be expressed as

where τ is the angle of twist per unit length along the z direction and φ(x，y) the warping function.The corresponding stresses are

where μ0is the shear modulus of the matrix，and μk(k=1，2，…，N) are the shear modulus of the k-th inclusion. For this problem，the equilibrium equa?tion becomes that［16］

where the body force is neglected. Substituting Eq.（2）into Eq.（3）leads to the Laplace equation

The solution of Eq.（4）is

where F(z) is called as the complex potential. Sub?stituting Eq.（5）into Eqs.（1—2），the displacement and stress can finally be expressed as［16］

Once F(z) is obtained，the torsional rigidity D can be calculated by［16］

To solve F(z)，we introduce the resultant trac?tion p on any boundary as

Inserting Eq.（7）into Eq.（9）leads to

where C0is a constant and it can be assumed to be zero without affecting the stresses. Since there is no external traction on the cylindrical surface，p=0 and the boundary condition on L0can be derived from Eq.（10）that

where the subscripts“0”denote the matrix. At the interface between the matrix and inclusions Lk(k=0，1，2，…，N)，the continuity conditions require that

where wkand pkstand for the displacement and re?sultant traction along the kth inclusion’s boundary，respectively. Substituting Eqs.（6，10）into Eqs.（12，13），we have

2 Theoretical Analyses

In this case，the matrix is a multiple-connected region containing N elliptic holes and enclosed by the circle，so the complex potential in the matrix has the form that

where f0(z) is an analytical function inside L0，and fk(z) is another analytical function outside the ellipti?cal hole Lk. Thus f0(z) can be expanded into the Taylor series as

where Ajare unknown coefficients. Introduce the following conformal mapping function

which conformably maps the region outside the ellip?tical hole Lkonto the external region of a unit circle of ξk=eiθin the ξk-plane，and thus fk(z) can be ex?panded into the Laurent series as

On the surface of the shaft L0，z=R0eiθ=R0σ and thus Eq.（21）has the form that

which can be used to determine the coefficients by Eq.（24）. Substituting Eq.（22）into Eq.（11），after some rearrangement，we can get that

On the other hand，for any inclusion Lp( p=1，2，…，N )，moving the origin of the global system x-y in to the point zp，that is，making the follow?ing coordinate translation：z-zp=z*，one can ex?press the complex potential in the matrix and inside the inclusions， in the local coordinate system xp-yp，as

From Eqs.（29—31），we can obtain the coeffi?cientsIn Eq.（27），the termis a given function that is analytic inside the inclusion Lp( p ≠k)，and thus it can be expanded into the Faber series as［20-21］

Inserting Eqs.（27，28）into Eqs.（14，15）yields that

where m=1，2，…，M and

In detail，a system of linear equations with re?spect to the unknown coefficientscan be obtained by equating the corresponding coefficients on the two sides of Eqs.（40—45）. In addition，Cm(m ≠0) can be determined by equaling the con?stant terms on the two sides of Eq.（39），but they have no influence on stresses，and thus are ignored.

Once all the complex coefficients are deter?mined from these linear equations，the stresses both in the matrix and the inclusions can be obtained by［16］

At the same time，the torsional rigidity can be easily determined by

where

If one assumes μk=0，and the short-axis radi?us of the ellipse bk=0，thus the inclusion Lkdegen?erates into a crack parallel to the x-axis. In this case one can calculate the stress intensity factor for the shaft with the crack by［16］

where ξ0is the point at the unit circle of ξ plane，and it is corresponding to the tip of the crack in the zplane. Substituting Eq.（18）into Eq.（51），we final?ly have

3 Numerical Results and Discussion

3.1 Torsion of a circular shaft containing a cir?cular inclusion

In order to show validity of the present meth?od，we make some comparisons with the related work in Ref.［11］，which considered a circular shaft of radius containing a circular inclusion，as shown in Fig.2. In the example we take the ratio of R1/R2=0.3，l/R0=0.6 and the shear modulus ratio g1=μ1/μ0=0.6. The numerical results for the non-di?mensional torsional rigidity D*=2D/() ver?sus the number of power series terms are shown in Fig. 3. It is shown that the results are accurate enough when the numbers of power series terms are above five. Furthermore，the variations of D*as functions of the shear modulus ratio g1are shown in Table 1，where results in Refs.［16，17］using the in?tegral formulation and in Ref.［11］based on the nullfield integral approach are also listed. It can be seen that the results obtained in the present work are well consistent with those in the previous works.

Fig.2 Torsion of a circular shaft containing a circular inclu?sion

Fig.3 Torsional rigidity versus the number of power series terms

Table 1 Torsional rigidity of a circular shaft containing a circular inclusion

3.2 Torsion of a circular shaft containing mul?tiple circular inclusions

Consider the case of four circular inclusions lo?cated in a rhombic array，as shown in Fig.4. For comparison with previous work，we take the same parameters as those in Refs.［3，4，11］as follows：l/R0=0.6， gk=29.4 and Rk/R0=0.25 (k=1，2，3，4). The results of the non-dimensional tor?sional rigidity D*are shown in Table 2，and it is shown that the present solutions are in very agree?ment with the work in Ref.［11］，and more accurate than the results in Refs.［3，4］.

Fig.4 Torsion of a circular shaft containing four circular in?clusions of equal radii

Table 2 Torsional rigidity of a circular shaft containing four circular inclusions

As special cases，the torsion of a circular shaft containing three circular holes of equal radii is also given，as shown in Fig.5，where l/R0=0.6，gk=0 and Rk/R0=0.2(k=1，2，3). When taking the series expansion as M =20，the shear stresses along the boundaries are shown in Table 3，and it can be found that the present solutions are well con?sistent with the results in Ref.［1］.

Fig.5 Torsion of a circular shaft containing three circular holes of equal radii

Table 3 Shear stresses σθ/μτR0 along the boundaries (M=20)

3.3 Torsion of a circular shaft containing mul?tiple elliptical inclusions

Fig.6 Torsion of a circular shaft containing two elliptical in?clusions

Fig.7 Non-dimensional torsional rigidity D* versus b1/a1

Consider the case of the circular shaft with two elliptical inclusions，as shown in Fig.6，where a1=a2=0.2R0，b1=b2as variables，l/R0=0.6，and μk=μ(k=1，2). When g=μ/μ0＞1，the inclu?sions are called as hard inclusions，and when g ＜1，they are called as soft inclusions. Specially，when g=0，the inclusion becomes an elliptic hole. The non-dimensional torsional rigidity D*for different values of b1/a1is shown in Fig.7，which shows that for the soft inclusions，the torsional rigidity greatly decreases as the inclusions become softer，and for the hard inclusions，the torsional rigidity greatly in?creases as the inclusions become harder. On the oth?er hand，the stresses at the interface between the matrix and the inclusions are shown in Figs.8—10 for three different cases. It is found from Fig.8 that when g=0，the greatest stress occurs at the points on the hole’s boundary nearest to the exterior sur?face of the shaft. Fig.9 is the case of the shaft with two soft elliptic inclusions (g=0.5)，and it is found that the maximum stress σθalong the surface be?tween the matrix and inclusions decreases compared with that in the case of g=0. However，for the case of two hard inclusions，the maximum value of the stress occurs on the exterior surface of the shaft rather than at the boundary of the inclusions，as shown in Fig.10.

Fig.8 Shaft with two elliptic holes:Stress along the bound?aries of the shaft and holes

Fig.9 Shaft with two soft elliptic inclusions:Stress along the boundaries of the shaft and inclusions

Fig.10 Shaft with two hard elliptic inclusions:Stress along the boundaries of the shaft and inclusions

3.4 Torsion of a circular shaft containing mul?tiple elliptical inclusions and cracks

Consider the case of the circular shaft with two elliptical inclusions and two cracks shown in Fig.11， where a1=a2as variables， a3=a4=0.1R0，b1=b2=0，b3=b4=0.2R0，l/R0=0.6，μ1=μ2=0，μ3=μ4=μ.The non-dimensional tor?sional rigidity D*for different values of a1/R0are shown in Fig.12，which shows that with the larger ratio of a1/R0，the loss of relative D*changes great?er. Owing to its geometrical symmetry，we just dis?cuss the stress intensity factors of the crack on the right of the circular. Therefore，the non-dimensional stress intensity factorsfor different values ofare shown in Fig.13. It is shown that with the increase of the g，the variation of the stress intensity factoris larger than that of，and as the crack length becomes large，thesharply increases while thedecreases slowly.

Fig.11 Torsion of a circular shaft containing two elliptical inclusions and cracks

Fig.12 Non-dimensional torsional rigidity D* versus a1/R0

Fig.13 Non-dimensional intensity factor andversus a1/R0

Fig.14 Torsion of a circular shaft containing a center circu?lar hole and a crack

As a special case，the solutions for the torsion of a circular shaft containing a center circular hole and a crack，as shown in Fig.14，are also given. In the example，we take R1/R0=0.2，and l/R0=0.5，0.6，0.7，respectively. The results of D*for dif?ferent values of a2/R0are given in Fig.15. It is found that with the larger ratio of l/R0，the loss of relative D*changes greater. The non-dimensional stress intensity factorsandfor differ?ent values of a2/R0are shown in Fig.16，which shows that when the crack length approaches the edge of the center hole and boundary of the shaft，andsharply increase due to the inter?action.

Fig.15 Non-dimensional torsional rigidity D* versus a2/R0

Fig.16 Non-dimensional intensity factor and versus a2/R0

4 Conclusions

We studied the Saint-Venant’s torsion of a cir?cular shaft with multiple elliptical inclusions with dif?ferent material constants from the matrix. Based on the complex variable method，the complex poten?tials are expressed in the form of Faber series，and then their unknown coefficients are solved by the continuous conditions at the interface between the matrix and inclusions. Solutions for the cases of mul?tiple circular inclusions/elliptic holes/cracks can be easily obtained as special cases of the present work，and they are compared with previous results ob?tained based on the integral equation method，the null-field integral approach or other numerical meth?ods. It is shown that the present work has advantag?es such as high accuracy and good convergence. Fur?thermore，several numerical examples for the inter?action between multiple elliptical inclusions/holes/cracks are presented to discuss the effects of the pa?rameters of these defects on the stress，torsional ri?gidity and the stress intensity factors，and it is found that the geometry size，material constants and loca?tions of the defects play a significant role in these variables of fields.