An Introduction to the

Boundary Element Method (BEM)

and Its Applications in
Modeling Composite Materials
Yijun Liu
Department of Mechanical, Industrial and Nuclear Engineering
University of Cincinnati, P.O. Box 210072
Cincinnati, Ohio 45221-0072, U.S.A.
E-mail: Yijun.Liu@uc.edu

September 23, 2004

1 CAE Research Lab
 Outline

• An introduction to the Boundary Element Method (BEM)

• Applications of the BEM in solving engineering problems

• BEM in large-scale modeling of fiber-reinforced composites

• Discussions

• References

• Acknowledgements

• Further Information

2 CAE Research Lab

 An Introduction to the BEM
- Two Different Approaches in Computational Mechanics

Engineering Problems

Mathematical Models

Differential Equation (ODE/PDE) (Boundary) Integral Equation (BIE)

Formulations Formulations

Analytical Solutions Numerical Solutions Analytical Solutions Numerical Solutions

FDM FEM EFM Others BEM Others

3 CAE Research Lab

 A Brief History of the BEM

Jaswon and Symm (1963)

– 2D Potential Problems

Integral equations T. A. Cruse and F. J. Rizzo (1968)

(Fredholm, 1903) – 2D elastodynamics
BEM emerged in 1980’s …

Modern numerical P. K. Banerjee (1975)

solutions of BIEs – coined the name “boundary element method”
(in early 1960’s) (this has been disputed by others)
F. J. Rizzo (1964, paper 1967)
– 2D Elasticity Problems

4 CAE Research Lab

Advantages of the BEM and the Mysteries
Advantages:
• Accuracy – due to the semi-analytical nature and use of integrals
• More efficient in modeling stage due to the reduction of dimensions
• Good for stress concentration and infinite domain problems
• Good for modeling thin shell-like structures/materials
• Neat …

Mysteries:
• BIEs are singular which are difficult to deal with (wrong!)
• BEM is slow and thus inefficient (not necessarily!)
• FEM can solve everything. Who needs BEM? (not exactly true!)

5 CAE Research Lab

 Formulation: The Potential Problem
• Governing Equation

∇ 2 u( P ) = 0, ∀P ∈V .
V S

• For 3D problems, the Green’s function is Po (xok )

r
1
G ( P , Po ) = , r = Po P .
4πr P(xk )
• BIE formulation n
E
 ∂u ∂G 
C ( Po )u( Po ) = ∫ G − u dS , ∀Po ∈ S .
S 
∂n ∂n 
• Discretization of the BIE using the boundary elements
∂u 
[ H ]{ u} = [ G ]  , or [ A]{ x} = {b} .
∂
 n

6 CAE Research Lab

 Singular or Non-Singular?
• Re-examine the BIE
 ∂u ∂G 
C ( Po )u( Po ) = ∫ G − u dS , ∀Po ∈ S .
S  ∂ n ∂ n
The second integral in the BIE is singular and is considered as a CPV integral

• However, the constant in the free term is also a CPV integral

∂G ( P , Po )
C ( Po ) = − ∫ dS ( P ).
S
∂ n
• Re-write the BIE to obtain the weakly-singular form of the BIE
∂G ∂u
∫S ∂n [ u ( P ) − u ( Po )]dS = ∫S ∂ndS , ∀Po ∈ S .
G

O (1/r2) O (r)
• Non-singular form also exists (Liu & Rudolphi, EABE, 1991 and CM, 1999)

7 CAE Research Lab

Example: Results for Heat Transfer in a Fuel Cell

Predicted Temperature Distributions Using the BEM and FEM

(a) The fuel cell model (b) BEM (max. temp. = 378.40 K) (c) FEM (max. temp. = 378.31 K)

8 CAE Research Lab

Example: Coupled Structural Acoustics Analysis
• Applications
¾ Acoustic radiation/scattering from elastic structures submerged in fluids
¾ Prediction of noises of an elastic structure in vibration
¾ Dynamics of fluid-filled elastic piping system
¾ Acoustic cavity analysis

shell structure (V)

Sb interior
exterior fluid (E) Sa
n

9 CAE Research Lab

The BIE Formulation for Structural Acoustics
• Governing Equations
¾ In elastic domain: (c12 − c22 )uk,ki ( P ) + c22ui,kk ( P ) + ω 2ui ( P ) = 0, ∀P ∈V
¾ In acoustic domain: ∇2φ ( P ) + k 2φ ( P ) = 0, ∀P ∈ E

