Y.C.Shiah,C.L.Tanand Li-Ding Chan

Boundary Element Analysis of Thin Anisotropic Structures by a Self-regularization Scheme

Y.C.Shiah1,C.L.Tan2,3and Li-Ding Chan1

In the conventional boundary element method(BEM),the presence of singular kernels in the boundary integral equation or integral identities causes serious inaccuracy of the numerical solutions when the source and field points are very close to each other.This situation occurs commonly in elastostatic analysis of thin structures.The numerical inaccuracy issue can be resolved by some regularization process.Very recently,the self-regularization scheme originally proposed by Cruse and Richardson(1996)for 2D stress analysis has been extended and modified by He and Tan(2013)to 3D elastostatics analysis of isotropic bodies.This paper deals with the extension of the technique developed by the latter authors to the elastostatics analysis of 3D thin,anisotropic structures using the self-regularized displacement boundary integral equation(BIE).The kernels of the BIE employ the double Fourier-series representations of the fundamental solutions as proposed by Shiah,Tan and Wang(2012)and Tan,Shiah and Wang(2013)recently.Numerical examples are presented to demonstrate the veracity of the scheme for BEM analysis of thin anisotropic bodies.

Self-regularization,boundary element method,nearly singular integrals,thin anisotropic bodies.

1 Introduction

Thin-walled structures are commonly seen in engineering applications,e.g.,pressure vessels and surface panels of flight vehicles,and composite materials are increasingly being used for their construction.Although the boundary element method(BEM)is well recognized as a powerful numerical tool for engineering analysis,itsusefortreatingthin-walledstructuresrequiresspecialattentionbecause of the presence of singular kernel functions in the boundary integral equation(BIE)which gives rise to nearly singular integrals.These nearly singular integrals occur because the source point on one surface becomes very close to the field points on the elements which lie on the the other side of the thin wall.Because of the rapid variations of the integrands due to this condition,the conventional numerical integration scheme in BEM will fail to yield accurate results unless a regularization scheme is adopted,without having to employ a very refined mesh.

Integral regularization schemes that have been implemented in BEM include those by,e.g.,Ma and Kamiya,2001；Niu,Wendland,Wang and Zhou(2005)；Zhou,Niu,Cheng and Guan(2008)；Xie,Zhang,Qin,and Li(2011)；Shiah(2014).Although these approaches have been shown to be very effective in treating 2D problems,they are either analytically complex and/or relatively difficult to implement for general 3D problems.It is worth noting too that the above-mentioned works deal only with isotropic analyses；to treat 3D generally anisotropic problems in elastostatics,these regularization schemes of the boundary integrals will be even more challenging due to the additional mathematical complexity of the anisotropic Green’s function.

For dealing with thin bodies in particular,the methods hitherto reported in the BEM literature are mostly to carry out the integration of the nearly singular integrals analytically[see,e.g.Ye and Liu(1985)；Liu(1987)；Krishnasamy,Rizzo,and Liu(1994)].In Liu(1998),a non-degeneracy approach to analytically evaluate the boundary integrals was proposed.Cruse and Aithal(1993)also proposed a scheme using regularization of integral operators and a theta integration algorithm.Again,these mathematically involved approaches are far from simple to implement,and even more so for 3D generally anisotropic elastic bodies.Another general scheme,often referred to as“self-regularization”,has also been proposed by several authors,see,e.g.Rudolphi(1991)；Sladek,Sladek,and Tanaka(1993)；Matsumoto and Tanaka(1993)；Cruse and Richardson(1996)；Richardson and Cruse(1999)；Frangi and Guigianni(2001)；Dong and Atluri(2012,2013)；He and Tan(2013).For treating 3D thin isotropic bodies,He and Tan(2013)proposed a self-regularized displacement-BIE,employing similar notions as for the self-regularized traction-BIE presented by Cruse and Richardson(1996)；Richardson and Cruse(1999).This proposed self-regularization scheme does not require the C1,αHolder continuity of the boundary element being used.It was found in that study that for 3D,this requirement must be met for the self-regularized traction-BIE if accuracy of the numerical solution is to be assured.The big advantage of the self-regularization procedure of the displacement-BIE as proposed by He and Tan(2013)is the relative simplicity in its implementation.The aim of the present study is to extend this scheme to treat 3D generally anistotropic thin elastic bodies.In the present work,the fundamental solutions in the displacement-BIE are represented by the formulations of double-Fourier series as proposed by Shiah,Tan,and Wang(2012)；Tan,Shiah and Wang(2013).The use of the double-Fourier series representation for the anisotropic Green’s function and its derivatives has been shown to significantly improve the numerical efficiency of the BEM algorithm.A brief review of this will be presented next,followed similarly of the self-regularization scheme of He and Tan(2013).The successful implementation and verification will then be demonstrated by some numerical examples.

