Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.019160

ARTICLE

Weakly Singular Symmetric Galerkin Boundary Element Method for Fracture Analysis of Three-Dimensional Structures Considering Rotational Inertia and Gravitational Forces

Shuangxin He1, Chaoyang Wang1, Xuan Zhou1,*, Leiting Dong1,* and Satya N. Atluri2

1School of Aeronautic Science and Engineering, Beihang University, Beijing, China
2Department of Mechanical Engineering, Texas Tech University, Lubbock, USA
*Corresponding Authors: Xuan Zhou. Email: zhoux@buaa.edu.cn; Leiting Dong. Email: ltdong@buaa.edu.cn
Received: 06 September 2021; Accepted: 26 November 2021

Nomenclature

1  Introduction

The Symmetric Galerkin Boundary Element Method (SGBEM) [1–3] has gained increasing popularity in fracture and crack-growth analysis of solid structures due to its attractive features of symmetric coefficient matrices, weak-singularity, and that only boundary & crack-surface elements are needed. The papers by Bonnet et al. [3–5] are devoted to the formulation, numerical evaluation and implementation of SGBEM. Atluri et al. [6–9] utilized a simple and straightforward methodology to develop regularized traction Boundary Integral Equations (tBIE) for two and three-dimensional linear-elastic solids containing cracks, and also developed weakly-singular SGBEMs for the fracture and fatigue analysis of various complex structures. However, for the fracture mechanics problems such as turbine discs and turbine blades of aircraft engines, concrete gravity dam, etc., SGBEM may lose its advantages, because evaluation of domain integral terms resulting from body forces such as rotational inertia and gravitational loads leads to the meshing of the interior of the domain. For this reason, a method to evaluate such domain integral terms using only boundary meshes, is desired to efficiently analyze cracked structures considering body forces with SGBEM.

For the conventional collocation boundary element method based on Somigliana’s identity for the displacement vector, a few methods were developed for this purpose. Considering centrifugal loads presented in rotating gas turbines, Cruse et al. [10,11] transformed domain integrals to boundary integrals by utilizing the divergence theorem. By making use of the Galerkin vector or the Green’s second identity, Danson [12] transformed the volume integral terms to boundary integral terms, for three kinds of body forces, i.e., gravitational loads, the rotational inertia and steady-state thermal loads. Gao [13] also developed a radial integration technique and applied it to deal with various body forces. Brebbia et al. [14] developed the dual reciprocity method [15] which converts the associated domain integrals into boundary integrals by using a series of basis functions to approximate the body force fields. Brebbia et al. [16] extended the idea of dual reciprocity and proposed another approach, multiple reciprocity method.

Different from the conventional collocation boundary element method [17–19] based on the Somigliana’s identity, formulations of SGBEM [5,8,20] result in weak-form displacement Boundary Integral Equations (dBIE) and weak-form traction Boundary Integral Equations (tBIE). As a matter of fact, the domain integrals caused by body forces appear both in dBIE and tBIE. Moreover, it is beneficial to use tBIE to derive weak-form equations on crack-surfaces, where displacement discontinuities are to be solved as unknowns [5]. Thus, if SGBEM is utilized for linear fracture analysis of cracked structures, while for domain integrals appearing in dBIE, one may refer to the above-mentioned transformation techniques, the treatment for domain integral terms appearing in tBIE needs further study.

This paper presents the weakly singular traction boundary integral equation for solids undergoing rotational inertia and gravitational Loads. By using the divergence theorem (div) or the radial integration method (RIM), domain integrals induced by rotational inertia or gravitational forces are transformed into boundary integrals correspondingly. The derived formulas show that these transformed boundary integral terms have no influence upon the coefficient matrix of SGBEM, but only affect the right-hand-side vector. The transformed boundary integral terms derived by the divergence theorem and radial integration method, possessing 1/r singularity, is weakly singular. Numerical examples demonstrate that only a few Gauss points are sufficient to evaluate boundary integrals. The developed SGBEM with only weakly-singular boundary integrals are thus applied to simulate various examples of 3D solids with/without considering rotational inertia and gravitational loads.

