Numerical integration to obtain moment of inertia of nonhomogeneous material

Engineering Analysis with Boundary Elements 101 (2019) 149–155

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Numerical integration to obtain moment of inertia of nonhomogeneous material Yoshihiro Ochiai Department of Mechanical Engineering, Kindai University, Kowakae 3-4-1, Higashi-Osaka-shi, 577-8502 Japan

a r t i c l e

i n f o

Keywords: Moment of inertia Boundary element method Interpolation Integration Triple-Reciprocity Method

a b s t r a c t The moment of inertia of a continuous object with an arbitrary shape made of a nonhomogeneous material is usually calculated by dividing it into small domains. However, it is a burdensome process to specify the density of the small domains. When the Monte Carlo method is used in the case of an arbitrary shape, the computation time increases. In this paper, a technique of easily calculating the moment of inertia of a 3D nonhomogeneous material using boundary integral equations is proposed. It is also shown how to calculate the mass, primary moment, and center of mass of an arbitrary object made of a nonhomogeneous material. A technique employed in the triplereciprocity boundary element method is used to evaluate integral. In this paper, a formulization of the boundary element method is utilized, and a technique for the direct numerical integration of the three-dimensional domain using a three-dimensional interpolation method without carrying out domain division is proposed. To investigate the efficiency of this technique, several numerical examples are given.

1. Introduction Research on rotating machine elements made of functionally graded materials has been reported [1,2]. The moment of inertia of an object with an arbitrary shape made of a nonhomogeneous material is usually calculated by dividing it into small domains. Apart from the Monte Carlo method, there is no good numerical integration method for an object with an arbitrary shape. If a domain is divided into cubes or a tetrahedra, refined numerical integration methods can be used [3]. These considerations resemble a comparison between the finite element method and the boundary element method. We can effectively use an automatic element dividing method in the finite element method for decomposition into elementary regions [4]. However, it is a burdensome process to specify the density of the small domains of a nonhomogeneous material. When the Monte Carlo method is used in the case of an arbitrary shape, the computation time increases. Gao proposed the radial integral method (RIM) without internal cells [5], and Mohammadi et al. proposed the Cartesian transformation (CTM) for the boundary element method [6]. However, these methods are not suitable for calculating the moment of inertia of an object made of functionally graded materials. Another application of this research is the estimation of the moment of inertia of human and animal bodies made of nonuniform materials [7]. In this paper, a technique for easily calculating the moment of inertia of a three-dimensional nonhomogeneous material using boundary integral equations is proposed. It is also shown how to calculate the mass, primary moment, and center of mass of an arbitrary object made

of a nonhomogeneous material. A technique employed in the boundary element method is used to evaluate the integral. In this paper, a formulization of the boundary element method is utilized [8], and a technique for the direct numerical integration of a three-dimensional domain without carrying out domain division is proposed. The numerical integration is performed on a polyharmonic function using a previously reported interpolation method [9–11]. Moreover, this numerical integration has also been generalized to the meshless boundary element method [12–15]. In the numerical integral calculation, the numerical integration of an arbitrary shape is possible, and three-dimensional integration is approximately changed into two-dimensional integration using Green’s theorem, which is possible even in the presence of a singularity. In the numerical integration, surface elements and internal points are used. This technique utilizes the same concept as that in the triple-reciprocity boundary element method [14]. After the given integral function is interpolated using boundary integral equations, the value of the numerical integration is obtained. As demonstrated in this paper, the numerical integral formula of the three-dimensional problem is applied in the same manner as various boundary element methods. 2. Theory 2.1. Moment of inertia The density distribution of a nonhomogeneous body in a domain Ω is denoted as 𝜌(q), where q is a point inside the body. The mass M of the

E-mail address: [email protected] https://doi.org/10.1016/j.enganabound.2019.01.001 Received 15 May 2018; Received in revised form 4 December 2018; Accepted 9 January 2019 Available online 16 January 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

Y. Ochiai

Engineering Analysis with Boundary Elements 101 (2019) 149–155

body is given by 𝑀=

∫Ω

𝜌(𝑞)𝑑 Ω.

(1)

The moment of inertia Ix (p) about the x-axis of a nonhomogeneous body is given by [16] 𝐼𝑥 (𝑝) =

𝜌(𝑞)(𝑦2 + 𝑧2 )𝑑 Ω.

∫Ω

(2)

The point p is on the x-axis. The moments of inertia about the y- and z-axes are given similarly, and the moment of inertia Iij about the xi and xj -axes is defined as 𝐼𝑖𝑗 (𝑝) =

∫Ω

𝜌(𝑞)𝑥𝑖 𝑥𝑗 𝑑 Ω.

(3)

The polar moment of inertia IPO (p) at zero starting point O(p) is defined as 𝐼𝑃 𝑂 (𝑝) =

∫Ω

𝜌(𝑞)𝑟2 𝑑 Ω.

Fig. 1. Polyharmonic function with spherical volume distribution.

(4)

where r is the distance from O. The primary moment Ji (p) of a nonhomogeneous material about the xi -axis is given by 𝐽𝑖 (𝑝) =

∫Ω