Figure 1:Unit circle on the oblique plane at the field point Q.

2 Displacement-BIE for 3D anisotropic elasticity

As has been well established in BEM for 3D elastostatics,the boundary integral equation relating the displacements ujand tractions tjat the surface S of the homogeneous elastic domain can be written in the indicial notation as

whereCij(P)is the free term of the source point P；and,Uij(P,Q)≡Uand Tij(P,Q)≡Trepresent the fundamental solutions of displacements and tractions,respectively,inthexi-directionatthefieldpointQduetoaunitloadappliedinthexj-direction at P in a homogeneous infinite elastic body.The Green’s function,U,for a 3D anisotropic medium was first derived by Lifschitz and Rosenzweig(1947)in the form of a line integral around a unit circle about the source point on an oblique plane,with the integrand containing the Christoffel matrix defined in terms of elastic constants.Since then,there have been numerous efforts to reformulate it into simpler or more explicit analytical forms.As discussed in Shiah,Tan,and Lee(2008)and Tan,Shiah,and Lin(2009),the fully explicit expression forUderived by Ting and Lee(1997)is very well suited indeed for implementation in BEM.Referring to Figure 1,letnandmbe two mutually perpendicular unit vectors on the oblique plane at Q normal to the position vector x；the vectors[n,m,x/r]form a right-angle triad.In the spherical coordinate system as shown,the Green’s function of an anisotropic material with stiffness?Cijksmay be expressed as

where r represents the radial distance between the source point P and the field point Q；all the other parameters are given by

In Equation(3a),the Stroh eigenvalues,pi,are the roots of the sextic equation,obtained from setting|κ |=0.They appear as three pairs of complex conjugates,expressed as pv= αv+iβv,where βv＞ 0 andi= √(-1).The conjugate of the eigenvalue is represented with an overbar,namely,ˉpi.Carrying out the appropriate algebraic operations,the quantity?Γ(n)in Equation(3b)can be shown to be as follows:

From the above,the computations for obtaining the numerical values ofU(x)can be seen to be relatively straightforward.The only numerical scheme required is to obtain the roots of the sextic equation for the Stroh’s eigenvalues.

Unlike the fundamental solution for displacements,the explicit expressions of the derivatives ofU(x)are significantly more complex in form and their computations are significantly more involved,as they contain very high order tensors.Taking advantage of the periodic nature of the fundamental solution when expressed in spherical coordinates,the present lead authors have proposed that the Green’s function and its derivatives be represented by a double-Fourier series[see,Tan,Shiah,and Wang(2013)；Shiah and Tan(2013)].This alternative representation ofU(x)and its derivatives has been demonstrated to be computationally more efficient than using the original form of Ting and Lee’s(1997)solution,besides being relatively simpler to implement.In the double-Fourier series form,U(x)can be written as follows:

Taking spatial derivatives of Equation(5)in the spherical coordinate system yields

Computations ofT≡ Tijare required in the boundary integral equation；they can be carried out using