This paper is organized as follows. In Section 2, transformation from domain integrals induced by gravitational and rotational inertia forces to the boundary integrals by div or RIM respectively is carried out. Some numerical examples for solids undergoing rotational inertia or gravitational loads are presented in Sections 3 and 4 with and without cracks correspondingly. In Section 5, we complete this paper with some concluding remarks.

2  Weakly Singular Galerkin Boundary Integral Equations and Boundary Element Method with Rota- tional Inertia and Gravitational Loads

Consider a linear elastic, homogenous and isotropic solid undergoing an infinitesimal elasto-static deformation, as shown in Fig. 1. Ω is the solution domain of the problem with the boundary ∂Ω. ξ represents the field point at a generic location in Cartesian coordinates. x is the source point of the 3D Kervin’s solution [21] where a unit load in an arbitrary direction p is applied. The displacement fundamental solution uj∗p in the j direction corresponding to this unit load and other kernel functions Gij∗p,φij∗p,Htbpq∗ derived by uj∗p are listed in the appendix. One may also refer to other forms of these kernel functions in [5,20].

Figure 1: A solution domain with source point x and field point ξ

The symmetric Galerkin formulations of displacement and traction Boundary Integral Equations (d & tBIE) for linear elastic solids can be found in [8]. The derivation of the conventional boundary element method and SGBEM [5,8,20] generally ignored body forces. Here, the domain integrals considering body forces are added in the equations:

−∫∂Ω12δtp(x)up(x)dSx=−∫∂Ωδtp(x)dSx∫Ωfj(ξ)uj∗p(ξ−x)dΩξ−∫∂Ωδtp(x)dSx∫∂Ωtj(ξ)uj∗p(ξ−x)dΓξ−∫∂Ωδtp(x)dSx∫∂ΩDi(ξ)uj(ξ)Gij∗p(ξ−x)dΓξ−∫∂Ωδtp(x)dSx∫∂Ωuj(ξ)ni(ξ)φij∗p(ξ−x)dΓξ, (1)

∫∂Ω12δub(x)tb(x)dSx=−∫∂Ωδub(x)na(x)dSx∫Ωfj(ξ)σab∗j(ξ−x)dΩξ−∫∂ΩDt(x)δub(x)dSx∫∂Ωtj(ξ)Gtb∗j(ξ−x)dΓξ+∫∂Ωδub(x)na(x)dSx∫∂Ωtj(ξ)φab∗j(ξ−x)dΓξ−∫∂ΩDt(x)δub(x)dSx∫∂ΩDp(ξ)uq(ξ)Htbpq∗(ξ−x)dΓξ. (2)

In the above two equations, if the domain integral or boundary integral is with respect to the field point, the integral domain is denoted by Ωξ or Γξ, respectively; otherwise for source point, the integral domain is denoted by Sx. up(x) and tp(x) are the displacement and the traction at the source point, respectively. δ is the variational symbol which is used to import the Galerkin weight function. fj(ξ) is the body force per unit volume. ni(ξ) is the component of outward unit normal at a field point on the boundary. Dt is a surface tangential operator:

Dt(ξ)=nr(ξ)erst∂∂ξs (3)

where erst is the permutation coefficient defined by e123=e231=e312=1; e321=e213=e132=−1; erst=0 if any two of the indices are identical.

In this paper, the domain integral:

I=∫Ωfj(ξ)σab∗j(ξ−x)dΩξ (4)

appearing in traction boundary integral Eq. (2) considering rotational inertia and gravitational loads is transformed into weakly singular boundary integral, using the divergence theorem or the radial integration method.

The radial integration method is introduced here briefly. For further details, one may refer to [13]. Domain integral on the left-hand-side of Eq. (5) with a general function f(ξ) may be written in Cartesian coordinate system (x1,x2,x3) or in spherical coordinate system (r,θ,ϕ) with the origin at the source point P shown in Fig. 2.

Figure 2: Cartesian and spherical coordinate systems

In the spherical coordinate system

∫Ωf(ξ)dΩ=∫02π∫0π∫0r(θ,ϕ)f(r,θ,ϕ)r2drsin⁡θdθdϕ=∫02π∫0πF(θ,ϕ)sin⁡θdθdϕ (5)

where

F(θ,ϕ)=∫0r(θ,ϕ)f(r,θ,ϕ)r2dr. (6)