• BIE Formulations
¾ In elastic domain: Cij (Po ) u j (Po ) = ∫ Uij ( P, Po ) t j ( P )dS ( P ) − ∫ Tij ( P, Po )u j ( P )dS ( P )
S S

¾ In acoustic domain: ∂ G( P, Po ) ∂ φ ( P) 
C( Po )φ ( Po ) = ∫  φ ( P) − G( P, Po )  dS( P) + φ I ( Po )
Sa
 ∂n ∂n 

• Interface Conditions
∂φ
¾ Velocity continuity condition across the interface: = ρ f ω 2un
∂n
¾ Stress equilibrium condition: ti = −φ ni

• Discretization of the BIE using boundary elements

10 CAE Research Lab

 Results for A Structural Acoustics Analysis
Radiated Sound Pressure from Steel Spherical Shells with Different Thickness and
Under an Internal Harmonic Pressure Load (r = 5a, M = 112)
y
r
Analytical (h/a = 0.5) BEM (h/a = 0.5)
0.25 Analytical (h/a = 0.05) BEM (h/a = 0.05)
a
Analytical (h/a = 0.01) BEM (h/a = 0.01)
θ

0.20 O
x
b

|φ/tb|0.15 z
Difficult for
FEM
0.10 y
bx
r

0.05 a
by

θ
O bz x
0.00
0 2 4 6 8 10 12

ka z

11 CAE Research Lab

Analysis of the Acoustic Fields of a Submarine

A simplified submarine model with BEM Sound pressure in the exterior domain
(Surface elements only)

12 CAE Research Lab

 BEM for Modeling Thin Layered Materials

Advantages:
• BEM is good for modeling thin shell-like materials/structures
• Much fewer elements are needed using the BEM than the FEM in the
modeling (no element connectivity and aspect-ratio restrictions).
• Accuracy.
S+
n−
S−
n+
Difficulties: Treatment of the nearly singular integrals in the BIEs.
• 3D elasticity case (Liu, IJNME, 1998)
• 2D elasticity case (Luo, Liu and Berger, CM, 1998)
• 2D piezoelectricity case (Liu and Fan, CMAME, 2002)

13 CAE Research Lab

 Analysis of Fiber-Reinforced Composites with
the Presence of the Interphases
A Unit Cell Model of Fiber-Reinforced Composites

fiber
matrix

matrix
fiber
Interphase

y
A fiber-reinforced
composite x
a unit cell
interphase

Interphases in fiber-reinforced composites are modeled using the BEM and FEM to investigate
their effects on the mechanical properties and interface failures of the material.

14 CAE Research Lab

Analysis of Fiber-Reinforced Composites with
the Presence of the Interphases (Cont.)

BEM (384 quadratic line elements) FEM (10,188 quadratic area elements)

15 CAE Research Lab

Analysis of Fiber-Reinforced Composites with
the Presence of the Interphases (Cont.)

8.0
Without interphase
7.0 With interphase (h=0.01a)
With interphase (h=0.03a)
6.0
Other BEM (Pan et al. 1998)
Self-consistent (Oshima and Watari 1992)
5.0

Ex / E(m)
4.0

3.0

2.0

1.0

0.0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Fiber volume fraction (Vf)

Stress distribution Effective Young’s modulus

16 CAE Research Lab

Analysis of Fiber-Reinforced Composites with
the Presence of the Interphases (Cont.)
A Circular-Arc Crack Between the Interphase and the Matrix
y
1.0

interphase 0.8
crack

Normalized SIF K1
h 0.6

α
R α x Material One (same as the matrix)
0.4
Material Two (E=3.0E+06 psi)
fiber
Material Three (E=5.25E+06 psi)

0.2
Material Four (E=8.5E+06 psi)
Material Five (same as the fiber)
matrix Sub-interface crack
No interphase
0.0
δ 0 15 30 45 60 75 90
Angle alpha (deg.)

Effects of the interphase materials on the stress intensity

factor (K1 / σ ave π Rα ) for the circular-arc interface crack

17 CAE Research Lab

 BEM for Thin Piezoelectric Solids
Applications of Piezoelectric Materials
• Thin piezo films and coatings as sensors/actuators in smart materials
• Micro-electro-mechanical systems (MEMS)
•…

The mechanical and electrical coupling effect in piezoelectric materials

18 CAE Research Lab

BEM for Thin Piezoelectric Solids (Cont.)
BIE for piezoelectricity (weakly-singular form):

∫ T( P, P )[u( P) − u( P )]dS ( P) = ∫ U( P, P )t( P)dS ( P)

S
0 0
S
0

∫
+ U( P, P )b( P )dV ( P ),
V
0 ∀P ∈ S ,
0