Equation(1)can be solved for the unknown displacements and tractions on the surface of the body by the usual collocation process of BEM analysis.The selfregularization of this equation by the scheme proposed by He and Tan(2013)will next be discussed as it is extended to 3D general anisotropy.

3 Self-regularization of the displacement-BIE

The quantitiesUandThave singularities of orders r-1and r-2,respectively.The element sub-division process for dealing with singular integrations still applies,when the source point belongs to the element being integrated.However,the problemofnearlysingularintegrationariseswhenthesourcepointonthesurfaceisvery near the opposite surface of a thin body.The self-regularization scheme proposed by He and Tan(2013)is now extended to BEM for 3D general anisotropy.In this scheme,some terms corresponding to a constant stress state need to be subtracted from the displacement-BIE in Equation(1).For a reference point?P on the surface of the body,the expression for the constant stress state can be expressed as

where uLi(P?,Q)and tiL(P?,Q)represent the displacements and tractions at Q under the constant-stress state atP?,respectively.In Equation(10),uLj(P?,PQ)andare given by

Applying the generalized Hooke’s law,Equation(11c)is rewritten as

Subtracting Equation(10)from Equation(1)yields the following self-regularized displacement-BIE:

When using the constant-stress terms defined by Equations(11)and(12),the displacement gradients atP? need fi rst be determined.In a similar manner as presented in He and Tan(2013),the tractions and displacement gradients in the local coordinate system,denoted by ui,j=?ui/?ξj,for general anisotropy can be expressed in the following matrix form:

where the non-zero elements of Amncan be shown to be given by

In Equation(15),Cijis the contracted notation ofrepresents?xi/?ξj,given by

By substituting Equations(17)and(18)into Equations(11a)and(12),the displacements and tractions corresponding to the constant-stress state can now be expressed in terms of the corresponding nodal quantities as follows:

4 Numerical examples

Three example problems are presented here to demonstrate the veracity of the formulation discussed in Section 3 and its successful implementation into a BEM computer program for 3D general anisotropy based on the Fourier series representations of the Green’s function and its derivatives.They are:(A)a thin plate subjected to direct stresses；(B)a thin-walled cylinder subjected to internal pressure；and(C)a thin-walled sphere subjected to internal pressure.For the analyses,crystal alumina(Al2O3)is chosen as the material in all these examples.It has the following stiffness coefficients defined in its principal directions:

To generate generally anisotropic properties for the purpose of the demonstration,the material principal axes are rotated counterclockwise about the x1-,x2-,and x3-axis in succession by 30°,45°and 60°,respectively,to yield

Quadratic isoparametric boundary elements are used throughout.The solutions to these problems were also obtained using the commercial code.ANSYS,based on thefiniteelementmethod(FEM),forthepurposeofverificationoftheBEMresults.Example(A)

The first example deals with a thin square(2L×2L)plate with thicknesst=0.005L as shown in Figure 2.As for the boundary conditions,all the three planes at x1=0,x2=0,and x3=0 are restrained from motion in their respective normal directions,while the remaining planes are subjected to a uniform tensile stress of magnitude P at x1=2L；a uniform tensile stress of magnitude 2P at x2=2L；and uniform compressive stress of magnitude 3P at x3=0.005L as depicted in the figure.Also shown is the BEM mesh employed which has 48 elements and 146 nodes.For the ANSYS-FEM analysis,6400 Shell-181 elements with 45603 nodes were used.To compare the results from the two different numerical methods,the resultant displacement,at points along the centerline at x1=L on the top surface(x3=0.005L)are calculated.Similarly,the von-Mises equivalent stress,σo,at points along the centerline at x1=L on the bottom surface(x3=0)is also determined.Figure 3 and Figure 4 show the variations of the normalized resultant displacementsandequivalentstressescalculatedalongtheselocations,respectively.To demonstrate also the inadequacy of the conventional BEM to treat thin bodies,this problem was also solved using the conventional,un-regularized displacement-BIE,i.e.Eq.(1),employing the same BEM mesh.The comparison of the results for the normalized values of the resultant displacement and equivalent stress at the nodes is shown in Table 1 and Table 2,respectively.From the results,it can be clearly seen that the solutions obtained from the self-regularized displacement-BIE are in excellent agreement with those obtained from the ANSYS analysis,while the un-regularized BIE fails to give reasonable values.