In the spherical coordinate system, the area of infinitesimal element dS on the spherical surface can be expressed as

dS=r2sin⁡θdθdϕ. (7)

If the field point is on the boundary Γ of domain Ω, geometric projective transformation can be established between the spherical surface infinitesimal element dS and the real boundary surface infinitesimal element dΓ shown in Fig. 3.

dS=rirnidΓ, (8)

where ni is the component of outward unit normal of field point on the real boundary surface dΓ, ri is the Cartesian component of r, i.e.,

ri=ξi−xi. (9)

Figure 3: Spherical surface dS and real boundary dΓ

By some derivations, the domain integral can be rewritten as

∫Ωf(ξ)dV=∫∂Ω1r2∂r∂nF(r)dΓ (10)

where F(r) is evaluated by a radial integration of R2f(ξ) on the segment linking the source point and the field point, i.e.,

F(r)=∫0rR2f(ξ)dR. (11)

∂r/∂n is the directional derivative at the field point on the boundary, which may be expressed as

∂r∂n=r,ini (12)

where (),i denotes the partial differentiation with respect to the Cartesian component of field point if not otherwise stated. And r,i can be expressed as

r,i=∂r∂ξi=rir=−∂r∂xi. (13)

Some useful formulas related to r (r≠0) are listed as follows:

r=riri (14)

r,ir,i=1(15)

r,ij=1r(δij−r,ir,j) (16)

r,ii=2r (17)

r,ijk=−1r2(r,iδjk+r,jδik+r,kδij−3r,ir,jr,k) (18)

(1r),i=−1r2r,i (19)

In Subsections 2.1 and 2.2, the domain integral terms with rotational inertia and gravitational loads in tBIE are transformed into weakly singular boundary integral terms by two methods of divergence theorem and radial integration method, respectively.

2.1 Transformation of Domain Integrals with Gravitational Loads to Boundary Integrals

Consider a solid body with a constant mass density ρ, and a constant gravitational field gi=const. The body force will also be constant, where

fi=ρgi=const. (20)

In this section, the body force fj(ξ) in Eq. (4) is defined as gravitational force. The purpose of this section is transforming the domain integral of Eq. (4) considering gravitational force into the boundary integral. Note that σab∗j(ξ−x) in Eq. (4) is the stress field of Kelvin’s solution:

σab∗j(ξ−x)=18π(1−ν)r2[(1−2ν)(δabr,j−δajr,b−δbjr,a)−3r,ar,br,j] (21)

where ν is the Poisson’s ratio; δab is the Kronecker Delta.

Thus, the constant gravity force fi can be taken outside the integral in Eq. (4). Then we get

∫Ωρgjσab∗j(ξ−x)dΩξ=ρgj∫Ω18π(1−ν)r2[(1−2ν)(δabr,j−δajr,b−δbjr,a)−3r,ar,br,j]dΩξ. (22)

2.1.1 Using Divergence Theorem to Transform Domain Integrals with Gravitational Forces

Substitution of Eq. (19) into Eq. (26), we have

∫Ωρgjσab∗j(ξ−x)dΩξ=ρgj18π(1−ν)∫Ω[2νδab(−1r2r,j)+2(1−ν)δaj(−1r2r,b)+2(1−ν)δbj(−1r2r,a)−r,abj]dΩξ. (23)

Substituting Eq. (22) into Eq. (27), we have

∫Ωρgjσab∗j(ξ−x)dΩξ=ρgj18π(1−ν)∫Ω[2νδab(1r),j+2(1−ν)δaj(1r),b+2(1−ν)δbj(1r),a−r,abj]dΩξ. (24)

Using divergence theorem and Eq. (16), we can get that

∫Ωρgjσab∗j(ξ−x)dΩξ=ρgj18π(1−ν)∫∂Ω1r[(2ν−1)δabnj(ξ)+2(1−ν)δajnb(ξ)+2(1−ν)δbjna(ξ)+r,ar,bnj(ξ)]dΓξ. (25)

Note that, a singularity of 1/r appears in the boundary integral of Eq. (29). This integral is weakly-singular [8], thus Cauchy principal value integral [22] does not need to be taken into account. The numerical integration method to evaluate this weakly-singular integral is stated briefly in Section 2.3.3.