for a finite piezoelectric solid, in which (for 2D case):

 u1   t1   f1 
     
u =  u2 , t =  t 2 , b =  f 2 ,
− φ  − ω  − q 
     
U 11 U 12 Φ1  T11 T12 Ω1 
U = U 21 U 22 Φ 2 , T = T21 T22 Ω 2 .
   
U 31 U 32 Φ 3  T31 T32 Ω 3 

19 CAE Research Lab

 BEM for Thin Piezoelectric Solids (Cont.)
A PZT-5 Strip Subjected to Pressure Load P and Voltage V0
x3
p 1.2E-03
Analyti cal (2h/L = 1.0)
Piezo BEM (2h/L = 1.0)
1.0E-03 Analytical (2h/L = 0.01)
Piezo BEM (2h/L = 0.01)
Analytical (2h/L = 0.0001)

D isplacem en t w (mm )
+V0 -V0 h 8.0E-04
Piezo BEM (2h/L = 0.0001)

x1 6.0E-04

L h
4.0E-04

2.0E-04

p
0.0E+00
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

x (mm)

Displacement component w along the bottom edge

of the strip (M = 24) for different thicknesses

20 CAE Research Lab

BEM for Thin Piezoelectric Solids (Cont.)
Analysis of Piezoelectric Parallel Bimorph
(Bending deformation with applied voltage)
x3
V = V0 V=0
1.00

Layer 1 h
Layer 2 h x1 0.0 2.0 4.0 6.0 8.0 10.0

Polarization -1.00

direction
L

A parallel bimorph The deformed shape

(Note that the thickness of the layers can be made arbitrarily small without the
need to use smaller and smaller elements in the BEM)

21 CAE Research Lab

 Current Status of the BEM Research
• Fast solvers that can solve problems beyond the reach of other methods
• Large-scale analyses with DOFs above 20M
• Multi-physics and multi-scales

MEMS

Electromagnetic wave scatterings from targets

(Chew, et al., 2004)
22 CAE Research Lab
Large-Scale Modeling of Fiber-Reinforced
Composites with a
Fast Multipole Boundary Element Method

In collaboration with:
Professor Naoshi Nishimura at Kyoto University

23 CAE Research Lab

 The Approach
• A model with elastic matrix and rigid inclusions for fiber-reinforced composites
is adopted (the rigid-inclusion model)
• This model is likely to be valid for short fibers or long fibers with much higher
stiffness than that of matrix
• This approach is the first step towards more general elastic matrix/elastic fiber
models
• The fast-multipole method is used to solve the large-scale BEM equations for
this problem

24 CAE Research Lab

 Boundary Integral Equation Formulation
Representation integral:
u( x ) = ∫ [U(x, y ) t ( y ) − T( x, y )u( y )]dS ( y ) + u ∞ (x ), ∀x ∈ V (1)
S

with S = U Sα
α

For each rigid inclusion Sα :

u( y ) = d + ω × p( y ) (2)
with d and ω being the rigid-body translation and rotation, respectively

It can be shown that:

∫Sα
T( x, y )u( y )dS ( y ) = 0 (3)

for each rigid inclusion Sα

25 CAE Research Lab

 Boundary Integral Equation Formulation (Cont.)

“Simplified” BIE formulation for rigid-inclusion problems:

u( x ) = ∫ U( x, y ) t ( y )dS ( y ) + u ∞ ( x ), ∀x ∈ S (4)
S

Both u and t are unknown. Need six more equations for each inclusion

Consider the equilibrium of each inclusion (6 equations):

∫
Sα
t ( y )dS ( y ) = 0; (5)

∫Sα
p( y ) × t ( y )dS ( y ) = 0; (6)

for α = 1, 2, …, n

Eqs. (4-6) provide enough equations for solving the rigid-inclusion problem

26 CAE Research Lab

 Fast Multipole Method (FMM)

• Ranked among the top ten algorithms

of the 20th century (with FFT, QR, …)
• Developed by Rokhlin and
Greengard (mid of 1980’s)
• For 3-D elasticity: Peirce and Napier
(1995); Rodin, et al. (1997); Popov
and Power (2001), and many others
• More research (more large-scale
applications)
• Education or re-education
• A review: Nishimura, Applied
Mechanics Review, July 2002

27 (From Yoshida, 2001) CAE Research Lab

 Fast Multipole Algorithm