Figure 2:BEM mesh for a thin plate subjected to direct stresses-Example(A).

Figure 3:Variation of the normalized resultant displacement along the centerline at x1=L on the top surface-Example(A).

Figure4:Variationofthenormalizedequivalentstressalongthecenterlineatx1=L on the bottom surface-Example(A).

Table 1:Numerical values of the normalized resultant displacement along the centerline at x1=L on the top surface-Example(A)

Table 2:Numerical values of the normalized equivalent stress along the centerline at x1=L on the bottom surface-Example(A)

Example(B)

Figure 5:A thin-walled cylinder subjected to internal pressure-Example(B).

The dimensions of the thin-walled alumina cylinder for the second example are shown in Figure 5.The cylinder is subjected to a uniform internal pressure P,with both its ends completely fixed.Figure 6 shows the BEM mesh employed for the analysis；it has 704 elements and 2112 nodes.Since the analytical solution to this typical 2D isotropic problem is available,the analysis was first conducted using the algorithm for anisotropy but with elastic constants corresponding to isotropy with Young’s modulus=1000(units)and Poisson’s ratio=0.3.Table 3 lists the numerical values of the normalized hoop stress at some sample nodes on the inner surface of the axial plane x3=L as calculated by the self-regularized BIE and also by the conventional,un-regularized BIE.Also shown in the table in parentheses are the percentage deviations of the BEM solutions from the exact analytical solution from Lame theory.As can be seen,the results obtained from the selfregularized BIE agree very well with the analytical solution,while those from the un-regularized scheme are totally erroneous.Next,the generally anisotropic problem with the elastic constants given in Equation(21)was analyzed.For comparison,the anisotropic analysis was also performed using ANSYS,the commercial FEM software.The ANSYS model is relatively refined with 5120 Shell-181 elements and 36480 nodes as shown in Figure 7.As in the previous analysis of the degenerate case of isotropy,the hoop stress at the nodes at the inner circumference of the cylinder at the mid-axial plane,were obtained,as listed in Table 4.The agreement of the results obtained using the self-regularized BEM algorithm and those from the FEM is very good indeed,the percentage differences between them are shown in parentheses in the table.The same cannot be said of those results obtained with the conventional BEM.

Example(C)

The problem considered is a thin-walled sphere with inner radius ri=L and thick-nesst=0.005L.AsinExample(B),theanalysiswasfirstperformedfortheisotropic case using Young’s modulus=1000(units)and Poisson’s ratio=0.3.The inner surface of the sphere is subjected to uniform pressure P,while the outside surface is traction free.Appropriate nodal constraints were applied to preclude rigid body motion of the numerical model.The exact analytical solution,based again on Lame theory,for the circumferential stress is also readily available.Figure 8 shows the BEM mesh for the analysis；there are 256 elements with 708 nodes.The hoop stresses at the nodes on the inner surface around the equator are computed by the both regularized and un-regularized BIE algorithms.Listed in Table 5 are the computed results for comparison with the exact analytical solution.As expected,consistent,accurate results are obtained by the self-regularized BIE but not for the un-regularized scheme.The anisotropic material properties in Equation(21)were next used in the analysis.For the boundary condition of this anisotropic analysis,the same uniform internal pressure P was applied,but the exterior surface of the sphere was constrained in all directions from displacements.For verification,the