𝜌(𝑞)𝑥𝑖 𝑑 Ω.

On the other hand, the quantity W1 (P) is obtained as follows using the Gauss divergence theorem twice and Eqs. (7) and (8):

(5) 𝑐 𝑊1 (𝑃 ) = −

−

Each domain integral appearing in Eqs. (1)–(5) is expressed by [12] ∫Ω

𝑊1 (𝑞 )𝐹1 (𝑝, 𝑞 )𝑑 Ω,

(6)

∇ 𝑊1 (𝑞) = −𝑊2 (𝑞), 2

𝑀 ∑ 𝑚=1

𝑐 𝑊2 (𝑃 ) =

𝑀 ∑ 𝑚=1

𝑊3𝑃𝐴 (𝑞𝑚 )𝛿(𝑞 − 𝑞𝑚 ).

(8)

𝑇1𝐴 (𝑝, 𝑞) =

𝐴3 3𝑟

(9)

𝐻1𝑖𝑗 = 𝐺1𝑖𝑗 = 𝐻2𝑖𝑗 = 𝐺2𝑖𝑗 =

𝑟< =𝐴 𝑟>𝐴

∫0 ∫0

∫Γ

{ 𝑇 1 (𝑃 , 𝑄)

𝑀 ∑ 𝜕 𝑊2 (𝑄) 𝜕 𝑇 1 (𝑃 , 𝑄) − 𝑊2 (𝑄)}𝑑Γ + 𝑇1𝐴 (𝑃 , 𝑞)𝑊3𝑃𝐴(𝑚) . 𝜕𝑛 𝜕𝑛 𝑚=1

𝜋

∫0

𝑇𝑓 (𝑃 , 𝑄)𝑎2 sin 𝜃𝑑 𝜃 𝑑 𝜙 𝑑 𝑎

(18)

𝜕 𝑇1 (𝑃 , 𝑄) 1 𝛿 + 𝑑 Γ𝑗 2 𝑖𝑗 ∫Γ𝑗 𝜕𝑛 ∫Γ𝑗 ∫Γ𝑗 ∫Γ𝑗

(19)

𝑇 1 ( 𝑃 , 𝑄 ) 𝑑 Γ𝑗

(20)

𝜕 𝑇2 (𝑃 , 𝑄) 𝑑 Γ𝑗 𝜕𝑛

(21)

𝑇 2 ( 𝑃 , 𝑄 ) 𝑑 Γ𝑗

(22) (23)

where the boundary element Γj involves the node Qj . The upper index I in Eq. (20) is related to W3 . W2 is obtained from Eq. (17) as

(13)

𝐇1 𝐖2 = 𝐆1 𝐕2 + 𝐆𝐼1 𝐖3 , where

𝑟4 − 15𝐴4 − 10𝑟2 𝐴2 𝑟< (14) =𝐴 120 TfA , which expresses the state of the uniformly distributed polyharmonic function Tf (P,Q) in a spherical region with radius A, as shown in Fig. 1, is introduced for a smooth interpolation [12]. 2𝜋

(16)

(12)

𝑇2𝐴 (𝑝, 𝑞) = −

𝐴

𝑇2𝐴 (𝑃 , 𝑞)𝑊3𝑃𝐴(𝑚) .

where H1 , G1 , H2 , G2 , and are the matrices with the following elements for a given boundary point i [12]:

(11)

𝐴3 2 𝐴2 (𝑟 + ) 6𝑟 5

𝑇𝑓 𝐴 (𝑃 , 𝑄) =

𝑚=1

𝐺2𝐼𝑖𝑗 = 𝑇2𝐴 (𝑃 , 𝑞 𝐼 ),

3𝐴2 − 𝑟2 𝑇1𝐴 (𝑝, 𝑞) = 6 𝑇2𝐴 (𝑝, 𝑞) =

𝑀 ∑

𝐆𝐼2

(10)

𝑟>𝐴

𝑊𝑓 (𝑄)}𝑑Γ −

𝜕𝑛

𝐇1 𝐖1 = 𝐆1 𝐕1 + 𝐇2 𝐖2 − 𝐆2 𝐕2 − 𝐆𝐼2 𝐖3 ,

W2 (q) = 0 is assumed on the boundary [12]. Eqs. (7) and (8) are solved using boundary integral equations. Then, a polyharmonic function Tf (P,Q) and a polyharmonic function with volume distribution TfA are introduced. 𝑟2𝑓 −3 𝑇𝑓 (𝑃 , 𝑄) = 4𝜋(2𝑓 − 2)!

𝜕𝑛

𝜕 𝑊𝑓 (𝑄)

In practice, W1 in Eq. (16) is given; however, W2 , 𝜕 Wf /𝜕 n, and W3 must be obtained numerically. Constant elements are used for the boundary. The upper index P corresponds to W3 . Replacing Wf and 𝜕Wf / 𝜕n by the vectors Wf and Vf , respectively, and discretizing Eq. (17), we obtain [11,12]

where 𝛿(q − qm ) is the Dirac delta function. From Eqs. (7) and (8), the following equation can be obtained: ∇4 𝑊1 (𝑞) =

𝜕 𝑇𝑓 (𝑃 , 𝑄)