2.1.2 Using the Radial Integration Method to Transform Domain Integrals with Gravitational Forces

Using radial integration method, Eq. (26) can be rewritten as

∫Ωρgjσab∗j(ξ−x)dΩξ=ρgj∫∂Ω∂r∂n1r2F(r)dΓξ (26)

where

F(r)=∫0rR218π(1−ν)R2[(1−2ν)(δabR,j−δajR,b−δbjR,a)−3R,aR,bR,j]dR. (27)

From Eq. (13) and Fig. 2, one can find that R,i is the cosine between r and coordinate axis i, i.e., R,i=r,i. Thus, R,i can be taken out of this radial integral Eq. (31) directly. Substitution of Eq. (31) into Eq. (30) gives

∫Ωρgjσab∗j(ξ−x)dΩξ=ρgj∫∂Ω18π(1−ν)1r∂r∂n[(1−2ν)(δabr,j−δajr,b−δbjr,a)−3r,ar,br,j]dΓξ. (28)

Note that, when the field point approaches the source point, ∂r/∂n→0. Singularity of the boundary integral in Eq. (32) therefore is weaker than that in Eq. (29).

2.2 Transform Domain Integrals with Rotational Inertia to Boundary Integrals

About an analytical expression of the rotational inertial force in detail, one may refer to [19]. Here we introduce it briefly. Consider a solid body of uniform mass density ρ rotating about one axis with angular velocity ωi. For simplicity and without loss of the generality, we consider that the axis of rotation passes through the origin of Cartesian coordinate system shown in Fig. 4.

Figure 4: The rotational axis passing through the origin of Cartesian coordinate system

By the D’Alembert’s principle, body force resulting from the rotational inertia is

f(ξ)=−ρω×(ω×ξ). (29)

Eq. (33) may be written in index notation as

fj(ξ)=−ρejqkωqekpiωpξi=hjiξi (30)

where

hji=−ρejqkωqekpiωp. (31)

Note that hji is constant and can be described in a more straightforward way:

[hji]=ρ[ω22+ω32−ω1ω2−ω3ω1−ω1ω2ω32+ω12−ω2ω3−ω3ω1−ω2ω3ω12+ω22]. (32)

Then this dynamic problem can be treated as an elastostatics problem. Using Eqs. (4), (25) and (34), we get

∫Ωhjiξiσab∗j(ξ−x)dΩξ=∫Ωhjiξi18π(1−ν)r2[(1−2ν)(δabr,j−δajr,b−δbjr,a)−3r,ar,br,j]dΩξ. (33)

2.2.1 Using Divergence Theorem to Transform Domain Integrals with Inertial Force

Similar to the derivation of Eq. (28), the inertial force domain integrals with the rotational inertia can be written as

∫Ωhjiξiσab∗j(ξ−x)dΩξ=∫Ω18π(1−ν)r2hji[2νδabξi(1r),j+2(1−ν)δajξi(1r),b+2(1−ν)δbjξi(1r),a−ξir,abj]dΩξ. (34)

Substituting Eqs. (39) and (40) into Eq. (38) and using Eq. (17),

(ξi1r),j=δij1r+ξi(1r),j (35)

(ξir,ab),j=δijr,ab+ξir,abj (36)

We get

∫Ωhjiξiσab∗j(ξ−x)dΩξ=∫Ωhji18π(1−ν)[2νδab(ξi1r),j−νδabδijr,kk+2(1−ν)δaj(ξi1r),b−(1−ν)δajδbir,kk+2(1−ν)δbj(ξi1r),a−(1−ν)δbjδair,kk+δijr,ab−(ξir,ab),j]dΩξ. (37)

Then using the divergence theorem, we get

∫Ωhjiξiσab∗j(ξ−x)dΩξ=∫∂Ω18π(1−ν)1r[2νδabnj(ξ)hjiξi−νδabgiirr,mnm(ξ)+2(1−ν)nb(ξ)haiξi+2(1−ν)na(ξ)hbiξi−(1−ν)rr,mnm(ξ)(hab+hba)+rhiir,anb(ξ)−(δab−r,ar,b)nj(ξ)hjiξi]dΓξ. (38)

Note that the boundary integrals in Eq. (42) have the property of 1/r weak-singularity.