problem was also analyzed by ANSYS with the FEM mesh shown in Figure 9,which has 4167 Shell-181 elements and 29140 nodes.The von Mises equivalent stress,σ0,atsomesamplenodesaroundtheoutersurfaceoftheequatorarecomputed.As before,the problem was also solved using the conventional un-regularized scheme of the BEM.Table 6 shows the comparison of the computed results from the numerical analyses.The results by the self-regularized BIE are again in excellent agreement with those from the ANSYS analysis；this is not the case with the un-regularized BIE.

Table 3:Normalized hoop stress around the inner circumference of the thinwalled isotropic cylinder at the midaxial plane-Example(B)

Table 4:Normalized hoop stress around inner circumference of the anisotropic thin-walledaluminacylinderatthemidaxial plane-Example(B)

Figure 6:BEM mesh the thin-walled cylinder-Example(B).

Figure 7:Figure 7:FEM mesh of the thin-walled cylinder by ANSYS-Example(B).

Table 5:Normalized circumferential stress at the outer circumference of the thin-walled isotropic sphere-Example(C)

Table 6:Normalized equivalent stress at the outer circumference around the equator of the thin-walled anisotropic sphere-Example(C)

Figure 8:BEM mesh the thin-walled sphere-Example(C).

Figure 9:FEM mesh of the thin-walled sphere by ANSYS-Example(C).

5 Conclusions

Thin-bodies are quite commonly encountered in engineering applications.For these problems,the conventional BEM is not a very suitable tool for the stress analysis because of the nearly-singular integrals that arise in the boundary integral equation(BIE).A self-regularization scheme proposed by He and Tan(2013)for 3D isotropic elastostatics has been extended to 3D generally anisotropic elasticity in the present study.The self-regularization is performed on the displacement-BIE,hence there is no requirement for C1,αHolder continuity between the boundary elements.The scheme has been found to be relatively simple to implement and it has been successfully carried out on an existing BEM code,based on the quadratic isoparametric formulation,for 3D general anisotropic elasticity which utilizes double-Fourier series representations of the Green’s function and its derivatives.This has been demonstrated by three example problems where the numerical results are verified by similar analysis performed using relatively refined meshes with the commercial FEM code ANSYS.

Acknowledgement:The authors gratefully acknowledge the financial support for this work by the Ministry of Science and Technology,Taiwan(No.102-2221-E-006-290-MY3)and the National Science and Engineering Research Council of Canada.

Cruse,T.A.；Aithal,R.(1993):Non-singular boundary integral equation implementation.International Journal for Numerical Methods in Engineering,vol.36,pp.237-254.

Cruse,T.A.；Richardson,J.D.(1996):Non-singular Somigliana stress identities in elasticity.International Journal for Numerical Methods in Engineering,vol.39,pp.3273-3304.

Dong,L.；Atluri,S.N.(2012):SGBEM(using non-hyper-singular traction-BIE)and super elements,for non-collinear fatigue-growth analyses of cracks in stiffened panels with composite-patch repairs.CMES-Computer Modeling in Engineering&Sciences,vol.89,pp.415-456.

Dong,L.；Atluri,S.N.(2013):SGBEM voronoi cells(svcs),with embedded arbitrary-shaped inclusions,voids,and/or cracks,for micromechanical modeling of heterogeneous materials.CMC:Computers,Materials&Continua,vol.33,pp.111-154.

Granados,J.J.；Gallego,R.(2001):Regularization of nearly hypersingular integrals in the boundary element method.Engineering Analysis with Boundary Elements,vol.25,pp.165-184.

He,M.G.；Tan,C.L.(2013):A self-regularization technique in boundary element method for 3-D stress analysis.CMES-Computer Modeling in Engineering&Sciences,vol.95,pp.317-349.

Krishnasamy,G.；Rizzo,F.J.；Li,Y.J.(1994):Boundary integral equations for thin bodies.International Journal for Numerical Method in Engineering,vol.37,pp.107-121.

Lifshitz,I.M.；Rozenzweig,L.N.(1947):Construction of the Green Tensor for the fundamental equation of elasticity theory in the case of unbounded elastically anisotropic medium,Zh.Eksp.Teor.Fiz.,vol.17,pp.783-791.