{𝑇𝑓 (𝑃 , 𝑄)

(17)

(7)

𝑊3𝑃𝐴 (𝑞𝑚 )𝛿(𝑝 − 𝑞𝑚 ),

∫Γ

Similarly, the following equation can be obtained from Eq. (8):

where 𝑊1𝑆 (𝑞) is an arbitrary distribution function, which is the density 𝜌(q) in this paper. The function F1 (p,q) is (y2 + z2 ), r2 , and so forth in Eqs. (1)–(5), where r is the distance between the observation point p and the source point q. The following equations can be used for the three-dimensional interpolation [12,17,18]:

∇2 𝑊2 (𝑞) = −

(−1)𝑓

𝑓 =1

2.2. Interpolation

𝐼(𝑝) =

2 ∑

𝐆𝐼1

(24)

is a matrix with the following elements:

𝐺1𝐼𝑖𝑗 = 𝑇1𝐴 (𝑃 , 𝑞 𝐼 ).

(25)

In the same manner, from Eq. (16), using the vector notation 𝐖𝐼1 for the function values W1 (pi ) at internal points with the value of the function W(pI ), we obtain [11,12] 𝐖(𝑝𝐼 ) = −𝐇3 𝐖1 + 𝐆3 𝐕1 + 𝐇4 𝐖2 − 𝐆4 𝐕2 − 𝐆𝐼3 𝐖3 ,

(15) 150

(26)

Y. Ochiai

Engineering Analysis with Boundary Elements 101 (2019) 149–155

where H3 , G3 , H4 , G4 and 𝐆𝐼3 are the matrices with the following elements: 𝐻3𝑖𝑗 = 𝐺3𝑖𝑗 = 𝐻4𝑖𝑗 = 𝐺4𝑖𝑗 =

∫Γ𝑗 ∫Γ𝑗 ∫Γ𝑗 ∫Γ𝑗

Differentiating Eqs. (38) and (39) with respect to the normal vector n on a boundary, we obtain

𝜕 𝑇1 (𝑝𝐼 , 𝑄) 𝑑 Γ𝑗 . 𝜕𝑛

(27)

𝜕 𝐹2 𝑟1−𝑘 𝜕𝑟 = 𝜕𝑛 (3 − 𝑘) 𝜕𝑛

(40)

𝑇 1 ( 𝑝 𝐼 , 𝑄 ) 𝑑 Γ𝑗

(28)

𝜕 𝐹3 𝑟3−𝑘 𝜕𝑟 = . 𝜕𝑛 (2 − 𝑘)(3 − 𝑘)(5 − 𝑘) 𝜕𝑛

(41)

𝜕 𝑇2 (𝑝𝐼 , 𝑄) 𝑑 Γ𝑗 𝜕𝑛

(29)

𝑇 2 ( 𝑝 𝐼 , 𝑄 ) 𝑑 Γ𝑗

(30)

Substituting Eq. (39) into Eq. (34), we obtain 𝐹3𝐴 =

+{(7 − 𝑘)𝐴 + 𝑟}(−𝐴 + 𝑟)7−𝑘

(31) 𝐹3𝐴 =

Assuming W2 = 0 to interpolate the distribution, the following equation is obtained from Eqs. (18), (24) and (26): −𝐆𝐼2 ⎤⎧ 𝐕1 ⎫ ⎪ ⎪ 𝐆𝐼1 ⎥⎨ 𝐕2 ⎬ ⎥⎪ ⎪ 𝐼 − 𝐆 3 ⎦ ⎩𝐖 3 ⎭

−𝐆2 𝐆1 −𝐆4

⎧ ⎫ 𝐇1 𝐖1 ⎪ ⎪ =⎨ 𝟎 ⎬. ⎪𝐇3 𝐖1 + 𝐖(𝑝𝐼 )⎪ ⎩ ⎭

(32) ∫Ω

(33)

The function FfA is introduced in the same manner as Eq. (15) as follows: 𝐹𝑓 𝐴 (𝑃 , 𝑄) =

2𝜋

∫0 ∫0

𝜋

∫0

𝐹𝑓 (𝑃 , 𝑄)𝑎2 sin 𝜃𝑑 𝜃 𝑑 𝜙 𝑑 𝑎.

2 ∑

(−1)𝑓

𝑓 =1

−

𝜕 𝐹𝑓 +1 (𝑝, 𝑄) 𝜕𝑛

∫Γ

{𝐹𝑓 +1 (𝑝, 𝑄)

𝑊𝑓 (𝑄)}𝑑Γ +

𝑀 ∑ 𝑚=1

𝜕𝑛 (35)

1 [ 𝑟2 𝐹𝑓 (𝑝, 𝑞)𝑑 𝑟]𝑑 𝑟. ∫ 𝑟2 ∫

1 𝑟𝑘

(𝑘 ≠ 2, 3, 4, 5),

𝐹3 (𝑝, 𝑞) =

1 , (2 − 𝑘)(3 − 𝑘)(4 − 𝑘)(5 − 𝑘)𝑟𝑘−4

(39)

𝑟2 6

(45)

𝐹3 (𝑝, 𝑞) =

𝑟4 120

(46)

𝜕 𝐹2 𝑟 𝜕𝑟 = 𝜕𝑛 3 𝜕𝑛

(47)

𝜕 𝐹3 𝑟3 𝜕𝑟 = 𝜕𝑛 30 𝜕𝑛

(48)

𝜋𝐴3 (3𝐴4 + 14𝐴2 𝑟2 + 7𝑟4 ). 630

(49)

𝜋 𝑎3 ⟨ 2𝛿𝑖𝑗 (3520𝑎6 − 19008𝑎2 𝑟4 − 7040𝑟6 ) 39916800 +55𝑟,𝑖 𝑟,𝑗 (1112𝑎6 + 8781𝑎4 𝑟2 + 6912𝑎2 𝑟4 +

⟩ 1920𝑟6 )

𝑟 ≤ 𝑎.

1 𝑟,𝑖 (𝑘 ≠ 1, 3, 4, 6 real number ). (52) 𝑟𝑘 The following high-order functions are respectively obtained by the partial differentiation of Eqs. (38)–(43): 1 𝑟,𝑖 (1 − 𝑘)(4 − 𝑘)𝑟𝑘−2 [ ] 𝑛𝑖 𝜕 𝐹2 1 𝜕𝑟 = + 𝑟,𝑖 𝑘 −1 𝜕𝑛 (1 − 𝑘) 𝜕𝑛 (4 − 𝑘)𝑟

𝐹2 (𝑝, 𝑞) =

where k is a real number. Substituting Eq. (37) into Eq. (36), we obtain (38)

(44)

𝐹1 (𝑝, 𝑞) =

(37)

1 (2 − 𝑘)(3 − 𝑘)𝑟𝑘−2

𝜌(𝑞)𝑑 Ω = 𝑀,

In the next case, the following fundamental solution has the partial differential term r, i = 𝜕 r/𝜕 xi = xi /r:

(36)

𝐹2 (𝑝, 𝑞) =

∫Ω

(51)

First, the following function is considered: 𝐹1 (𝑝, 𝑞) =

(43)

𝐹2 (𝑝, 𝑞) =

𝐸3𝐴 =

Thus, the domain integral is replaced by boundary integrals and a sum, and the domain discretization procedure is completely eliminated. Using Eq. (33), the function Ff +1 in three dimensions is obtained as 𝐹𝑓 +1 (𝑝, 𝑞) =

𝑟 ≤ 𝐴.

(42)

Eqs. (42) and (43) become Eq. (49). Even if the position of the observation point p changes, the integration value M does not change. Setting k = −2 in Eq. (37) and using F1 = r2 , the polar moment of inertia IPO in Eq. (4) can be calculated, and Eqs. (42) and (43) respectively become 𝜋 𝑎3 ⟨ 𝐸3𝐴 = 19958400 −𝛿𝑖𝑗 (14080𝑎6 + 59345𝑎4 𝑟2 + 19008𝑎2 𝑟4 + 7040𝑟6 ) (50) ⟩ +880𝑟,𝑖 𝑟,𝑗 𝑟2 (133𝑎4 + 266𝑎2 𝑟2 + 60𝑟4 ) 𝑟>𝑎

(34)

𝜕 𝑊𝑓 (𝑄)

𝐹3𝐴 (𝑝, 𝑞 )𝑊3𝐼𝐴(𝑚) (𝑞 ).

𝑊1 (𝑞)𝑑 Ω =

𝐹3𝐴 (𝑝, 𝑞) =

Using Eq. (33) and the Gauss divergence theorem, the following equation is obtained because the function W1 satisfies Eqs. (7) and (8): 𝐹1 (𝑝, 𝑞 )𝑊1 (𝑞 )𝑑 Ω =

1 𝓁−𝑘

𝑟>𝐴

which enables us to obtain the mass M in Eq. (1). In this case, the integral value M is independent of the observation point p, which can be out of the domain [14]. From Eqs. (38)–(43), we obtain

Numerical integration in an arbitrary domain is presented. After interpolation using Eqs. (16) and (17), numerical integration in domain Ω is performed using the interpolated values. A new function Ff +1 is introduced as follows: ∇2 𝐹𝑓 +1 = 𝐹𝑓 .

1 𝓁−𝑘

If k = 0, Eq. (35) becomes

2.3. Numerical integration

𝐴

2𝜋 [ {(7 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)7−𝑘 𝑟 8 ]∏ −{(7 − 𝑘)𝐴 + 𝑟}(𝐴 − 𝑟)7−𝑘 𝓁=2

From Eq. (32), we can obtain V1 , V1 and W3 . In the same manner, V2 = 0 can be assumed on a symmetric axis to interpolate a symmetric distribution. When using constant elements and dividing the boundary into N0 elements and M internal points, simultaneous linear algebraic equations with 2N0 + M unknowns must be solved.

∫Ω

8 ]∏ 𝓁=2

𝐺3𝐼𝑖𝑗 = 𝑇2𝐴 (𝑝𝐼 , 𝑞 𝐼 ),

⎡𝐆 1 ⎢𝟎 ⎢ ⎣𝐆 3

2𝜋 [ {(7 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)7−𝑘 𝑟

𝐹3 (𝑝, 𝑞) = 151

𝑟4−𝑘 𝑟, (1 − 𝑘)(3 − 𝑘)(4 − 𝑘)(6 − 𝑘) 𝑖

(53)

(54)

(55)

Y. Ochiai

Engineering Analysis with Boundary Elements 101 (2019) 149–155

𝜕 𝐹3 1 𝜕𝑟 = [𝑛𝑖 + (3 − 𝑘)𝑟,𝑖 ] 𝜕𝑛 𝜕𝑛 (6 − 𝑘)(4 − 𝑘)(3 − 𝑘)(1 − 𝑘)𝑟𝑘−3

𝐹3𝐴 =

𝜋(9−𝑘) ⟨ 𝛿𝑖𝑗 [−𝐴{(𝐴 + 𝑟)9−𝑘 − (𝐴 − 𝑟)9−𝑘 } + 𝑟{(9 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)8−𝑘 𝐹3𝐴 = −2 (2−𝑘)𝑘𝑟3 +𝑟{(9 − 𝑘)𝐴 + 𝑟}(𝐴 − 𝑟)8−𝑘 ] + 𝑟,𝑖 𝑟,𝑗 [3𝐴{(𝐴 + 𝑟)9−𝑘 + (𝐴 − 𝑟)9−𝑘 } −𝑟{3(9 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)8−𝑘 − 𝑟{3(9 − 𝑘)𝐴 + 𝑟}(𝐴 − 𝑟)8−𝑘 ⟩ +𝑟2 {(9 − 𝑘)𝐴 − 𝑟}(8 − 𝑘)(𝐴 + 𝑟)7−𝑘 − 𝑟2 {(9 − 𝑘)𝐴 + 𝑟}(8 − 𝑘)(𝐴 − 𝑟)7−𝑘 ] 8 ∏ 1 2𝜋𝛿 + 𝑘𝑟𝑖𝑗 [{(7 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)7−𝑘 − {(7 − 𝑘)𝐴 + 𝑟}(𝐴 − 𝑟)7−𝑘 ] × 𝓁−𝑘+2

(56)

2𝜋(8 − 𝑘)𝑟,𝑖 [ −𝐴{(𝐴 + 𝑟)8−𝑘 + (−𝐴 + 𝑟)8−𝑘 } + 𝑟{(8 − 𝑘)𝐴 − 𝑟} (1 − 𝑘)𝑟2

× (𝐴 + 𝑟)7−𝑘 + 𝑟{(8 − 𝑘)𝐴 + 𝑟}(−𝐴 + 𝑟)7−𝑘

8 ]∏ 𝓁=2

1 𝓁−𝑘+1

𝓁=2 8

𝑟>𝐴

∏

×

𝓁=2

1 𝓁−𝑘

𝑟 ≤ 𝐴.

(67)

(57)

𝐹3𝐴 =

Eqs. (2) and (3) are calculated by setting k = −2 in Eq. (60) and using the relation F1 = r, i r,j r2 = xi xj . In addition, Eqs. (66) and (67) respectively become 𝜋 𝑎3 ⟨ 𝐸3𝐴 = 19958400 −𝛿𝑖𝑗 (14080𝑎6 + 59345𝑎4 𝑟2 + 19008𝑎2 𝑟4 + 7040𝑟6 ) (68) ⟩ 𝑟>𝑎 +880𝑟,𝑖 𝑟,𝑗 𝑅2 (133𝑎4 + 266𝑎2 𝑟2 + 60𝑟4 )

2𝜋(8 − 𝑘)𝑟,𝑖 [ −𝐴{(𝐴 + 𝑟)8−𝑘 − (𝐴 − 𝑟)8−𝑘 } + 𝑟{(8 − 𝑘)𝐴 − 𝑟} (1 − 𝑘)𝑟2

× (𝐴 + 𝑟)7−𝑘 + 𝑟{(8 − 𝑘)𝐴 + 𝑟}(𝐴 − 𝑟)7−𝑘

8 ]∏ 𝓁=2

1 𝓁−𝑘+1

𝑟 ≤ 𝐴,

𝐸3𝐴 =

(58) where ni is the component of the outward unit vector from the surface. Setting k = −1 in Eq. (52) and using the relationship F1 = r, i r = xi , the primary moment Ji in Eq. (5) is calculated. In addition, Eqs. (57) and (58) become 𝐹3𝐴

𝑟, 𝜋 𝐴3 𝑟 = 𝑖 (5𝐴4 + 14𝐴2 𝑟2 + 5𝑟4 ). 1050

1 𝑟,𝑖 𝑟,𝑗 𝑟𝑘

(𝑘 ≠ 0 , 2 , 3 , 4 , 5 , 7

2.4. Monte Carlo method To show the validity of the proposed technique, the values obtained by this technique are compared with those calculated by the Monte Carlo method [19]. Although the Monte Carlo method has a drawback of a long computation time, the moment of inertia of an arbitrarily shaped object can be calculated by this method for comparison with the proposed technique. Moreover, although the Monte Carlo method is inapplicable to numerical integration with a singularity, there are no singularities in the integration in this paper. The volume of the rectangular parallelepiped surrounding domain Ω of the nonhomogeneous body under consideration is set to V, and NT is the number of random points. The number of points in domain Ω of the nonhomogeneous body among the random points is NI . The integration values in Eqs. (1)–(5) can be acquired by the following formula:

(60)

The following function is obtained if the second-order partial differentiation of the Eq. (37) is carried out: −𝑘𝛿𝑖𝑗 𝑘(𝑘 + 2) 𝜕2 𝐹1 = + 𝑟,𝑖 𝑟,𝑗 . 𝜕 𝑥𝑖 𝜕 𝑥𝑗 𝑟𝑘+2 𝑟𝑘+2

(61)

Therefore, the following functions are obtained from the first term of the function by carrying out the second partial differentiation of Eqs. (33) and (60): [ ] 2𝛿𝑖𝑗 1 𝐹2 (𝑝, 𝑞) = − 𝑟,𝑖 𝑟,𝑗 (62) 𝑟𝑘−2 𝑘(5 − 𝑘) (2 − 𝑘)(3 − 𝑘) [ ] 2𝛿𝑖𝑗 𝜕𝑟 𝜕 𝐹2 𝑟1−𝑘 𝜕𝑟 = − 𝑛𝑗 𝑟,𝑖 − 𝑛𝑖 𝑟,𝑗 + 𝑘 𝑟,𝑖 𝑟,𝑗 𝜕𝑛 𝑘(5 − 𝑘) 3 − 𝑘 𝜕𝑛 𝜕𝑛 𝐹3 (𝑝, 𝑞) =

[ ] 4𝛿𝑖𝑗 𝑟4−𝑘 − 𝑟,𝑖 𝑟,𝑗 𝑘(2 − 𝑘)(5 − 𝑘)(7 − 𝑘) (3 − 𝑘)(4 − 𝑘)

𝐼𝑀𝑂𝑁 =

(63)

(64)

To confirm the accuracy of the numerical integral, the moment of inertia of the spherical domain Ω of radius R shown in Fig. 2(a) was obtained. For the density 𝜌0 = 1 and the radius R = 10 , the numerically computed value obtained using the proposed technique was Ix = 0.1672252 × 106 , compared with the exact value of 8𝜋 5 𝑅 = 0.1675516 × 106 . (72) 15 In addition, the computation time was 107 s and linear boundary elements were used in the calculation. When quadrilateral boundary elements were used, the numerically computed value obtained using the proposed technique was Ix = 0.1675486 × 106 . When the spherical domain is homogeneous, the internal points of Fig. 2(b) and interpolation are not necessary, and Eq. (35) becomes

𝜋(9−𝑘) ⟨ 𝛿𝑖𝑗 [−𝐴{(𝐴 + 𝑟)9−𝑘 + (−𝐴 + 𝑟)9−𝑘 } + 𝑟{(9 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)8−𝑘 𝐹3𝐴 = −2 (2−𝑘)𝑘𝑟3 + 𝑟{(9 − 𝑘)𝐴 + 𝑟}(−𝐴 + 𝑟)8−𝑘 ] + 𝑟,𝑖 𝑟,𝑗 [3𝐴{(𝐴 + 𝑟)9−𝑘 + (−𝐴 + 𝑟)9−𝑘 } − 𝑟{3(9 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)8−𝑘 − 𝑟{3(9 − 𝑘)𝐴 + 𝑟}(−𝐴 + 𝑟)8−𝑘 ⟩ + 𝑟2 {(9 − 𝑘)𝐴 − 𝑟}(8 − 𝑘)(𝐴 + 𝑟)7−𝑘 +𝑟2 {(9 − 𝑘)𝐴 + 𝑟}(8 − 𝑘)(−𝐴 + 𝑟)7−𝑘 ] 8 ∏ 1 2𝜋𝛿 + 𝑘𝑟𝑖𝑗 [{(7 − 𝑘)𝐴 − 𝑟}(𝐴 + 𝑟)7−𝑘 + {(7 − 𝑘)𝐴 + 𝑟}(−𝐴 + 𝑟)7−𝑘 ] × 𝓁−𝑘+2

∏

𝓁=2

1 𝓁−𝑘

(70)

3. Numerical examples

The following functions are obtained similarly:

𝓁=2 8

𝑁𝐼 𝑉 ∑ 𝜌(𝑞𝑖 )𝐹 (𝑝, 𝑞𝑖 ). 𝑁𝑇 𝑖=1

However, only when point q is inside the domain of an object, the sum in Eq. (70) is added. Point p is judged inside or outside the domain using Eq. (71). When H(p) is 1, it is inside the domain, and when H(p) is 0, it is outside the domain. 𝜕 𝑇1 (𝑝, 𝑞𝑗 ) ∑ 𝐻(𝑝) = 𝑑 Γ𝑗 (71) ∫ 𝜕𝑛 Γ𝑗 𝑗

[ ] 4𝑛𝑚 𝑟,𝑚 𝛿𝑖𝑗 𝜕 𝑟3−𝑘 𝐹3 = − (2 − 𝑘)𝑛𝑚 𝑟,𝑚 𝑟,𝑖 𝑟𝑗 − 𝑛𝑖 𝑟,𝑗 − 𝑟,𝑖 𝑛𝑗 . 𝜕𝑛 𝑘(2 − 𝑘)(5 − 𝑘)(7 − 𝑘) (3 − 𝑘) (65)

×

𝑟 < 𝑎. (69)

(59)

real number).

⟩ 1920𝑟6 )

Linear boundary elements are used in this paper.

Next, the following function, in which two partial differential terms exist, is considered: 𝐹1 (𝑝, 𝑞) =

𝜋 𝑎3 ⟨ 2𝛿𝑖𝑗 (3520𝑎6 − 19008𝑎2 𝑟4 − 7040𝑟6 ) 39916800 +55𝑟,𝑖 𝑟,𝑗 (1112𝑎6 + 8781𝑎4 𝑟2 + 6912𝑎2 𝑟4 +

𝐼𝑥 =

𝑟>𝐴

(66)

∫Ω 152

𝐹1 (𝑝, 𝑞)𝑑 Ω =

∫Γ

𝜕 𝐹2 (𝑝, 𝑄) 𝑑Γ. 𝜕𝑛

(73)

Y. Ochiai

Engineering Analysis with Boundary Elements 101 (2019) 149–155

Table 1 CPU time for different numbers of calculation points in Monte Carlo method. Number of points

Mass M (Exact 4188.79)

Moment of inertia Ix

CPU time (s)

27,000 125,000 343,000 1,000,000

4147.556 4191.936 4190.134 4187.64

165407.7014 167175.3023 167359.0904 167352.2021

150 641 1744 4225

the number of calculation points is set to 106 , Ix can be calculated with satisfactory accuracy. In the next case, the distribution of density in a nonhomogeneous sphere of radius R is given as 𝜌1 (𝑟) = 𝜌0 [2 − (𝑟∕𝑅)2 ],

(74)

where r is the distance from the center of the sphere and 𝜌0 = 1. The quadrilateral elements shown in Fig. 2(a) and the internal points at random positions shown in Fig. 2(b) were used for the numerical integration. The numerically computed value obtained using the proposed technique was Ix = 0.1072931 × 106 , compared with the exact value of 𝐼𝑥 =

12𝜋 5 𝑅 = 0.1077117 × 106 . 35

(75)

Moreover, the numerically computed value of the polar second moment about the starting point O was IPO = 0.3217785 × 106 , and a value near the exact solution was acquired as follows: 𝐼𝑃 𝑂 =

∫

𝜌1 (𝑟)𝑟2 𝑑𝑉 = 𝜋𝜌0 [

36 5 𝑅 ] = 0.3231352 × 106. 35

(76)

Next, we obtained the moment of inertia of a nonhomogeneous cube of side L = 10 shown in Fig. 3(a) whose density was given by 𝜌1 (𝑥) = 𝜌0 exp(𝛼𝑥)

(77)

where 𝛼 = 0.1. The exact moment of inertia is given by 𝐼𝑦 (𝑥) =

∫

𝜌0 exp(𝛼𝑥)(𝑥2 + 𝑧2 )𝑑Ω = 0.129120 × 106 .

(78)

The moment of inertia about the y-axis was Iy = 0.129157 × 106 , which was obtained using the internal points shown in Fig. 3(b). When the density was given by [ ( ) (𝑦𝜋 ) ( )] 𝑥𝜋 𝑧𝜋 sin sin . 𝜌1 (𝑥, 𝑦, 𝑧) = 𝜌0 1 + sin (79) 𝐿 𝐿 𝐿 The moment of inertia about the y-axis was Iy = 0.8200686 × 105 . In addition, the value calculated by the Monte Carlo method was 0.8208332 × 105 . The moment of inertia of the torus shown in Fig. 4 with the internal radius Ri = 4 and the external radius Re = 10 was next obtained. For a homogeneous case, the numerically computed value obtained using the proposed technique was Iz = 0.36171 × 106 , compared with the exact value of 1 𝐼𝑧 = 𝜋 2 𝑅𝑒 𝑅2𝑖 (4𝑅2𝑒 + 3𝑅2𝑖 ) = 0.3537266217 × 106 (80) 2 When quadrilateral boundary elements were used, the numerically computed value obtained using the proposed technique was Iz = 0.353723 × 106 . When the spherical domain was nonhomogeneous and the density was given by [ ] 𝑧 , 𝜌1 ( 𝑧 ) = 𝜌 0 2 + (81) 𝑅

Fig. 2. Spherical object. (a) Boundary elements (b) Internal points.

the value calculated using the proposed technique was Iz = 0.3752332 × 106 , compared with Iz = 0.3765598 × 106 obtained by the Monte Carlo method. The position (0, 0, z) of the center of gravity was obtained from the primary moment Jz and the mass M using the proposed technique as

When linear boundary elements are used, a curved surface is approximated as many small flat planes and the calculation accuracy is reduced. Table 1 shows the mass M obtained by the Monte Carlo method, the CPU time for different numbers of calculation points, and the calculated moment of inertia Ix . Since there are no singularities in this calculation, if

𝑧= 153

𝐽𝑧 0.3113763 × 104 = = 0.498731. 𝑀 6261.237

(82)

Y. Ochiai

Engineering Analysis with Boundary Elements 101 (2019) 149–155

Fig. 4. Nonhomogeneous torus.

Fig. 3. Cubic region. (a) Boundary elements (b) Internal points.

In the case of the Monte Carlo method, z = 0.5025418. The model of Venus de Milo shown in Fig. 5 was used to verify the validity of the proposed method for an arbitrary shape. The height of the model was 80 mm and the distribution of the density was assumed to be [ ] 𝑧 𝜌1 ( 𝑧 ) = 𝜌 0 1 + (83) . 80

Fig. 5. Venus de Milo.

The calculated mass using the proposed technique was M = 0.1955411 × 106 , compared with M = 0.1928176 × 106 for the Monte Carlo method. Moreover, the calculated moment of inertia was Iy = 0.3752332 × 106 , compared with Iy = 0.3765598817 × 106 for the Monte Carlo method. As an example of an arbitrary density distribution, the points in Venus de Milo shown in Fig. 6(a) were used. The interpolated density distribution is shown in Fig. 6(b).

The calculated moment of inertia using the proposed technique was Iy = 0.9261561 × 108 . Furthermore, a high-order element is required to increase the calculation accuracy. It is not necessary to use internal points of the lattice pattern, and points at arbitrary positions can be used. The proposed method was found to have a shorter calculation time than the Monte Carlo method. 154

Y. Ochiai

Engineering Analysis with Boundary Elements 101 (2019) 149–155

4. Conclusion A simple technique for calculating the moment of inertia of an object made of a nonhomogeneous material by the boundary integral approach was proposed. Since it is not necessary to divide the domain, it is easy to perform the integration for an object with an arbitrary shape and density distribution. Moreover, the calculation time was shorter than that for the Monte Carlo method. Since the integral formula is a general formula, it can be used in many other fields. Also in this technique, although inner points are used, it is easy to specify the values of inner points and it is possible to use inner points at an arbitrary positions in this technique. In the future, research will be carried out on decreasing the computation time and increasing the accuracy through the use of high-order elements. References [1] Rizwan UZK, Sarda A. Finite element analysis of concave thickness FGM rotating disk subject to thermo-mechanical load. Int J Eng Manag Res 2016;6(3):261–5. [2] Mazzei AJJ. On the effect of functionally graded materials on resonances of rotating beams. Shock Vib 2012;19:1315–26. [3] Davis P, Rabinowitz P. Methods of numerical integration. 2nd ed. Academic Press; 1984. p. 341–417. [4] Cheng SW, Dey TK, Shewchuk JR. Delaunay mesh generation. CRC Press; 2013. p. 301–32. [5] Gao XW. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng Anal Bound Elem 2002;26:905–16. [6] Mohammadi M, Hematiyan MR, Aliabadi MH. Boundary element analysis of thermo-elastic problems with non-uniform heat sources. Strain Anal 2010;45:606–27. [7] Krishnan RH, Devanandh V, Brahma AK. Estimation of mass moment of inertia of human body, when bending forward, for the design of a self-transfer robotic facility. J Eng Technol 2016;11(2):166–76. [8] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques −theory and applications in engineering. Springer-Verlag; 1984. p. 46–70. [9] Ochiai Y, Sekiya T. Steady heat conduction analysis by improved multiple-reciprocity boundary element method. Eng Anal Bound Elem 1996;18:111–17. [10] Ochiai Y. Steady heat conduction analysis in orthotropic bodies by triple-reciprocity BEM. Comput Model Eng Sci 2001;2(4):435–45. [11] Ochiai Y. Multidimensional numerical integration for meshless BEM. Eng Anal Bound Elem 2003;27(3):241–9. [12] Ochiai Y, Sladek V. Numerical treatment of domain integrals without internal cells in three-dimensional BIEM formulations. CMES(Comput Model Eng Sci) 2005;6(6):525–36. [13] Ochiai Y. Meshless large plastic deformation analysis considering with a friction coefficient by triple-reciprocity boundary element method. J Comput Mech Exp Measur 2018;6(6):989–99. [14] Ochiai Y. Three-dimensional heat conduction analysis of inhomogeneous materials by triple-reciprocity boundary element method. Eng Anal Bound Elem 2015;51(1):101–8. [15] Ochiai Y, Sekiya T. Generation of free-form surface in CAD for dies. Adv Eng Softw 1995;22:113–18. [16] Juvinall RC, Marshek KM. Fundamentals of machine component design. Wiley; 2006. [17] Micchelli CA. Interpolation of scattered data. Construct Approx 1986;2:12–22. [18] Dyn N. Interpolation of scattered data by radial functions, In: Topics in multivariate approximation, editors. Chui C. K., Schumaker L. L. and Utreras F. I., pp. 47–61, (1987), Academic Press, London. [19] Kroese DP, Taimre T, Botev ZI. Handbook of Monte Carlo methods. Wiley; 2011.

Fig. 6. Arbitrary density distribution in Venus de Milo. (a) Points used to assign density (b) Interpolated density distribution (y = 0).

155