• The entire boundary is
divided into multi-level cells
• Each boundary element is
placed in a cell, which contains
a specified number of elements
• A tree structure of the
boundary elements is obtained
• Interactions (integrations) of
element-to-element is replaced
by those of cell-to-cell
• Expansions are employed to
accelerate the evaluations of
these interactions
(Nishimura, 2002)
28 CAE Research Lab
 Fast Multipole Expansions
Apply the following expansion:
∞ n
1
= ∑ ∑ S n ,m (Ox) Rn ,m (Oy ) (7)
r ( x, y ) n = 0 m = − n
where O represents a third point, Rn ,m and S n ,m are solid harmonic functions

Displacement kernel is written as:

∑ ∑ [F ]
∞ n
1 (8)
U ij ( x, y ) = (Ox) Rn ,m (Oy ) + Gi ,n ,m (Ox)(Oy ) j Rn ,m (Oy )
8πµ
ij ,n ,m
n =0 m = − n

which is in the form: U ~ kn(1) (Ox)kn( 2 ) (Oy )

The FMM expansion:

∑ ∑ [F ]
∞ n
1
∫So
U ij ( x, y )t j ( y )dS ( y ) =
8πµ n =0 m = − n
ij ,n ,m (Ox) M j ,n ,m (O) + Gi ,n ,m (Ox) M n ,m (O) (9)

where the four multipole moments are given by:

M j ,n ,m (O) = ∫ Rn ,m (Oy )t j ( y )dS ( y ); M n ,m (O) = ∫ (Oy ) j Rn ,m (Oy )t j ( y )dS ( y ) (10)
So So

29 CAE Research Lab

Fast Multipole Algorithm (Cont.)

x
y

(Nishimura, 2004)

30 CAE Research Lab

 A Rigid Sphere in Elastic Medium

S a 2

V
1

A sphere with tri-axial loading

A BEM mesh with 1944 constant elements

31 CAE Research Lab

 A Rigid Sphere in Elastic Medium (Cont.)
1.60

Radial stress - Analytical

Radial stress - BEM

1.40
Tangential stress - Analytical

Tangential stress - BEM

1.20
stresses

1.00

0.80

0.60
1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00
radius r/a

Radial and tangential stresses obtained by a BEM model with 120 elements
32 CAE Research Lab
A Rigid Sphere in Elastic Medium (Cont.)

Contour plot for stress on the surface of the sphere

33 CAE Research Lab
 Study of Fiber-Reinforced Composites:
The Representative Volume Element (RVE)

Matrix

Fiber (rigid inclusion)

A data-collection surface

34 CAE Research Lab

 A BEM Mesh

A BEM mesh used for the short fiber inclusion (with 456 constant elements)

35 CAE Research Lab

 Load Transfer Studies

Location of
maximum stress

A model with 216 “randomly” distributed and oriented short fibers

36 CAE Research Lab
Efficiency of the Fast Multipole BEM
100,000

10,000
CPU (sec.)

1,000

CPU time (Uniform case)

100
CPU time ("Random" case)

10
10,000 100,000 1,000,000 10,000,000
Total DOF

CPU time used for solving the BEM models for the short-fiber cases

37 CAE Research Lab

 Modeling of CNT-Based Composites (Cont.)

A small RVE containing 2,000 CNT fibers with the total DOF = 3,612,000 (CNT
length = 50 nm, volume fraction = 10.48%). A larger model with 16,000 CNT fibers
and 28.9M DOFs was solved successfully on a FUJITSU HPC2500 supercomputer
(at the Kyoto University) within 34 hours.
38 CAE Research Lab
Modeling of CNT-Based Composites (Cont.)
100.0

90.0 Odegard, et al. [23]

Aligned uniform case - small RVE

80.0
Aligned "random" case - small RVE
Effective Young's moduli E_eff (GPa)

70.0 Aligned uniform case - large RVE

60.0

50.0

40.0

30.0

20.0

10.0

0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
CNT volume fraction (%)

Computed effective moduli of CNT/polymer composites using three RVEs and

compared with NASA’s multi-scale results
39 CAE Research Lab
 Discussions

• BEM is a very efficient numerical tool for many problems

in engineering
• Computational mechanics can play a significant role in
the development of composite materials
• Multi-scale, multiphysics and large-scale approaches are
urgently needed for the development of new materials
• There are plenty of opportunities for the computational
mechanics (FEM/BEM/BNM/Meshfree methods) in
material modeling, bio-engineering and many other fields

40 CAE Research Lab

 A Bigger Picture of
Computational Solid Mechanics
FEM: Large-scale structural, nonlinear, BEM: Large-scale continuum, linear,
and transient problems and steady state (wave) problems

Meshfree: Large deformation, fracture