2.2.2 Using Radial Integration method to Transform Domain Integrals with Inertial Force

Using radial integration method, Eq. (37) may be rewritten as

∫Ωhjiξiσab∗j(ξ−x)dΩξ=∫∂Ω1r2∂r∂nF(r)dΓξ (39)

F(r)=∫0rR2hjiξi18π(1−ν)R2[(1−2ν)(δabR,j−δajR,b−δbjR,a)−3R,aR,bR,j]dR. (40)

As is mentioned above, R,j can be taken outside the integral directly. Note that, F(r) is the radial integral about the field point ξ. Substitution of Eq.(9) into Eq. (44) gives

F(r)=hji18π(1−ν)[(1−2ν)(δabr,j−δajr,b−δbjr,a)−3r,ar,br,j]∫0rR(R,i+xiR)dR. (41)

Note that, for radial integral F(r), source point x is constant. We can directly compute this radial integral. Substitution of Eq. (45) into Eq. (43) gives

∫Ωhjiξiσab∗j(ξ−x)dΩξ=∫∂Ω116π(1−ν)1r∂r∂nhji(ξi+xi)[(1−2ν)(δabr,j−δbjr,a−δajr,b)−3r,ar,br,j]dΓξ. (42)

Eq. (46) is the boundary integral form with the rotational inertia force obtained by the radial integration method.

2.3 Weakly-Singular SGBEM with Numerical Implementation

We have obtained weakly singular boundary integrals transformed from domain integrals considering rotational inertia and gravitational loads by the divergence theorem or radial integration method. In this section, the displacement and traction boundary integral equations considering crack surfaces and rotational inertia and gravitational loads are given. Then numerical evaluation of weakly singular double layer surface integrals by using quadrilateral elements is introduced briefly.

2.3.1 Traction and Displacement BIEs Considering Rotational Inertia and Gravitational Loads by Divergence Theorem

Consider a crack embedded in the domain Ω shown in Fig. 5. The crack surfaces are denoted as SC+ and SC− which are geometrically coincident. The outward normal direction of SC+ is opposite to that of SC−. With the assumption that the traction acting on crack surfaces satisfies that tj++tj−=0, the boundary of the domain Ω can be defined as

∂Ω=Su+St+SC  (43)

where Su is the part of boundary where displacement is known and St is the part of boundary where traction is known. The displacement discontinuity on crack surfaces may be defined as

Δu=u+(x+)−u−(x−) (44)

where u+(x+) is the displacement of point x+ on SC+; u−(x−) is the displacement of point x− on SC−; Δu must be zero around the crack front. Points x+ and x− are geometrically coincident.

Figure 5: Displacement discontinuity in domain Ω

If the weak-form traction boundary integral equation is applied on St, we may get that

12∫Stδub(x)tb(x)dSx+∫Stδub(x)na(x)dSx∫Su+Stρgj18π(1−ν)1r[(2ν−1)δabnj(ξ)+2(1−ν)δajnb(ξ)+2(1−ν)δbjna(ξ)+r,ar,bnj(ξ)]dΓξ+∫Stδub(x)na(x)dSx∫Su+St18π(1−ν)1r[2νδabnj(ξ)hjiξi−νδabhiirr,mnm(ξ)+2(1−ν)nb(ξ)haiξi+2(1−ν)na(ξ)hbiξi−(1−ν)rr,mnm(ξ)(hab+hba)+rhiir,anb(ξ)−(δab−r,ar,b)nj(ξ)hjiξi]dΓξ=−∫StDtδub(x)dSx∫Su+Sttj(ξ)Gtb∗j(ξ−x)dΓξ+∫Stδub(x)na(x)dSx∫Su+Sttj(ξ)φab∗j(ξ−x)dΓξ−∫StDtδub(x)dSx∫Su+StDpuq(ξ)Htbpq∗(ξ−x)dΓξ−∫StDtδub(x)dSx∫ScDpΔuq(ξ)Htbpq∗(ξ−x)dΓξ. (45)

And if the weak-form traction boundary integral equation is applied on the crack surfaces Sc, we may get that

