images
Computer Modeling in Engineering & Sciences
images

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

Abstract: The Symmetric Galerkin Boundary Element Method is advantageous for the linear elastic fracture and crack-growth analysis of solid structures, because only boundary and crack-surface elements are needed. However, for engineering structures subjected to body forces such as rotational inertia and gravitational loads, additional domain integral terms in the Galerkin boundary integral equation will necessitate meshing of the interior of the domain. In this study, weakly-singular SGBEM for fracture analysis of three-dimensional structures considering rotational inertia and gravitational forces are developed. By using divergence theorem or alternatively the radial integration method, the domain integral terms caused by body forces are transformed into boundary integrals. And due to the weak singularity of the formulated boundary integral equations, a simple Gauss-Legendre quadrature with a few integral points is sufficient for numerically evaluating the SGBEM equations. Some numerical examples are presented to verify this approach and results are compared with benchmark solutions.

Keywords: Symmetric Galerkin boundary element method; rotational inertia; gravitational force; weak singularity; stress intensity factor

Nomenclature

E : | Young’s modulus
ν : | Poisson’s ratio
μ : | Shear modulus
ρ : | Density
gi : | Gravitational acceleration
ωi : | Angular velocity

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)

12∫Scδub(x)tb(x)dSx+∫Scδ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Γξ+∫Scδ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Γξ=−∫ScDtδub(x)dΓx∫Su+Sttj(ξ)Gtb∗j(ξ−x)dΓξ+∫Scδwb(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Γξ. (49)

−12∫Suδtp(x)up(x)dSx+∫Suδtp(x)dSx∫Su+St1+ν16πE(1−ν)∂r∂n[(3−4ν)ρgp(ξ)+r,pr,jρgj(ξ)]dΓξ+∫Suδtp(x)dSx∫Su+St1+ν24πE(1−ν)r,mnm(ξ)[(3−4ν)δpj+r,pr,j]hji(ξi+12xi)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Γξ. (50)

By the same discretization procedure mentioned above for Eqs. (52)–(54), the SGBEM equations obtained by radial integration method can be obtained, and we denote this as SGBEM-RIM in this paper.

It can be seen that, for Eqs. (49), (50) using the divergence theorem, there exists 1/r singularity in boundary integral terms containing rotational inertia and gravitational loads; while for Eqs. (52), (53) using the radial integration method, there exists 1/r⋅∂r/∂n in boundary integral terms containing rotational inertia and gravitational loads. As is mentioned above, when the field point approaches the source point, ∂r/∂n→0. In other words, by the radial integration method, the obtained boundary integral terms may have weaker singularity compared with those obtained by the divergence theorem.

2.3.3 Numerical Evaluation of Weakly-Singular Double Surface Integrals Using Quadrilateral Elements

In this paper, 8-noded quadrilateral isoparametric elements are selected for the numerical implementation, and quarter-point singular quadrilateral elements with two mid-side nodes shifted towards the crack front as shown in Fig. 6 are adopted at the crack front. For the numerical evaluation of double surface integrals by quadrilateral isoparametric elements in detail, one may refer to [3], here it is introduced briefly.

Figure 6: A quarter-point singular quadrilateral element

As shown in Fig. 7, there are four quadrilateral elements A, B, C, D. In the computation of the double layer surface (Sx & Γξ) integrals, two elements will form a pair. One appears in the Sx, while the other appears in Γξ. There exist four kinds of cases: coincident elements, e.g., Ax&Aξ; adjacent elements sharing one edge, e.g., Ax&Bξ sharing edge pq; adjacent elements sharing one vertex, e.g., Ax&Cξ sharing vertex p; distinct elements, e.g., Ax&Dξ. Numerical integral for a pair of distinct elements do not need special treatment. But for the first three cases, a coordinate transformation is used for numerical integration, which can introduce a Jacobian exploited to cancel singularity of the boundary integral.

Figure 7: Cases of pairs of quadrilateral elements

For a pair of distinct elements, standard isoparametric coordinate transformation is used together with the standard Gauss-Legendre quadrature. As an example, the double layer surface integral considering gravitational loads obtained by the divergence theorem in Eq. (49) is considered at here.

I=∫Stδub(x)na(x)dSx∫Su+Stρgj18π(1−ν)1r[(2ν−1)δabnj(ξ)+2(1−ν)δajnb(ξ)+2(1−ν)δbjna(ξ)+r,ar,bnj(ξ)]dΓξ. (51)

For simplicity, we rewrite it as

I=∫SdSx∫SB(x,ξ)dΓξ=∫01∫01∫01∫01B′[x(x1′,x2′),ξ(ξ1′,ξ2′)]dx1′dx2′dξ1′dξ2′ (52)

where x1′,x2′,ξ1′,ξ2′ are isoparametric coordinates corresponding to Cartesian coordinates x1,x2,ξ1,ξ2. It should be noted that, in Eq. (56), B′[x(x1′,x2′),ξ(ξ1′,ξ2′)] includes the Jacobians of the coordinate transformation.

For cases of coincident elements, adjacent elements sharing one edge, adjacent elements sharing one vertex, further coordinate transformations are given in below to cancel the singularity caused by 1/r appearing in Eq. (56).

For a pair of coincident elements, local isoparametric coordinates are shown in Fig. 8. The boundary integral domain is partitioned into 8 subdomains. For each case we may implement a further transformation of variables listed in Table 1.

Figure 8: Isoparametric coordinates for a pair of coincident elements

In Table 1, v1,v2,v3,v4 are defined as follows:

v1=w1v2=w1w2v3=w3(1−w1)v4=w4(1−w1w2) with 0≤w1≤10≤w2≤10≤w3≤10≤w4≤1. (53)

The Jacobian for such a variable transformation can be used to cancel the singularity in Eq. (56):

J=w1(1−w1)(1−w1w2). (54)

For a pair of coincident elements, Eq. (56) can be rewritten as

I=∑case=18∫01∫01∫01∫01B′[x(x1′,x2′),ξ(ξ1′,ξ2′)]w1(1−w1)(1−w1w2)dw1dw2dw3dw4. (55)

For a pair of common-edge elements, local isoparametric coordinates are shown in Fig. 9.

Figure 9: Local isoparametric coordinates for a pair of common-edge elements

This boundary integral domain is partitioned into 6 subdomains. For each case we may implement a transformation of variables listed in Table 2.

In Table 2, v1,v2,v3,v4,v5 and J1,J2 are defined as follows:

v1=w1v2=w1w2v3=w1w3v4=w4(1−w1)v5=w4(1−w1w2) with 0≤w1≤10≤w2≤10≤w3≤10≤w4≤1. (56)

Jacobians of the variable transformation are

J1=w12(1−w1)J2=w12(1−w1w2). (57)

For a pair of elements with a common vertex, local isoparametric coordinates is shown in Fig. 10.

Figure 10: Isoparametric coordinates for a pair of elements with a common vertex

This boundary integral domain is partitioned into 4 subdomains. For each case, a transformation of variables listed in Table 3 is implemented.

Variables v1,v2,v3,v4 are defined as follows:

v1=w1v2=w1w2v3=w1w3v4=w1w4with0≤w1≤10≤w2≤10≤w3≤10≤w4≤1. (58)

The Jacobian of the variable transformation can be used to cancel the singularity in Eq. (56):

J3=w13. (59)

3  Numerical Examples without Cracks

In this section and the next section some examples without or with crack are implemented respectively to verify SGBEM-div or SGBEM-RIM developed in Section 2.

3.1 Numerical Test of the Effect of the Number of Integration Points

In this section, the double surface integral term in Eq. (55), for a pair of coincident square elements, is evaluated using the quadrature method given in Section 2.3.3, considering the problem of a cube of two kinds of meshes undergoing gravity given in Section 3.2. Fig. 11 shows the logarithmic value of the absolute value of relative errors for the numerical integration of both a pair of square elements and a pair of distorted elements. The error is very small when the number of Gauss integration points is larger than 6. Thus, 8 gauss points are used for the evaluation of double layer surface integrals in the following examples except for the cube undergoing gravitational loads in Section 3.2.

The effect of the number of integration points is shown in Fig. 12, where the relative error is defined as follows: relative error=[I(n)−I(48)]/I(48)×100%, where I(n) is evaluated double surface integral with n Gauss points.

Figure 11: Relative errors for the evaluated weakly-singular boundary integral

3.2 A Cube Undergoing Gravitational Loads

We consider a cube with dimensions of 10 mm×10 mm×10 mm [13], which is discretized into 96 quadratic boundary elements with 290 boundary nodes. Two kinds of meshes of the cube is presented in Fig. 12. The surface z=0 is completely fixed. The elastic constants are chosen to be the Young’s modulus E=1000Mpa and the Poisson’s ratio ν=0.

Figure 12: Mesh of a cube (a) elements being square, (b) elements being distorted

The gravitational force ρg3=−10 Mpa/mm is considered. And the analytical solution for the vertical displacement is

uz=ρg3Ez(L−z2). (60)

Because the analytical solution is only quadratic with respect to z coordinate, 3 Gauss points are used for the evaluation of vertical displacements along the direction z shown in Table 4. The computational results of both square and distorted elements are in excellent agreement with the exact solution.

3.3 A Rotating Disk

In the second example, a disk with inner radius of 0.1m and outer radius of 0.2 m, rotating at a constant angular speed ω=10000rpm (Fig. 13), is considered. The thickness of this disk is t=0.02m. The elastic constants are chosen to be the Young’s modulus E=7000Mpa and the Poisson’s ratio ν=0.3; density ρ=2800 kg/m3. All the boundary surface of this disk is free from traction. The distribution of displacement in a rotating elastic disk

uR=3+ν8ρω2[(Ri2+Ro2)(1−ν)+Ri2Ro2(1+ν)1R2−1−ν23+νR2]RE (61)

can be found in [23] where R is the radial coordinate, E the Young’s modulus and ν the Poisson’s ratio. The boundary of the disk is discretized with 3 elements in radial direction, 16 elements in circumferential direction, and 1 element in axial direction (Fig. 14). 4 nodes on the x-y plane highlighted in Fig. 14 are fixed in z direction; 2 nodes on the x-z plane are fixed in y direction; and 2 nodes on the y-z plane are fixed in x direction.

Figure 13: A rotating disk

Figure 14: SGBEM dense mesh of the rotating disk

Table 5 shows the computed radial displacements with the mesh shown in Fig. 14. “Exact” denotes exact solutions by the Eq. (61). For each point, “Maximum error” of SGBEM-div and SGBEM-RIM is computed with the exact solution as the reference. As can be seen, computational results by SGBEM-div and SGBEM-RIM are in excellent agreement with the exact solutions.

4  Numerical Examples with Cracks

In this section, numerical examples with cracked solids considering body forces are given. In each example, after obtaining the displacement discontinuities for the quarter-point node using the developed SGBEM method, displacement extrapolation is used to calculate the stress intensity factors.

4.1 A Cuboid Hanging under Its Own Weight with a Through-Thickness Crack

Consider a solid cuboid with a crack of length 2a (see Fig. 15) under gravitational loads [24], where l = 4, b = 1, h = 0.5l, a = 0.1, t = 0.2, ρg = −10. The elastic constants are chosen to be E = 1000 and ν = 0.

Figure 15: A cracked cuboid hanging under its own weight

Computed stress intensity factors are presented in Table 6, in which “Error” means the relative error between SGBEM--div and FEM solution. For this through-thickness crack, KI results computed by SGBEM-RIM are in better agreement with the FEM solution.

4.2 A Rotating Disk with a through-Thickness Crack

A rotating disk with a through-thickness crack (a = 0.03 m) is computed shown in Fig. 16. The rotating disk is identical to the disk in Section 3.3. Again, excellent agreement between the computed SGBEM results and FEM results are shown in Tables 7 and 8.

Figure 16: SGBEM mesh of a cracked rotating disk

4.3 A Rotating Disk with Semi-Elliptic Surface Cracks

This section gives a series of results for a cracked disk in Fig. 17 with various semi-elliptic surface cracks, shown in Fig. 18. All the parameters of this disk are identical to that of disk in Section 3, except for the semi-elliptic cracks. Various semi-elliptic cracks with a fixed depth (a = 0.004 m), and various semi-elliptic cracks with a fixed length/depth ratio (b/a = 2), are computed using both SGBEM-div and SGBEM-RIM.

Figure 17: A disk with a semi-elliptic surface crack

Figure 18: Various semi-elliptic cracks with a fixed depth (a = 0.004 m), and various semi-elliptic cracks with a fixed length/depth ratio (b/a = 2)

For simplicity, we give the stress intensity factor KI at point P, i.e., the deepest point of various semi-elliptic cracks, as shown in Figs. 19, and 20. These results can be used for the benchmark solutions for future studies.

Figure 19: KI at the deepest point of semi-elliptic cracks with a fixed depth (a = 0.004 m)

Figure 20: KI at the deepest point of semi-elliptic cracks with a fixed length/depth ratio (b/a = 2)

5  Conclusions

In this paper, weakly-singular SGBEM for fracture analysis of three-dimensional structures considering rotational inertia and gravitational forces is developed. By using the divergence theorem (div) or the radial integration method (RIM), rotational inertia or gravitational forces induced domain integrals are transformed into boundary integrals. The derived boundary integral terms with the gravitational and inertial forces are weakly-singular, which only influence the SGBEM right-hand-side vector.

Several numerical examples of solids with and without cracks undergoing body forces are studied. The calculated stress intensity factors and displacements show high accuracy compared with reference solutions. The test of numerical integration also shows that only a small number of quadrature points are needed.

The symmetric Galerkin boundary element method considering gravity and inertia loads presented in this paper appears promising in the fracture analysis of structural components with body forces, such as dams and rotating machineries. Furthermore, with some effort, the methodology given in this study can also be extended to deal with domain integrals for SGBEM with thermoelastic problems, which will be given in a subsequent work.

Funding Statement: The first four authors acknowledge the support of the National Natural Science Foundation of China (12072011).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Sutradhar, A., Paulino, G. H., Gray, L. J. (2008). Symmetric galerkin boundary element method. Berlin, Heidelberg: Springer.

 2.  Bonnet, M., Maier, G., Polizzotto, C. (1998). Symmetric galerkin boundary element methods. Applied Mechanics Reviews, 51(11), 669–704. DOI 10.1115/1.3098983.

 3.  Novati, G., Frangi, A. (2002). Symmetric galerkin BEM in 3D elasticity: Computational aspects and applications to fracture mechanics. In: Selected topics in boundary integral formulations for solids and fluids. Vienna: Springer.

 4.  Bonnet, M. (1999). Boundary integral equation methods for solids and fluids. Meccanica, 34(4), 301–302. DOI 10.1023/A:1004795120236.

 5.  Li, S., Mear, M. E. (1998). Singularity-reduced integral equations for displacement discontinuities in three-dimensional linear elastic media. International Journal of Fracture, 93(1), 87–114. DOI 10.1023/A:1007513307368.

 6.  Okada, H., Rajiyah, H., Atluri, S. N. (1988). A novel displacement gradient boundary element method for elastic stress analysis with high accuracy. Journal of Applied Mechanics, 55(4), 786–794. DOI 10.1115/1.3173723.

 7.  Okada, H., Rajiyah, H., Atluri, S. N. (1989). Non-hyper-singular integral-representations for velocity (displacement) gradients in elastic/plastic solids (small or finite deformations). Computational Mechanics, 4(3), 165–175. DOI 10.1007/BF00296664.

 8.  Han, Z. D., Atluri, S. N. (2003). On simple formulations of weakly-singular traction & displacement BIE, and their solutions through Petrov-Galerkin approaches. Computer Modeling in Engineering & Sciences, 4(1), 5–20. DOI 10.1.1.610.8720.

 9.  Dong, L. T., Atluri, S. N. (2012). SGBEM (Using non-hyper-singular traction BIEand super elements, for non-collinear fatigue-growth analyses of cracks in stiffened panels with composite-patch repairs. Computer Modeling in Engineering & Sciences, 89(5), 417–458. DOI 10.3970/cmes.2012.089.417.

10. Cruse, T. A. (1975). Boundary-integral equation method for three-dimensional elastic fracture mechanics analysis. AFOSR-TR-75-0813 Interim Report.

11. Cruse, T. A., Snow, D. W., Wilson, R. B. (1977). Numerical solutions in axisymmetric elasticity. Computers & Structures, 7(3), 445–451. DOI 10.1016/0045-7949(77)90081-5.

12. Danson, D. J. (1981). A boundary element formulation of problems in linear isotropic elasticity with body forces. In: Boundary element methods. Berlin, Heidelberg: Springer.

13. Gao, X. W. (2002). The radial integration method for evaluation of domain integrals with boundary-only discretization. Engineering Analysis with Boundary Elements, 26(10), 905–916. DOI 10.1016/S0955-7997(02)00039-5.

14. Nardini, D., Brebbia, C. A. (1983). A new approach to free vibration analysis using boundary elements. Applied Mathematical Modelling, 7(3), 157–162. DOI 10.1016/0307-904X(83)90003-3.

15. Partridge, P. W., Brebbia, C. A., Wrobel, L. C. (1992). The dual reciprocity boundary element method. Southampton Boston: Computational Mechanics Publications.

16. Nowak, A. J., Brebbia, C. A. (1989). The multiple-reciprocity method. A new approach for transforming BEM domain integrals to the boundary. Engineering Analysis with Boundary Elements, 6(3), 164–167. DOI 10.1016/0955-7997(89)90032-5.

17. Aliabadi, F. M. H. (2018). Boundary element methods. In: Encyclopedia of continuum mechanics. Berlin, Heidelberg: Springer.

18. Brebbia, C. A. (2017). The birth of the boundary element method from conception to application. Engineering Analysis with Boundary Elements, 77(4), iii–x. DOI 10.1016/j.enganabound.2016.12.001.

19. Brebbia, C. A. (1983). Progress in boundary element methods, vol. 2. New York: Springer.

20. Frangi, A., Novati, G., Springhetti, R., Rovizzi, M. (2002). 3D fracture analysis by the symmetric Galerkin BEM. Computational Mechanics, 28(3–4), 220–232. DOI 10.1007/s00466-001-0283-x.

21. Fung, Y., Tong, P., Chen, X. (2017). Classical and computational solid mechanics. New Jersey: World Scientific.

22. de Klerk, J. H. (2011). Building a body of knowledge: Cauchy principal value and hypersingular integrals. AIP Conference Proceedings, 1389(1), 456–459. DOI 10.1063/1.3636762.

23. Timoshenko, S. P., Goodier, J. N. (1970). Theory of elasticity. New York: McGraw-Hill.

24. Ostanin, I. A., Mogilevskaya, S. G., Labuz, J. F., Napier, J. (2011). Complex variables boundary element method for elasticity problems with constant body force. Engineering Analysis with Boundary Elements, 35(4), 623–630. DOI 10.1016/j.enganabound.2010.11.008.

Appendix

Kernel functions listed here are utilized in the numerical implementation of the SGBEM. Kernel functions (A1)–(A3) appear in the displacement boundary integral equation; kernel functions (A2)–(A4) appear in the traction boundary integral equation. μ is the shear modulus.

ui∗p(x,ξ)=116πμ(1−υ)r[(3−4υ)δip+r,ir,p] (A1)

Gij∗p(x,ξ)=18π(1−υ)r[(1−2υ)eipj+eikjr,kr,p] (A2)

φij∗p(x,ξ)=δpj14πr2r,i (A3)

Hijpq∗(x,ξ)=μ8π(1−υ)r[4υδiqδjp−δipδjq−2υδijδpq+δijr,pr,q+δpqr,ir,j−2δipr,jr,q−δjqr,ir,p] (A4)

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.