- Email: [email protected]
Fast Multipole BEM for Simulation of 2-D Solids Containing Large Numbers of Cracks
TSINGHUA SCIENCE AND TECHNOLOGY ISSN 1007-0214 11/18 pp76- 81 Volume 10, Number 1, February 2005
Fast Multipole BEM for Simulation of 2-D Solids Containing Large Numbers of Cracks WANG Pengbo (王朋波), YAO Zhenhan (姚振汉) **, WANG Haitao (王海涛) Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China Abstract:
The fast multipole method was used to solve the traction boundary integral equation for 2-D
crack analysis. The use of both multipole and local expansions reduces both the computational complexity and the memory requirement to O(N). The multipole expansion uses a complex Taylor series expansion to reduce the number of multipole moments. The generalized minimum residual method solver (GMRES) was selected as the iterative solver. An improved preconditioner for GMRES was developed which uses less CPU time and less memory. A new initial candidate vector for the iterative solver was developed to further improve the efficiency. The numerical examples apply the method to the analysis of cracks in infinite 2-D space with the largest model having 900 000 degrees of freedom. Key words: fast multipole method; boundary element method (BEM); cracks
Introduction The boundary element method (BEM) is recognized as a powerful method for the study of crack problems because of its semi-analytical nature and boundary-only discretization[1]. Compared with other numerical methods, the reduction in dimensionality dramatically reduces the initial data preparation and remeshing work during crack growth. Additionally, since there are no internal approximations, the singular stress field around the crack tip can be analyzed more accurately and more effectively. However, intrinsic basis for the displacement boundary integral equation (DBIE) causes the equation degeneration when the two surfaces of one crack are considered co-planar. The degeneration caused by the co-planar surfaces can be overcome by using the traction boundary integral equation (TBIE) on one of the two surfaces[2]. The TBIE can be used to develop a single-region formulation for general crack analysis. Received: 2004-08-16; revised: 2004-08-30
﹡﹡ To whom correspondence should be addressed. E-mail: [email protected]; Tel: 86-10-62772913
Conventional BEM is not efficient for large-scale problems because of the dense and asymmetric matrix. The memory requirements and the CPU time have been reduced through the application of the fast multipole method (FMM) to the BEM. Barnes and Hut[3] first used a tree data structure and the concept of multipole expansion to calculate the matrix-vector product without explicitly forming the matrix. The computational complexity and memory requirements for the matrix-vector multiplication are reduced from O(N2) to O(NlogN). Greengard and Rokhlin’s fast multipole method[4] leads to a further reduction to O(N) by introducing the concept of local expansion. The fast multipole BEM has been investigated by many authors including Yamada et al.[5] and Peirce et al.[6] for 2-D elastostatics, Fu et al.[7] and Popov et al.[8] for 3-D elasticity problems, and Yoshida et al.[9] for crack problems for the 3-D Laplace equation. This paper describes the application of FMM based on complex Taylor series expansions to TBIE for 2-D crack analysis. An improved preconditioner and a new initial candidate vector were also developed for the
77
WANG Pengbo (王朋波) et al：Fast Multipole BEM for Simulation of 2-D Solids ……
generalized minimum residual method solver (GMRES). Large numbers of cracks in infinite space were efficiently and accurately simulated. The largest model has 900 000 degrees of freedom (DOFs), and the total analysis time on one personal computer (PC) is only 6 h.
point x . For a smooth boundary, cαβ ( x ) = δ αβ / 2 .
1 Boundary Integral Equations for Crack Analysis
tion gives the following traction boundary integral equation.
For a 2-D elastic solid with a single ideal crack, the following integral equations can be obtained:
uα ( X ) = ∫ (Uαβ ( X , y ) tβ ( y )－Tαβ ( X , y ) uβ ( y ) )i Γo
dΓ ( y ) + ∫ − Tαβ ( X , y ) ∆uβ ( y ) dΓ ( y ), Γc
X ∈D
(1)
where Γ o is the outer boundary of the 2-D solid, Γ c− is one of the two co-planar surfaces of the crack Γ c , and U αβ ( X , y ), Tαβ ( X , y ) are the kernel functions for the 2-D elastostatics. X is the source point in the interior of the solid and y is the field point on the boundary, uβ and t β are the displacement and the traction on outer boundary Γ o . ∆uα is the crack opening displacement (COD) of the crack. The analysis assumes that the crack does not experience pressure loads, so the traction on the crack surfaces is zero. As Γ o tends to infinity, Eq. (1) becomes uα ( X ) = uα∞ ( X ) + ∫ − Tαβ ( X , y ) ∆uβ ( y) dΓ ( y), Γc
X ∈D
(2)
Sγαβ ( x , y ) is a hyper-singular kernel defined by
Sγαβ ( x , y ) = Cαβ kl
∂Tkγ ( x , y ) ∂xl
(5)
Finally, the assumption that Γ c is free from trac-
tβ ( x ) = nα ( x)∫ − Sγαβ ( x, y ) ∆uγ ( y ) dΓ ( y ) + Γc
tβ∞ ( x ) ,
(6)
x∈S
Note that the kernel Sγαβ ( x , y ) has a hypersingularity of O( r −2 ). The difficulty caused by hypersingular integrals was overcome using the concept of the finite part integral[10]. The finite part integral requires that the displacement fields be Hölder continuous at each node, so this analysis uses discontinuous elements. Because all the nodes are in the interior of the discontinuous elements, the Hölder continuity is naturally satisfied. With this technique, the hyper-singular integrals are accurately and efficiently calculated.
2 Fast Multipole BEM Formulation 2.1
Multipole expansion
In the multipole expansion, the boundary integral of the kernels is expanded into a series of products of functions. Consider the following equation,
∫Γ
∞
where uα is the so-called “no crack solution”, which is the displacement at point X when the existence of the crack is not taken into consideration. Now considering the relationship between the traction and the displacement ti ( x ) = n j ( x )C jikl uk ,l ( x ) (3)
∫
Γj
Sγαβ ( x , y )∆uγ ( y )dΓ ( y ) = j
Cαβµν
Cαβµν
∂Tµγ ( x , y ) ∂xν ( x )
∫
−
Γc
∆uγ ( y )dΓ ( y ) =
Tµγ ( x , y )∆uγ ( y )dΓ ( y ) ∂xν ( x )
(7)
and letting the inner source point X approach x , which is a point on Γ c− , give the following traction boundary integral equation:
We do not expand the relatively complicated Sγαβ ( x, y ) directly but expand the kernel Tαβ ( x, y ) .
cαβ ( x )tα ( x ) = nα ( x ) ∫ − Sγαβ ( x, y ) ∆uγ ( y ) dΓ ( y ) +
complex Taylor series around a selected point y0 as
Γc
tβ∞ ( x ) ,
x∈S
(4)
where t β∞ is the traction associated with uβ∞ and
cαβ ( x ) is related to the geometry at the boundary
The kernel function Tαβ ( x, y ) can be expanded into ∞
r r Tαβ ( x, y) = ∑ Re ⎡⎣ f ( ) ( x − y0 , k ) g ( ) ( y − y0 , k ) ⎤⎦ + k =0
∞
∑ Re ⎡⎣ f ( ) ( x − y , k ) g ( ) ( y − y , k )⎤⎦ k =0
i
i
0
0
(8)
78
Tsinghua Science and Technology, February 2005, 10(1): 76–81
where x should be far enough from y to satisfy
where l ( ) ( x, k ) = Re( x ) x k ,
y − y0 ≤ x − y0 / 2 .
r
From Eq. (8), the boundary integral of kernel Tαβ ( x, y ) can be given as
∫Γ
Tαβ ( x , y )∆u β dΓ ( y ) =
∑ Re ⎡⎣ f ( ) ( x − y , k )C ( r
fr )
0
k =0
∞
fi
0
k =0
(9)
0
f
( li )
−k
( x , k ) = Re( x ) x , −k
f ( x , k ) = Im( x ) x ,
( y0 , k ) and C
( fi )
(10)
k = 1, 2,...
( y0 , k ) are multipole moments
centered at y0 . In the fast multipole BEM for 2-D elastostatics, Yamada et al. used 6 groups of multipole moments[3]. But this analysis uses only 2 groups of multipole moments ( C ( fr ) , C ( fi ) ). The numerical examples show that the reduction of the multipole moments obviously improves the computational efficiency. 2.2
from the multipole moments centered at y0 as p
∑C ( x0 − y0 )
−l
k −1 l + k −1
−k
( fi)
C ( y0, k ) +
k =1
(r)
(i)
C
li
D ( x0 , l ) = ( y0 − x0 )
where
( fr )
lr D( ) ( x0 , k ) and D( ) (x0,k) are local moments centered
Eq. (12). The local moments centered at x0 can be obtained
( y0 , k ) ⎤⎦ +
∑ Re ⎡⎣ f ( ) ( x − y , k )C ( ) ( y , k )⎤⎦ i
(13)
i
at x0 . Only two groups of local moments are used in
j
∞
l ( ) ( x , k ) = Im( x ) x k , k = 0,1, 2,...
i ( y0 − x0 )
− l −1
p
− y0 ) Re ( x0 − y0 ) C
∑C ( x k −1 l +k
−k
0
( fr )
( y ,k) + 0
k =1
i ( y0 − x0 )
− l −1
p
∑C ( x k −1 l +k
0
− y0 ) Im ( x0 − y0 ) C −k
( fi )
( y ,k) 0
(14)
k =1
This translation is called multipole moment to local expansion translation (M2L). 2.4
Local expansion to local expansion translation
If the center of the local expansion center in Eq. (12) is shifted to x1 from x0 , the new local moments can be obtained from the old ones as
Multipole moment to multipole moment translation
D
( li )
p
( x1 , l ) = ∑ C lk ( x1 − x0 ) k − l D ( ) ( x0 , k ) + li
k =l
If the multipole expansion shifts to y1 from y0 in Eq. (9), the new multipole moments can be obtained from the original ones as C
( fi )
k −1
( y1 , k ) = ∑ C C q k −1
( fi )
i ∑ C qk − 2 ⎡ Re( y1 − y0 )C q =0
⎣
( fr )
k −2
i ∑ C qk − 2 ⎡ Im( y1 − y0 )C
⎣
( y0 , k − q )( y0 − y1 ) +
q =0
( y0 , k − q − 1) ⎤⎦ ( y0 − y1 ) q +
( fi )
( y0 , k − q − 1) ⎤⎦ ( y0 − y1 ) q (11)
This translation between multipole moments, which is applied when the multipole-expansion center shifts, is called multipole moment to multipole moment translation (M2M). 2.3
Local expansion
Expansion of the left-hand side of Eq. (9) with respect to the source point x around a selected point x0 gives:
∫Γ Tαβ ( x, y)∆uβ dΓ ( y) = j
∞
∑ Re ⎡⎣l ( ) ( x − x , k ) D( ) ( x , k ) +l ( ) ( x − x , k ) D( ) ( x , k )⎤⎦ (12) k =0
r
lr
0
i
0
li
0
i ∑ C ( x1 − x0 ) k − l −1 ⎡ Re( x1 − x0 ) D l +1 k
k = l +1
0
⎣
p
( lr )
i ∑ C lk+1 ( x1 − x0 ) k − l −1 ⎡ Im( x1 − x0 ) D k = l +1
q
q =0
k −2
p
⎣
( x0 , k ) ⎤⎦ +
( li )
( x0 , k ) ⎤
⎦
(15)
This translation between local moments, which is applied when the local-expansion center shifts, is called local expansion to local expansion translation (L2L).
3 Numerical Implementation The first step of the fast multipole algorithm is the construction of a quad-tree structure. The second step is the iteration process, which includes computation of the multipole moments (Upward) and the local moments (Downward). The Upward and Downward computations are executed iteratively until the defined accuracy is achieved[4,9]. This analysis used the GMRES as the iterative solver. An improved preconditioner and a new initial candidate vector were used to further improve the computational efficiency. A block diagonal matrix was used as the preconditioner for GMRES. For each crack, the preconditioner
WANG Pengbo (王朋波) et al：Fast Multipole BEM for Simulation of 2-D Solids ……
79
omitted the effect of other cracks. Considering only a single crack results in an equation system in matrix form as follows: M∆u = Τ ∞
(16)
where M is the coefficient matrix, ∆u is the COD solution when omitting the effects of other cracks and ∞
is called the “single crack solution,” and T is the traction vector associated with the “no crack solution” ∞
uα (see Eqs. (2) and (4)).
For each crack, ∆u and the inverse of M can be calculated. M −1 and ∆u of all the cracks compose the preconditioner and the initial vector, respectively. Cracks with identical shapes and identical mesh distributions have the same M −1 in local coordinates. Thus the computational cost and memory requirements for the preconditioner and the initial vector are so small that they can be essentially ignored. The numerical examples show that the preconditioner and the initial vector can reduce the CPU time and memory cost by 50%-70%.
Fig. 1 Comparison of the CODs from fast multipole BEM and elastic theory
4.2
Example 2 (example with four cracks)
Consider an infinite 2-D elastic solid with four line cracks as shown in Fig. 2. The crack length 2a = 16 mm , the material properties are G = 100 MPa and v = 0.3, and the stresses at infinity are
σ 11∞ = σ 22∞ = 10 MPa .
4 Numerical Examples 4.1
Example 1 (single crack)
The computational model contains a single line crack in an infinite 2-D elastic solid. For the crack length 2a = 16 mm, the material properties are G = 100 MPa and v = 0.3, and the stresses at infinity are
σ 11∞ = σ 22∞ = 20 MPa . From the elastic theory, for a crack along the x-axis, if the origin is set at the center of the crack, the COD can be analytically expressed as ∆u 2 ( x ) =
2(1 − v)a ∞ x2 σ 22 1 − 2 G a
Fig. 2
(17)
In the numerical analysis, the crack was divided into 16 quadratic elements. The order of finite series was taken as p=25. Figure 1 compares the variation of the COD on the crack boundary predicted by the present fast multipole BEM with the elastic theory. The numerical results agree well with the elastic analytical solution. The maximum error of the fast multipole BEM results is less than 1.0%.
Four cracks in an infinite 2-D elastic solid
The example was analyzed using the present fast multipole BEM and the commercial FEM software MSC/Marc. In the fast multipole BEM analysis, each crack was divided into 32 quadratic elements and the order of the finite series was taken as p=20. The FEM analysis using MSC/Marc considered a finite square plate with four cracks in the center where the plate width was 20 times the crack length. In total 13 000 six-node triangular elements were used. Comparison of the COD from the two numerical analyses in Fig. 3 shows that the fast multipole BEM results agree well with the FEM results.
80
Tsinghua Science and Technology, February 2005, 10(1): 76–81
The cracks were distributed in a square region with 900 000 DOFs. The order of the finite series was 20. The total CPU time on a desktop PC was 6.15 h. Figure 6 displays the COD results for part of the computational region.
Fig. 3 Comparison of the CODs from the FM-BEM and MSC/Marc analysis
4.3
Example 3 (comparison of the conventional BEM with the fast multipole BEM)
Various numbers of regularly distributed line cracks in an infinite 2-D elastic solid subjected to the uniform ∞ tension σ 11∞ = σ 22 were analyzed using the conven-
Fig. 5
3000 irregularly distributed cracks
tional BEM and the fast multipole BEM. DOFs were varied from 1000 to 20 000. Figure 4 compares the solution speeds of the conventional BEM (ordinary Gauss elimination) and the fast multipole BEM for various problem sizes. In the fast multipole BEM analysis, the order of the finite series was p=20. The results show that fast multipole BEM is much more efficient than the conventional BEM for large-scale problems.
Fig. 6
COD for part of the region
5 Conclusions
Fig. 4 Comparison of CPU time between the conventional BEM and the fast multipole BEM
4.4
Example 4 cracks)
(3000 irregularly distributed
This example simulated 3000 irregularly distributed cracks in an infinite 2-D elastic solid as shown in Fig.5.
The fast multipole BEM was applied to 2-D crack analysis to significantly reduce both the computational complexity and the memory requirement. Numerical examples show that the technique is very efficient for large-scale crack problems. A new form of complex Taylor series expansion was used in the multipole expansions which not only reduce the programming complexity but also improve computational efficiency. The improved preconditioner and the new initial vector for the GMRES solver lead
WANG Pengbo (王朋波) et al：Fast Multipole BEM for Simulation of 2-D Solids ……
to further reduction of the CPU time and memory requirements. This analysis considered only line cracks in an infinite domain. However, the fast multipole BEM scheme can be easily extended to finite domains and cracks of various shapes. Further investigation will simulate interface cracks in composite materials and the fracture process with the fast multipole BEM. References [1] Cruse T A. BIE fracture mechanics analysis 25 years of developments. Comput. Mech., 1996, 18(1): 1-11. [2] Portela A, Aliabadi M H, Rooke D P. Dual boundary element incremental analysis of crack propagation. Computers and Structures, 1993, 46(2): 237-247. [3] Barnes J, Hut P. A hierarchical O(NlogN) force calculation algorithm. Nature, 1986, 324: 446-449. [4] Greengard L, Rokhlin V. A fast algorithm for particle simulation. J. Comput. Phys., 1997, 135: 280-292. [5] Yamada Y, Hayami K. A multipole boundary element method for two-dimensional elastostatics. Technical
81
Report, Department of Mathematical Engineering and Information Physics, The University of Tokyo, 1995, METR 95-07: 1-20. [6] Peirce A P, Napier J A L. A spectral multipole method for efficient solutions of large scale boundary element models in elastostatics. Int. J. Numer. Methods Engng., 1995, 38: 4009-4034. [7] Fu Y H, Klimkowski K J, Rodin G J, et al. A fast solution method for three-dimensional many-particle problems of linear elasticity. Int. J. Numer. Methods Engng., 1998, 42: 1215-1229. [8] Popov V, Power H. An O(N) Taylor series multipole boundary element method for three-dimensional elasticity problems. Engineering Analysis with Boundary Elements, 2001, 25: 7-18. [9] Yoshida K, Nishimura N, Kobayashi S. Application of new fast multipole boundary integral equation method to crack problems in 3D. Engineering Analysis with Boundary Elements, 2001, 25: 239-247. [10] Guiggiani M. A general algorithm for the numerical solution of hypersingular boundary integral equations. J. Appl. Mech., 1992, 59: 604-614.
Fast Multipole BEM for Simulation of 2-D Solids Containing Large Numbers of Cracks WANG Pengbo (王朋波), YAO Zhenhan (姚振汉) **, WANG Haitao (王海涛) Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China Abstract:
The fast multipole method was used to solve the traction boundary integral equation for 2-D
crack analysis. The use of both multipole and local expansions reduces both the computational complexity and the memory requirement to O(N). The multipole expansion uses a complex Taylor series expansion to reduce the number of multipole moments. The generalized minimum residual method solver (GMRES) was selected as the iterative solver. An improved preconditioner for GMRES was developed which uses less CPU time and less memory. A new initial candidate vector for the iterative solver was developed to further improve the efficiency. The numerical examples apply the method to the analysis of cracks in infinite 2-D space with the largest model having 900 000 degrees of freedom. Key words: fast multipole method; boundary element method (BEM); cracks
Introduction The boundary element method (BEM) is recognized as a powerful method for the study of crack problems because of its semi-analytical nature and boundary-only discretization[1]. Compared with other numerical methods, the reduction in dimensionality dramatically reduces the initial data preparation and remeshing work during crack growth. Additionally, since there are no internal approximations, the singular stress field around the crack tip can be analyzed more accurately and more effectively. However, intrinsic basis for the displacement boundary integral equation (DBIE) causes the equation degeneration when the two surfaces of one crack are considered co-planar. The degeneration caused by the co-planar surfaces can be overcome by using the traction boundary integral equation (TBIE) on one of the two surfaces[2]. The TBIE can be used to develop a single-region formulation for general crack analysis. Received: 2004-08-16; revised: 2004-08-30
﹡﹡ To whom correspondence should be addressed. E-mail: [email protected]; Tel: 86-10-62772913
Conventional BEM is not efficient for large-scale problems because of the dense and asymmetric matrix. The memory requirements and the CPU time have been reduced through the application of the fast multipole method (FMM) to the BEM. Barnes and Hut[3] first used a tree data structure and the concept of multipole expansion to calculate the matrix-vector product without explicitly forming the matrix. The computational complexity and memory requirements for the matrix-vector multiplication are reduced from O(N2) to O(NlogN). Greengard and Rokhlin’s fast multipole method[4] leads to a further reduction to O(N) by introducing the concept of local expansion. The fast multipole BEM has been investigated by many authors including Yamada et al.[5] and Peirce et al.[6] for 2-D elastostatics, Fu et al.[7] and Popov et al.[8] for 3-D elasticity problems, and Yoshida et al.[9] for crack problems for the 3-D Laplace equation. This paper describes the application of FMM based on complex Taylor series expansions to TBIE for 2-D crack analysis. An improved preconditioner and a new initial candidate vector were also developed for the
77
WANG Pengbo (王朋波) et al：Fast Multipole BEM for Simulation of 2-D Solids ……
generalized minimum residual method solver (GMRES). Large numbers of cracks in infinite space were efficiently and accurately simulated. The largest model has 900 000 degrees of freedom (DOFs), and the total analysis time on one personal computer (PC) is only 6 h.
point x . For a smooth boundary, cαβ ( x ) = δ αβ / 2 .
1 Boundary Integral Equations for Crack Analysis
tion gives the following traction boundary integral equation.
For a 2-D elastic solid with a single ideal crack, the following integral equations can be obtained:
uα ( X ) = ∫ (Uαβ ( X , y ) tβ ( y )－Tαβ ( X , y ) uβ ( y ) )i Γo
dΓ ( y ) + ∫ − Tαβ ( X , y ) ∆uβ ( y ) dΓ ( y ), Γc
X ∈D
(1)
where Γ o is the outer boundary of the 2-D solid, Γ c− is one of the two co-planar surfaces of the crack Γ c , and U αβ ( X , y ), Tαβ ( X , y ) are the kernel functions for the 2-D elastostatics. X is the source point in the interior of the solid and y is the field point on the boundary, uβ and t β are the displacement and the traction on outer boundary Γ o . ∆uα is the crack opening displacement (COD) of the crack. The analysis assumes that the crack does not experience pressure loads, so the traction on the crack surfaces is zero. As Γ o tends to infinity, Eq. (1) becomes uα ( X ) = uα∞ ( X ) + ∫ − Tαβ ( X , y ) ∆uβ ( y) dΓ ( y), Γc
X ∈D
(2)
Sγαβ ( x , y ) is a hyper-singular kernel defined by
Sγαβ ( x , y ) = Cαβ kl
∂Tkγ ( x , y ) ∂xl
(5)
Finally, the assumption that Γ c is free from trac-
tβ ( x ) = nα ( x)∫ − Sγαβ ( x, y ) ∆uγ ( y ) dΓ ( y ) + Γc
tβ∞ ( x ) ,
(6)
x∈S
Note that the kernel Sγαβ ( x , y ) has a hypersingularity of O( r −2 ). The difficulty caused by hypersingular integrals was overcome using the concept of the finite part integral[10]. The finite part integral requires that the displacement fields be Hölder continuous at each node, so this analysis uses discontinuous elements. Because all the nodes are in the interior of the discontinuous elements, the Hölder continuity is naturally satisfied. With this technique, the hyper-singular integrals are accurately and efficiently calculated.
2 Fast Multipole BEM Formulation 2.1
Multipole expansion
In the multipole expansion, the boundary integral of the kernels is expanded into a series of products of functions. Consider the following equation,
∫Γ
∞
where uα is the so-called “no crack solution”, which is the displacement at point X when the existence of the crack is not taken into consideration. Now considering the relationship between the traction and the displacement ti ( x ) = n j ( x )C jikl uk ,l ( x ) (3)
∫
Γj
Sγαβ ( x , y )∆uγ ( y )dΓ ( y ) = j
Cαβµν
Cαβµν
∂Tµγ ( x , y ) ∂xν ( x )
∫
−
Γc
∆uγ ( y )dΓ ( y ) =
Tµγ ( x , y )∆uγ ( y )dΓ ( y ) ∂xν ( x )
(7)
and letting the inner source point X approach x , which is a point on Γ c− , give the following traction boundary integral equation:
We do not expand the relatively complicated Sγαβ ( x, y ) directly but expand the kernel Tαβ ( x, y ) .
cαβ ( x )tα ( x ) = nα ( x ) ∫ − Sγαβ ( x, y ) ∆uγ ( y ) dΓ ( y ) +
complex Taylor series around a selected point y0 as
Γc
tβ∞ ( x ) ,
x∈S
(4)
where t β∞ is the traction associated with uβ∞ and
cαβ ( x ) is related to the geometry at the boundary
The kernel function Tαβ ( x, y ) can be expanded into ∞
r r Tαβ ( x, y) = ∑ Re ⎡⎣ f ( ) ( x − y0 , k ) g ( ) ( y − y0 , k ) ⎤⎦ + k =0
∞
∑ Re ⎡⎣ f ( ) ( x − y , k ) g ( ) ( y − y , k )⎤⎦ k =0
i
i
0
0
(8)
78
Tsinghua Science and Technology, February 2005, 10(1): 76–81
where x should be far enough from y to satisfy
where l ( ) ( x, k ) = Re( x ) x k ,
y − y0 ≤ x − y0 / 2 .
r
From Eq. (8), the boundary integral of kernel Tαβ ( x, y ) can be given as
∫Γ
Tαβ ( x , y )∆u β dΓ ( y ) =
∑ Re ⎡⎣ f ( ) ( x − y , k )C ( r
fr )
0
k =0
∞
fi
0
k =0
(9)
0
f
( li )
−k
( x , k ) = Re( x ) x , −k
f ( x , k ) = Im( x ) x ,
( y0 , k ) and C
( fi )
(10)
k = 1, 2,...
( y0 , k ) are multipole moments
centered at y0 . In the fast multipole BEM for 2-D elastostatics, Yamada et al. used 6 groups of multipole moments[3]. But this analysis uses only 2 groups of multipole moments ( C ( fr ) , C ( fi ) ). The numerical examples show that the reduction of the multipole moments obviously improves the computational efficiency. 2.2
from the multipole moments centered at y0 as p
∑C ( x0 − y0 )
−l
k −1 l + k −1
−k
( fi)
C ( y0, k ) +
k =1
(r)
(i)
C
li
D ( x0 , l ) = ( y0 − x0 )
where
( fr )
lr D( ) ( x0 , k ) and D( ) (x0,k) are local moments centered
Eq. (12). The local moments centered at x0 can be obtained
( y0 , k ) ⎤⎦ +
∑ Re ⎡⎣ f ( ) ( x − y , k )C ( ) ( y , k )⎤⎦ i
(13)
i
at x0 . Only two groups of local moments are used in
j
∞
l ( ) ( x , k ) = Im( x ) x k , k = 0,1, 2,...
i ( y0 − x0 )
− l −1
p
− y0 ) Re ( x0 − y0 ) C
∑C ( x k −1 l +k
−k
0
( fr )
( y ,k) + 0
k =1
i ( y0 − x0 )
− l −1
p
∑C ( x k −1 l +k
0
− y0 ) Im ( x0 − y0 ) C −k
( fi )
( y ,k) 0
(14)
k =1
This translation is called multipole moment to local expansion translation (M2L). 2.4
Local expansion to local expansion translation
If the center of the local expansion center in Eq. (12) is shifted to x1 from x0 , the new local moments can be obtained from the old ones as
Multipole moment to multipole moment translation
D
( li )
p
( x1 , l ) = ∑ C lk ( x1 − x0 ) k − l D ( ) ( x0 , k ) + li
k =l
If the multipole expansion shifts to y1 from y0 in Eq. (9), the new multipole moments can be obtained from the original ones as C
( fi )
k −1
( y1 , k ) = ∑ C C q k −1
( fi )
i ∑ C qk − 2 ⎡ Re( y1 − y0 )C q =0
⎣
( fr )
k −2
i ∑ C qk − 2 ⎡ Im( y1 − y0 )C
⎣
( y0 , k − q )( y0 − y1 ) +
q =0
( y0 , k − q − 1) ⎤⎦ ( y0 − y1 ) q +
( fi )
( y0 , k − q − 1) ⎤⎦ ( y0 − y1 ) q (11)
This translation between multipole moments, which is applied when the multipole-expansion center shifts, is called multipole moment to multipole moment translation (M2M). 2.3
Local expansion
Expansion of the left-hand side of Eq. (9) with respect to the source point x around a selected point x0 gives:
∫Γ Tαβ ( x, y)∆uβ dΓ ( y) = j
∞
∑ Re ⎡⎣l ( ) ( x − x , k ) D( ) ( x , k ) +l ( ) ( x − x , k ) D( ) ( x , k )⎤⎦ (12) k =0
r
lr
0
i
0
li
0
i ∑ C ( x1 − x0 ) k − l −1 ⎡ Re( x1 − x0 ) D l +1 k
k = l +1
0
⎣
p
( lr )
i ∑ C lk+1 ( x1 − x0 ) k − l −1 ⎡ Im( x1 − x0 ) D k = l +1
q
q =0
k −2
p
⎣
( x0 , k ) ⎤⎦ +
( li )
( x0 , k ) ⎤
⎦
(15)
This translation between local moments, which is applied when the local-expansion center shifts, is called local expansion to local expansion translation (L2L).
3 Numerical Implementation The first step of the fast multipole algorithm is the construction of a quad-tree structure. The second step is the iteration process, which includes computation of the multipole moments (Upward) and the local moments (Downward). The Upward and Downward computations are executed iteratively until the defined accuracy is achieved[4,9]. This analysis used the GMRES as the iterative solver. An improved preconditioner and a new initial candidate vector were used to further improve the computational efficiency. A block diagonal matrix was used as the preconditioner for GMRES. For each crack, the preconditioner
WANG Pengbo (王朋波) et al：Fast Multipole BEM for Simulation of 2-D Solids ……
79
omitted the effect of other cracks. Considering only a single crack results in an equation system in matrix form as follows: M∆u = Τ ∞
(16)
where M is the coefficient matrix, ∆u is the COD solution when omitting the effects of other cracks and ∞
is called the “single crack solution,” and T is the traction vector associated with the “no crack solution” ∞
uα (see Eqs. (2) and (4)).
For each crack, ∆u and the inverse of M can be calculated. M −1 and ∆u of all the cracks compose the preconditioner and the initial vector, respectively. Cracks with identical shapes and identical mesh distributions have the same M −1 in local coordinates. Thus the computational cost and memory requirements for the preconditioner and the initial vector are so small that they can be essentially ignored. The numerical examples show that the preconditioner and the initial vector can reduce the CPU time and memory cost by 50%-70%.
Fig. 1 Comparison of the CODs from fast multipole BEM and elastic theory
4.2
Example 2 (example with four cracks)
Consider an infinite 2-D elastic solid with four line cracks as shown in Fig. 2. The crack length 2a = 16 mm , the material properties are G = 100 MPa and v = 0.3, and the stresses at infinity are
σ 11∞ = σ 22∞ = 10 MPa .
4 Numerical Examples 4.1
Example 1 (single crack)
The computational model contains a single line crack in an infinite 2-D elastic solid. For the crack length 2a = 16 mm, the material properties are G = 100 MPa and v = 0.3, and the stresses at infinity are
σ 11∞ = σ 22∞ = 20 MPa . From the elastic theory, for a crack along the x-axis, if the origin is set at the center of the crack, the COD can be analytically expressed as ∆u 2 ( x ) =
2(1 − v)a ∞ x2 σ 22 1 − 2 G a
Fig. 2
(17)
In the numerical analysis, the crack was divided into 16 quadratic elements. The order of finite series was taken as p=25. Figure 1 compares the variation of the COD on the crack boundary predicted by the present fast multipole BEM with the elastic theory. The numerical results agree well with the elastic analytical solution. The maximum error of the fast multipole BEM results is less than 1.0%.
Four cracks in an infinite 2-D elastic solid
The example was analyzed using the present fast multipole BEM and the commercial FEM software MSC/Marc. In the fast multipole BEM analysis, each crack was divided into 32 quadratic elements and the order of the finite series was taken as p=20. The FEM analysis using MSC/Marc considered a finite square plate with four cracks in the center where the plate width was 20 times the crack length. In total 13 000 six-node triangular elements were used. Comparison of the COD from the two numerical analyses in Fig. 3 shows that the fast multipole BEM results agree well with the FEM results.
80
Tsinghua Science and Technology, February 2005, 10(1): 76–81
The cracks were distributed in a square region with 900 000 DOFs. The order of the finite series was 20. The total CPU time on a desktop PC was 6.15 h. Figure 6 displays the COD results for part of the computational region.
Fig. 3 Comparison of the CODs from the FM-BEM and MSC/Marc analysis
4.3
Example 3 (comparison of the conventional BEM with the fast multipole BEM)
Various numbers of regularly distributed line cracks in an infinite 2-D elastic solid subjected to the uniform ∞ tension σ 11∞ = σ 22 were analyzed using the conven-
Fig. 5
3000 irregularly distributed cracks
tional BEM and the fast multipole BEM. DOFs were varied from 1000 to 20 000. Figure 4 compares the solution speeds of the conventional BEM (ordinary Gauss elimination) and the fast multipole BEM for various problem sizes. In the fast multipole BEM analysis, the order of the finite series was p=20. The results show that fast multipole BEM is much more efficient than the conventional BEM for large-scale problems.
Fig. 6
COD for part of the region
5 Conclusions
Fig. 4 Comparison of CPU time between the conventional BEM and the fast multipole BEM
4.4
Example 4 cracks)
(3000 irregularly distributed
This example simulated 3000 irregularly distributed cracks in an infinite 2-D elastic solid as shown in Fig.5.
The fast multipole BEM was applied to 2-D crack analysis to significantly reduce both the computational complexity and the memory requirement. Numerical examples show that the technique is very efficient for large-scale crack problems. A new form of complex Taylor series expansion was used in the multipole expansions which not only reduce the programming complexity but also improve computational efficiency. The improved preconditioner and the new initial vector for the GMRES solver lead
WANG Pengbo (王朋波) et al：Fast Multipole BEM for Simulation of 2-D Solids ……
to further reduction of the CPU time and memory requirements. This analysis considered only line cracks in an infinite domain. However, the fast multipole BEM scheme can be easily extended to finite domains and cracks of various shapes. Further investigation will simulate interface cracks in composite materials and the fracture process with the fast multipole BEM. References [1] Cruse T A. BIE fracture mechanics analysis 25 years of developments. Comput. Mech., 1996, 18(1): 1-11. [2] Portela A, Aliabadi M H, Rooke D P. Dual boundary element incremental analysis of crack propagation. Computers and Structures, 1993, 46(2): 237-247. [3] Barnes J, Hut P. A hierarchical O(NlogN) force calculation algorithm. Nature, 1986, 324: 446-449. [4] Greengard L, Rokhlin V. A fast algorithm for particle simulation. J. Comput. Phys., 1997, 135: 280-292. [5] Yamada Y, Hayami K. A multipole boundary element method for two-dimensional elastostatics. Technical
81
Report, Department of Mathematical Engineering and Information Physics, The University of Tokyo, 1995, METR 95-07: 1-20. [6] Peirce A P, Napier J A L. A spectral multipole method for efficient solutions of large scale boundary element models in elastostatics. Int. J. Numer. Methods Engng., 1995, 38: 4009-4034. [7] Fu Y H, Klimkowski K J, Rodin G J, et al. A fast solution method for three-dimensional many-particle problems of linear elasticity. Int. J. Numer. Methods Engng., 1998, 42: 1215-1229. [8] Popov V, Power H. An O(N) Taylor series multipole boundary element method for three-dimensional elasticity problems. Engineering Analysis with Boundary Elements, 2001, 25: 7-18. [9] Yoshida K, Nishimura N, Kobayashi S. Application of new fast multipole boundary integral equation method to crack problems in 3D. Engineering Analysis with Boundary Elements, 2001, 25: 239-247. [10] Guiggiani M. A general algorithm for the numerical solution of hypersingular boundary integral equations. J. Appl. Mech., 1992, 59: 604-614.