Liu,Y.J.(1987):Elastic stability analysis of thin plate by the boundary element method-a new formulation.Engineering Analysis with Boundary Elements,vol.4,pp.160-164.

Liu,Y.J.(1998):Analysis of shell-like structures by the boundary element method based on 3-D elasticity:formulation and verification.International Journal for Numerical Methods in Engineering,vol.41,pp.541-558.

Ma,H.；Kamiya,N.(2001):A general algorithm for accurate computation of field variables and its derivatives near the boundary in BEM.Engineering Analysis with Boundary Elements,vol.25,pp.833-841.

Matsumoto,T.；Tanaka,M.(1993):Boundary stress calculation using regularized boundary integral equation for displacement gradients.International Journal for Numerical Methods in Engineering,vol.36,pp.783-797.

Rudolphi,T.J.(1991):The use of simple solutions in the regularization of hypersingular boundary integral equations.Mathematical and Computer Modeling,vol.15,pp.269-278.

Shiah,Y.C.；Tan,C.L.；Lee,V.G.(2008):Evaluation of explicit-form fundamental solutions for displacements and stresses in 3D anisotropic elastic solids.CMES-Computer Modeling in Engineering&Sciences,vol.34,pp.205-226.

Shiah,Y.C.；Tan,C.L.；Wang,C.Y.(2012):Efficient computations of the Green’s function and its derivatives for three-dimensional anisotropic elasticity in BEM analysis.Engineering Analysis with Boundary Elements,vol.36,pp.1746-1755.

Shiah,Y.C.；Tan,C.L.；Chen,Y.H.(2013):Efficient BEM stress analysis of 3D generally anisotropic elastic solids with stress concentrations and cracks,CMESComputer Modeling in Engineering&Sciences,vol.96,no.4,pp.243-257.

Shiah,Y.C.(2014):3D elastostatic boundary element analysis of thin bodies by integral regularizations,Journal of Mechanics,vol.31,issue 5,pp.533-543.

Sladek,V.；Sladek,J.；Tanaka,M.(1993):Regularization of hypersingular and nearly singular integrals in potential theory and elasticity.International Journal for Numerical Methods in Engineering,vol.36,pp.1609-1628.

Tan,C.L.；Shiah,Y.C.；Lin,C.W.(2009):Stress analysis of 3D generally anisotropic solids usin the boundary element method.CMES-Computer Modeling in Engineering&Sciences,vol.41,pp.195-214.

Tan,C.L.；Shiah,Y.C.；Wang,C.Y.(2013):Boundary element elastic stress analysis of 3D generally anisotropic solids using fundamental solutions based on Fourier series.International Journal of Solids&Structures,vol.50,pp.2701-2711.

Ting,T.C.T.；Lee,V.G.(1997):The three-dimensional elastostic Green’s function for general anisotropic linear elastic solid.Q.J.Mech.Appl.Math.,vol.50,pp.407-426.

Niu,Z.；Wendland,W.L.；Wang,X.；Zhou,H.(2005):A semi-analytical algorithm for the evaluation of the nearly singular integrals in three-dimensional boundary element methods.Computer Methods in Applied Mechanics and Engineering,vol.194,no.9-11,pp.1057-1074.

Ye,T.Q.；Liu,Y.J.(1985):Finite def l ection analysis of elastic plate by the boundary element method.Applied Mathematical Modeling,vol.9,pp.183-188.

Zhou,H.L.；Niu,Z.R.；Cheng,C.Z.；Guan,Z.W.(2008):Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for anisotropic potential problems.Computers and Structures,vol.86,pp.1656-1671.

1Department of Aeronautics and Astronautics,National Cheng Kung University,Tainan 701,Taiwan,ROC.

2Corresponding author；Email:choonlai.tan@carleton.ca

3Department of Mechanical&Aerospace Engineering,Carleton University,Ottawa,Canada K1S 5B6