∫ScδΔub(x)tb(x)dSx+∫ScδΔub(x)na(x)dSx∫Su+Stρgj18π(1−ν)1r[(2ν−1)δabnj(ξ)+2(1−ν)δajnb(ξ)+2(1−ν)δbjna(ξ)+r,ar,bnj(ξ)]dΓξ+∫ScδΔub(x)na(x)dSx∫Su+St18π(1−ν)1r[2νδabnj(ξ)hjiξi−νδabhiirr,mnm(ξ)+2(1−ν)nb(ξ)haiξi+2(1−ν)na(ξ)hbiξi−(1−ν)rr,mnm(ξ)(hab+hba)+rhiir,anb(ξ)−(δab−r,ar,b)nj(ξ)hjiξi]dΓξ=−∫ScDtδΔub(x)dΓx∫Su+Sttj(ξ)Gtb∗j(ξ−x)dΓξ+∫ScδΔub(x)na(x)dSx∫Su+Sttj(ξ)φab∗j(ξ−x)dΓξ−∫ScDtδΔub(x)dSx∫Su+StDpuq(ξ)Htbpq∗(ξ−x)dΓξ−∫ScDtδΔub(x)dSx∫ScDpΔuq(ξ)Htbpq∗(ξ−x)dΓξ. (46)

Finally, the weak-form displacement boundary integral equation is applied on the prescribed displacement boundary surfaces Su, we may get that

−12∫Suδtp(x)up(x)dSx+∫Suδtp(x)dSx∫Su+St1+ν4πE[∂r∂nρgp−12(1−ν)ρgjnj(ξ)r,p]dΓξ+∫Suδtp(x)dSx∫Su+St1+ν4πE[hpiξir,jnj(ξ)−rhpini(ξ)−12(1−ν)r,pnj(ξ)hjiξi+12(1−ν)rnp(ξ)hjj]dΓξ=−∫Suδtp(x)dSx∫Su+Sttj(ξ)uj∗p(x,ξ)dΓξ−∫Suδtp(x)dSx∫Su+StDiuj(ξ)Gij∗p(x,ξ)dΓξ−∫Suδtp(x)dSx∫Su+Stuj(ξ)ni(ξ)φij∗p(x,ξ)dΓξ−∫Suδtp(x)dSx∫ScDiΔuj(ξ)Gij∗p(x,ξ)dΓξ−∫Suδtp(x)dSx∫ScΔuj(ξ)ni(ξ)φij∗p(x,ξ)dΓξ. (47)

Eqs. (49)–(51) are the weakly-singular traction and displacement boundary integral equations considering rotational inertia and gravitational loads obtained by using divergence theorem. E and ν are Young’s modulus and Poisson’s ratio of the isotropic solid, respectively. Then we may discretize boundary surfaces ∂Ω into boundary elements. Traction field functions can be written in terms of nodal shape functions as ti=timNm at Su, ti=t¯imNm at St; similarly displacement field functions can be written as ui=u¯imNm at Su, ui=uimNm at St, where an overline denotes that the nodal variables are known. In this way, the discretized traction and displacement SGBEM equations are obtained, and we denote this method as SGBEM-div in this paper.

2.3.2 Traction and Displacement BIEs Considering Rotational Inertia and Gravitational Loads by the Radial Integration Method

Similar to Eqs. (49)–(51), the weakly-singular traction and displace BIE considering rotational inertia and gravitational loads by radial integration method can be written as follows:

12∫Stδub(x)tb(x)dSx+∫Stδub(x)na(x)dSx∫Su+Stρgj18π(1−ν)1r∂r∂n[(1−2ν)(δabr,j−δajr,b−δbjr,a)−3r,ar,br,j]dΓξ+∫Stδub(x)na(x)dSx∫Su+St116π(1−ν)1r∂r∂nhji(ξi+xi)[(1−2ν)(δabr,j−δbjr,a−δajr,b)−3r,ar,br,j]dΓξ=−∫StDtδub(x)dΓx∫Su+Sttj(ξ)Gtb∗j(ξ−x)dΓξ+∫Stδwb(x)na(x)dSx∫Su+Sttj(ξ)φab∗j(ξ−x)dΓξ−∫StDtδub(x)dSx∫Su+StDpuq(ξ)Htbpq∗(ξ−x)dΓξ−∫StDtδub(x)dSx∫ScDpΔuq(ξ)Htbpq∗(ξ−x)dΓξ. (48)