and moving boundary problems “If the only tool
you have is a
hammer, then
every problem
you can solve
looks like a nail!”

41 CAE Research Lab

 Future of Computational Mechanics
Large scale, multiscale, instant and visual!

An Example:
Virtual Reality (VR)
with large scale MD
simulations of a
fractured ceramic
nanocomposite
(Spheres with different
colors represent atoms
of different materials
in the nanocomposite)

(A. Nakano, et al., 2001)

42 CAE Research Lab
 References
(on BIE/BEM)

• F. J. Rizzo, "An integral equation approach to boundary value problems of classical elastostatics,"
Quart. Appl. Math., 25, 83-95 (1967).
• T. A. Cruse, "Numerical solutions in three dimensional elastostatics," Int. J. of Solids & Structures,
5, 1259-1274 (1969).
• P. K. Banerjee and R. Butterfield, Boundary Element Methods in Engineering Science (McGraw
Hill, London, 1981).
• S. Mukherjee, Boundary Element Methods in Creep and Fracture (Applied Science Publishers, New
York, 1982).
• T. A. Cruse, Boundary Element Analysis in Computational Fracture Mechanics (Kluwer Academic
Publishers, Dordrecht, The Netherlands, 1988).
• C. A. Brebbia and J. Dominguez, Boundary Elements - An Introductory Course (McGraw-Hill, New
York, 1989).
• G. S. Gipson, Boundary Element Fundamentals - Basic Concepts and Recent Developments in the
Poison Equation (Computational Mechanics Publications, Southampton, 1987).
• J. H. Kane, Boundary Element Analysis in Engineering Continuum Mechanics (Prentice Hall,
Englewood Cliffs, NJ, 1994).

43 CAE Research Lab

 References
(on composites/nanocomposites and their modeling)
• E. T. Thostenson, Z. F. Ren, and T.-W. Chou, "Advances in the science and technology of carbon
nanotubes and their composites: a review," Composites Science and Technology, 61, 1899-1912 (2001).
• G. M. Odegard, T. S. Gates, K. E. Wise, et al., "Constitutive modeling of nanotube-reinforced polymer
composites," Composites Science & Technology, 63, 1671-1687 (2003).
• M. Griebel and J. Hamaekers, "Molecular dynamics simulations of the elastic moduli of polymer-carbon
nanotube composites," Computer Methods in Appl. Mech. Engrg., 194, 1773-1788 (2004).
• F. T. Fisher, R. D. Bradshaw, and L. C. Brinson, "Fiber waviness in nanotube-reinforced polymer
composites--I: Modulus predictions using effective nanotube properties," Composites Science and
Technology, 63, 1689-1703 (2003).
• Y. J. Liu and X. L. Chen, "Continuum models of carbon nanotube-based composites using the boundary
element method," Electronic Journal of Boundary Elements (http://tabula.rutgers.edu/EJBE/), 1, 316-335
(2003).
• Y. J. Liu and X. L. Chen, "Evaluations of the effective materials properties of carbon nanotube-based
composites using a nanoscale representative volume element," Mechanics of Materials, 35, 69-81 (2003).
• X. L. Chen and Y. J. Liu, "Square representative volume elements for evaluating the effective material
properties of carbon nanotube-based composites," Computational Materials Science, 29, 1-11 (2004).
• N. Nishimura and Y. J. Liu, "Thermal analysis of carbon-nanotube composites using a rigid-line inclusion
model by the boundary integral equation method," Computational Mechanics, in press (2004).
• Y. J. Liu, N. Nishimura, Y. Otani, T. Takahashi, X. L. Chen and H. Munakata, "A fast boundary element
method for the analysis of fiber-reinforced composites based on a rigid-inclusion model," ASME Journal
of Applied Mechanics, in press (2004).
44 CAE Research Lab
 Acknowledgments
• U.S. National Science Foundation (NSF)
• Japan Society for the Promotion of Science (JSPS)
• Academic Center for Computing and Media Studies, Kyoto
University, Japan
• Chunhui Fellowship and Tsinghua University, China
• Prof. N. Nishimura at Kyoto University
• Graduate Students at the University of Cincinnati, Kyoto University
and Tsinghua University

45 CAE Research Lab

 Further Information

Website:
http://urbana.mie.uc.edu/yliu

Contact:
Dr. Yijun Liu
CAE Research Laboratory
Department of Mechanical, Industrial and Nuclear Engineering
P.O. Box 210072
University of Cincinnati
Cincinnati, OH 45221-0072
Tel.: (513) 556-4607
E-Mail: Yijun.Liu@uc.edu

46 CAE Research Lab