Journal of Natural Disaster Science, Volume 35, Number 2, 2014, pp43-53

SEISMIC BEHAVIOR ANALYSIS OF MASONRY STRUCTURES USINGTHE DISTINCT ELEMENT METHOD CONSIDERING ELEMENT DEFORMABILITY
Aiko FURUKAWA*, Risa HORIKAWA**, Junji KIYONO***, and Kenzo TOKI*****Associate Professor, Graduate School of Global Environmental Studies, Kyoto University**Graduate Student, Dept. of Urban Management, Kyoto University***Professor, Graduate School of Global Environmental Studies, Kyoto University,****Professor, Research Center for Disaster Mitigation of Urban Cultural Heritage, Ritsumeikan University(Received March 5, 2014 Accepted August 3, 2014)
ABSTRACTIn this study, we propose a new distinct element method (DEM) that considers element
deformability for failure analysis of masonry structures. Many people in developing
countries live in masonry structures. In earlier DEM schemes, a structure is modeled
as an assembly of rigid elements, but element deformability cannot be considered. The
deformation of a structure can be modeled by overlapping between elements, but Poisson’s
effect cannot be modeled. However, bricks used in developing countries can readily be
deformed due to their low stiffness. Therefore, it is preferable to also consider element
deformability in the DEM. In the new DEM, each element is divided into two parts: an
inner part that considers deformation of the element itself, and an outer part that deals
with contact between elements. Deformation of a structure can be modeled by overlapping
between elements and deformation of the elements themselves. The validity of the method
is confirmed through a comparison of the elastic deformation with a finite element model.
It was found that the original DEM and the proposed method show different failure patternsof seismic behavior due to Poisson’s effect.
Keywords: DEM, Element deformability, Masonry structure, Poisson’s effect.

1. Introduction
The city of Bam in southeastern Iran was hit by astrong earthquake on December 26, 2003. About
43,200 people lost their lives [1], many due to the collapse of masonry buildings. Masonry structures are
low-cost and construction practices are simple, somany people in developing countries live in suchbuildings. However, masonry structures are vulnerable to earthquakes. Therefore, the seismic collapse
mechanisms of such structures need to be understoodin order to develop effective reinforcement measures.
To this end, the development of a numerical methodthat can accurately simulate the seismic behavior ofmasonry structures is needed.
There are two types of numerical simulation
methods that can handle seismic behavior. One isbased on continuum modeling, and the other, on discontinuum modeling.
Among numerical simulation methods, the finite
element method (FEM) [2] is the most common technique for the analysis of a continuum. However, theFEM is based on the mechanics of the continuum anduses a continuous shape function, so it has difficulty
in solving failure and collapse phenomena.
43

A. FURUKAWA, R. HORIKAWA, J. KIYONO, K. TOKI

A new FEM, referred to as FEM-β, that can be
regarded as a FEM with a special displacement field
was proposed by Oguni et al. in 2004 [3]. In FEM-β,
the continuous shape function is replaced by a discontinuous shape function, which enables it to solve failure phenomena. However, it does not model the re-establishment of contact between elements, so it cannot
simulate collapse behavior.
A method based on discontinuum modeling ismore suitable for simulating failure and collapse behaviors. The distinct element method (DEM) [4] is anestablished numerical simulation method for a discontinuum. In the DEM, a structure is modeled as an assembly of rigid elements, and a spring and dashpot
are generated between two elements in contact. In theoriginal DEM, there are no forces resisting tension
between elements; therefore, it was impossible tosimulate the behavior of a continuum. To give continuity to this discrete model, a spring that resists tension was added and this improved method is called
the modified DEM (MDEM) [5] or the extended
DEM (EDEM) [6]. In the MDEM and EDEM, thestructure is modeled as an assembly of circular elements.
However, the elements are modeled as a rigid
body, so deformation of the elements is not considered. Deformability of the structure can be modeledby overlapping between elements. Failure is modeledby breaking the spring. A method of determining thetheoretical spring constant has not been established,so the value is generally determined experimentallyor simply based on one-dimensional (1D) wave propagation equations. Therefore, the reliability of the results is not high.
A refined DEM [7] was then proposed by Furukawa et al. in 2011. In this scheme, a restoring spring
resisting tension is included, the structure is modeled
as an assembly of hexahedral elements, and the surface of each element is divided into many segments towhich springs and dashpots are attached. This segmentation allows theoretical spring constants to be
determined using three-dimensional (3D) stress-strain
relationships. However, the elements are still modeled
as a rigid body in the same way as in the original
DEM, and the deformability of elements and Poisson’s effect cannot be expressed. Because the bricks
used in developing countries that we target in thisstudy readily deform due to their low stiffness, it ispreferable to consider the deformability of the elements themselves in the DEM. Therefore, we havedeveloped a new DEM that considers the deformability of elements. The results of the new DEM are compared with those of the refined DEM to investigate
the need to consider element deformability.
2. Refined DEM
2.1 Basic concept
Like the original DEM, the refined DEM models astructure as an assembly of rigid elements. The surface of an element is divided into many segments(Fig. 1 (a)), and a spring and a dashpot are attached toeach segment (Fig. 1 (b)). Therefore, the elements
have multiple springs on the surface. This segmentation of the surface allows the spring constant to be defined theoretically, which is the main improvement
over earlier methods. For more details, see [7]. A restoring spring, which models the elasticity of the elements, is set between continuous elements (Fig. 1 (c)).
Structural failure is modeled as breakage of the restoring spring, at which time the restoring spring is
replaced with a contact spring and dashpot (Fig. 1(d)). If two elements are initially continuous, the restoring spring is set between the elements until failure
occurs. If these two elements recontact each otheragain after the failure, a contact spring and dashpot
are set between the elements.
A contact spring and dashpot are used for simulating the contact, separation, and recontact between
elements. The contact dashpots, which are set in parallel with the contact springs, are introduced to express energy dissipation due to the contact. Therefore,
in the refined DEM, the spring constant can be derived theoretically, which makes the method suitable
for simulating large displacement behaviors such as
failure and collapse.
2.2 Analytical parameters(1) Spring constant of each spring
A spring is attached to each surface segment. There
are two types of springs: restoring springs and contact
springs. Each spring constant is taken to have thesame value. The springs are set for both the normal
and shear direction of the surface. The spring constants per area, in both the normal and shear directions, are written as follows using the equilibrium offorce for each pyramid (Fig. 1 (e)) composed of a
44

SEISMIC BEHAVIOR ANALYSIS OF MASONRY STRUCTURES USING THE DISTINCT ELEMENT METHOD CONSIDERING ELEMENT DEFORMABILITY

(a) Segments(b) Multiple springs(c) Restoring spring(d) Contact spring and contact dashpot(e) Pyramid for deriving the spring constantFigure 1. Modeling of a spring and dashpot in the refined DEM%ULFN

0RUWDU

-RLQW

8QLW

(a) Masonry structure (b) Analytical model
Figure 2. Modeling of a masonry structure

segment and the center of gravity for each element:
kn =

E
2

(1 − ν )" , andks =E
2(1 + ν )" ,Young’s moduli are given by EA and EB, and Poisson’s
ratios by νA and νB for elements A and B. The distancesfrom each center of gravity to each segment are lA andlB. For the mortar, Young’s modulus and Poisson’s ratio are EM and νM, respectively, and tM is the mortar
thickness. Assuming that the springs of elements A
and B are connected in series, the spring constants between elements per area, kn and ks, are found to be(1)

where E is Young’s modulus and ν is Poisson’s ratio
for the element; l is the distance from the surface tothe center of gravity.

1" A − tM / 2tM" −t /2++ B M 22
2
E A /(1 − ν A ) E M /(1 − ν M ) E B /(1 − ν B ) , and1
ks =" A − tM / 2tM" B − tM / 2++
E A / 2(1 + ν A ) E M / 2(1 + ν M ) E B / 2(1 + ν B ) . (2)(2)(2) Spring constant between elements
In the refined DEM, individual components of themasonry structure shown in Fig. 2 (a) are modeled ina simple manner, as shown in Fig. 2 (b). The bricks
themselves are modeled with rigid elements whereas
the mortar joints between elements are modeled withmultiple springs and dashpots. The size of one unit is
the sum of the brick size and the mortar thickness.
Each unit is modeled with one rigid element. The
joints have zero thickness and are expressed bysprings between the elements.
Two elements, A and B, are assumed to be continuous and connected to each other with mortar.

kn =In the case where elements A and B are continuous and without mortar, the spring constants are written as
kn =

1
"A
E A /(1 − ν

45

2
A)

+

"B
E B /(1 − ν B2 ) ,and

A. FURUKAWA, R. HORIKAWA, J. KIYONO, K. TOKI

ks =

1
"A
"B
+
E A / 2(1 + ν A ) E B / 2(1 + ν B ) .f 2 (σ ) =| Ĳ | +σ tan φ − c ,where c is the bond strength andφis the friction angle.

(3)(3) Contact damping coefficient
If a segment of an element is in contact or in collision
with another element with which the segment is notcontinuous via the restoring spring, a contact springand dashpot are set between the elements. The contact
dashpot is introduced to express the energy dissipation of the contact. The damping coefficient per area
is calculated as follows:c n = 2hn m ave k n

and(3) Compression failure mode
For the compression failure mode, an ellipsoid cap
model is used. The yield function has the following
form that makes use of compressive strength fm andmaterial model parameter Cs (Fig. 3):2

f 3 (σ ) = σ 2 + C s | Ĳ | 2 − f m .

, c s = 2hs m ave k s ,m ave = ρ A " A + ρ B " B ,(8)

Cs is based on past research on masonry structures [8].
(4)Shear stress ττwhere hn and hs are the damping constants for the normal and shear directions, mave is the equivalent mass
per area relevant to this contact, and ρA and ρB are themass densities of elements A and B.

Shearμ=WDQφ
failureφφCompression failure

Fc

Tensileσ failureσ

IWft Normal stressIPfm

2.3 Modeling of failure phenomena
Let σ and τ be the normal stress vector and shear
stress vector acting at the contact point, and let un andus be the relative displacement vectors in the normal
and shear directions between the segments. Stresses σand τ are expressed as
σ = k n | u n | and Ĳ = k s u s .

(7)

Figure 3. Modeling of failure phenomena
2.4 Nonlinearity of contact force
When a segment of an element is not continuous, but
in contact (or recontact) with another element, a contact spring and contact dashpot are attached to eachsegment. The contact force is generated only while
there is active compression force. The contact force inthe shear direction is bounded by the friction limit.
The friction limit is expressed by(5)

When the stresses acting on the spring reach
their elastic limit, failure phenomena are expressedusing the breakage of the restoring spring. The elastic
limits are modeled using the criteria of tension, shear,and compression failure.

| Ĳ |= σ tan φ ,(9)whereφis the friction angle.
(1) Tension failure mode
When normal stress σ exceeds tensile strength ft, therestoring spring is assumed to be broken by tension
failure. The yield function has the following form
(Fig. 3):f 1 (σ ) = σ − f t .

2.5 Equations of motion(1) Translational motion of center of gravity
The forces acting on the center of gravity of an element are the sum of the restoring and contact forces
between the elements and external forces, such as
gravitational force and inertial force due to an earthquake. The restoring or contact force by the spring
between two certain elements in contact is obtainedby multiplying the spring constant per area (Eq. (3))
by the relative displacement between the relevant elements and the area of the segment. The contact force
by the dashpot between two certain elements in con-

(6)(2) Shear failure mode
For assessment of shear failure, the Coulomb friction
envelope is used. The yield function has the following
form (Fig. 3):

46

SEISMIC BEHAVIOR ANALYSIS OF MASONRY STRUCTURES USING THE DISTINCT ELEMENT METHOD CONSIDERING ELEMENT DEFORMABILITY

tact is obtained by multiplying the damping coefficient per area (Eq. (4)) by the relative velocity between the relevant elements and the area of thesegment. The total sum of the restoring and contact
force for a certain element is obtained by summing up
the restoring or contact force by springs and dashpots
for all elements with which the element of interest isin contact. For more details, see [7]. The equation of
motion for translational motion is written asmx g (t ) = −mg − mz(t ) +¦ F(t ) ,x p (t ) = x g (t ) + x gp (t )(10)Δt < m / K n

,

.

(14)

From Eqs. (2), (4), and (14), the stability condition for translational motion is determined by [9].
Δt ≤ ρ" 2 (1 − ν 2 ) / E

.

(15)

3.Proposed method (Deformable DEM)
3.1 Basic concepts
This study proposes a new element model to consider
the deformation of an element itself. In the new element model (Fig. 4), an element is divided into twoparts. An inner part (deformation part) considers deformation of the element itself, and an outer part(contact part) deals with contact between elements.
The deformation part is modeled as an elastic FEM
element. The stiffness matrix of an eight-node solid
element is used to obtain the restoring force due todeformation. The contact part is assumed to movealong the deformation part with a constant thickness.
Then, the coordinates of the vertices of an element are
determined by the coordinates of nodes of a FEM element. In the new DEM, the behavior of each node is
tracked step by step. The concept of the refined DEM
is applied so as to consider the contact between elements. Then, a restoring spring is set between continuous elements prior to mortar failure, and a contact
spring and a contact dashpot are attached between elements when they are in contact or recontact. It is desirable that the thickness of the contact part is set tobe as thin as possible, to obtain high accuracy.
In the new DEM, structural deformation is expressed by the sum of the deformation of elements
and the contact between elements. Then, the forces
acting on a node are the sum of external forces (e.g.,gravity and inertia), restoring and contact forces be-(11)where I is the tensor of the inertial moment, r(t) is thevector in the absolute coordinate system between thecenter of gravity and the point at which F(t) acts, andR(t) is the matrix representing the transformationfrom the absolute coordinate system to the inertial coordinate system.
The vector from center of gravity xg(t) to arbitrary point xp(t), xgp(t), is obtained by integrating thevelocity obtained in the following equation:x gp (t ) = (R T (t )Ȧ(t )) × x gp (t )(13)

2.6 Stability condition of solution
In this study, the leap-frog method is applied as atime integration scheme to solve the equation of motion. In the leap-frog method, the stability conditionof the solution is written as(2) Rotational motion around the center of gravity
The angular velocity vector in an inertial coordinate
system is obtained by solving the following Eulerian
equation of motion:¦ R(t )r(t ) × R(t )F(t ) ,

The coordinates of the arbitrary point on the element can be traced by solving Eqs. (11)-(13) step by
step.

where xg(t) is the displacement vector of the center ofgravity of an element at time t, m is the mass of theelement, g is the gravitational acceleration vector, isthe ground acceleration vector, and ΣF(t) is the sumof the restoring and contact force vectors at time t.
From Eq. (10), the acceleration of the center ofgravity of an element is obtained, and the velocity anddisplacement are then obtained by integrating the acceleration. Therefore, the coordinates of the center ofgravity of an element can be traced.

 (t ) + Ȧ(t ) × IȦ(t ) =
IȦ

.

(12)

where ω(t) on the right-hand side is the angular velocity vector in the inertial coordinate system, obtainedby solving Eq. (11).
The coordinates of arbitrary point xp(t) are obtained as follows by using xgp(t) in Eq. (12):47

A. FURUKAWA, R. HORIKAWA, J. KIYONO, K. TOKI

ter of gravity and rotation around the center of gravity
are obtained by solving two equations: Newton’s
equation for translational motion and Euler’s equation
for rotational motion around a center of gravity. Then,the coordinate axes of the element coordinate system
are updated. Reaction force against deformation is notconsidered at this point.

tween elements, and reaction forces that oppose element deformation. The coordinates of each node are
obtained by solving the equation of motion for eachnode. Therefore, the behavior of the structure is traced
by tracking the behaviors of each node, step by step.

Deformation
part

Contact part(2) Motion of the node of deformation
The forces acting on each node are the sum of external forces, restoring and contact forces between elements, and reaction force against deformation. The
coordinates of each node are obtained by solving theequations of motion.
Let xin and min be the relative displacement vector
and the mass of a node. The subscript “in” indicates“inner node” since it is the node of the inner part thatconsiders deformation of the element. The equation of
motion is

NodeFigure 4. New element model
3.2 Analytical parameters(1) Contact part
The analytical parameters in the contact part are thespring constants of the restoring and contact springs
and the damping constant of the contact dashpot.
The values used for the restoring and contact
springs are the same as for the refined DEM. The
spring constants of the spring between elements, A
and B, are written as follows by referring to Eq. (2):kn =min x in (t ) = − min g − min z(t ) + ¦ f contact (t ) + ¦ f deform (t ),(18)

where fcontact(t) is the sum of restoring and contact
forces acting on a node at time t, and fdeform(t) is the reaction forces against deformation acting on a node.
The acceleration of a node is determined by solving the equation of motion of a node. The coordinatesof a node at the next time step are obtained by integrating the acceleration. The coordinates of a vertex
are updated by considering the thickness of the part incontact. Then, the forces acting on a node at the nexttime step are also updated. In this way, the behaviorof a structure can be traced by calculating the equation of motion step by step.

1tM

r −t / 2+ B M
E M /(1 −ν M2 ) E B /(1 −ν B2 ) , and1
ks =rA − t M / 2tM
r −t / 2++ B M
E A / 2(1 + ν A ) E M / 2(1 +ν M ) E B / 2(1 +ν B ) , (16)rA − t M / 2

E A /(1 −ν A2 )+where rA and rB are the thicknesses of the contact partsin elements A and B.
The damping constant is determined by substituting k and k into Eq. (4).
ns(2) Deformation partA stiffness matrix for the 3D FEM is introduced in the
deformation part. This study applies an eight-node
hexahedral element. Letting K be the element stiffness matrix, u(t) be the relative displacement vector
at time t, and freact(t) be the vector of reaction force
against deformation, freact(t) is then written as
f react (t ) = Ku .

4.Numerical analysis
4.1 Verification of the proposed method(1) Basic conceptThis section investigates whether the proposed method can adequately model the deformation of individual elements. A vertical load is applied to an element,
and the horizontal displacement is computed to confirm the appropriate simulation of Poisson’s effect.
The results are compared with the results obtained using the FEM.

(17)

3.3 Equations of motion(1) Translational motion of the center of gravity androtational motion around the center of gravity
As with the refined DEM, the coordinates of the cen48

SEISMIC BEHAVIOR ANALYSIS OF MASONRY STRUCTURES USING THE DISTINCT ELEMENT METHOD CONSIDERING ELEMENT DEFORMABILITY(2) Analytical model
Figure 5 illustrates analytical models of the proposed
method and the FEM. In the proposed method, a highdensity rigid element (weight W) was placed on a deformable element to apply the vertical load. The element is a cube measuring 0.1 m on each side and thethickness of the contact part of the element is
0.001 m.
In the FEM, a cube also measuring 0.1 m oneach side was used. A load of W/4 was applied at eachnode on the upper surface of the element.
The displacement of point A (Fig. 5) in the x-direction was computed for both cases to compare thecalculated displacement due to Poisson’s effect.

ement. Table 2 compares the results using the proposed method and the FEM.
In Table 2, the results obtained using the proposed method are found to agree well with the resultsof the FEM. Figure 6 shows the deformation of theelement and the displacement in the x-direction.
It is confirmed that Poisson’s effect, which cannot be modeled by the existing DEM, was successfully expressed by the proposed method. It is also confirmed that the deformation computed by the proposed
method agreed well with the FEM. In masonry buildings, bricks support the vertical gravity load of thebricks stacked above, and this vertical gravity loadcauses deformation in the horizontal direction due toPoisson’s effect. In the following sections, the effectof this Poisson’s effect is examined by comparing theresults of the proposed method and the refined DEM.

0.1m
0.1m
0.1m

0.1m

z

Highdensity
element

Deformableelement

W/4

W/4
W/4

W/4

FEM
element

Table 2 Analysis results
Proposed method
FEM
Case 1 ( v = 0.25) 6.804 × 10-6 (m) 6.802 × 10-6 (m)Case 2 ( v = 0.40) 9.280 × 10-6 (m) 9.142 × 10-6 (m)Point Ayx(a) Proposed method(b) FEM

Figure 5. Analytical model for verificationz(3) Material properties
Table 1 presents the material properties used in thesimulations. The density and Young’s modulus werechosen on the basis of experimental results of adobe
bricks obtained by Ghannad et al. [10] The simulationin this analysis is elastic, so the strengths of mortar(i.e., the tensile, shear, and compressive strengths)
were set to have a value that was large enough not to
fail. Poisson’s ratio in Case 1 is 0.25, and in Case 2, itis 0.40.

   x

Figure 6. Result of this analysis (Poisson’s ratio 0.40)

4.2 Seismic behavior analysis of a masonry wall
(1) Basic concept
Seismic analyses using the proposed method and therefined DEM were conducted by modeling the propagation of a ground motion source into a simple masonry wall. From these analyses, the effect of considering deformation is investigated.

Table 1. Material properties
Brick Mortar Rigid element
Density p (kg/m3) 1.8 × 103 1.8 × 103 2.0 × 104
Young’s modulus E
9.8 × 107
(N/m2)
Poisson’s ratio v
0.25 (Case 1), 0.40 (Case 2)(2) Analytical model
Figure 7 illustrates the masonry wall models analyzed
in this section. The wall has a width of 1.2 m in the xdirection, a depth of 0.2 m in the y-direction, and aheight of 1.2 m in the z-direction. The x-direction isregarded as the in-plate direction, and the y-direction
is regarded as the out-of-plate direction.
Two types of models were applied in these anal-(4) Results
The displacement of point A (Fig. 5) was measured
when the load in the vertical direction acted on the el-49

A. FURUKAWA, R. HORIKAWA, J. KIYONO, K. TOKI

yses. Figure 7 (a) is a straight joint wall model consisting of straight joints in both the horizontal andvertical directions. Figure 7 (b) is referred to as an
English bond wall model. The joint is straight in thehorizontal direction, but not in the vertical direction.
The straight joint wall model consists of cubic elements with edges of 0.2 m long. The English bondwall model consists of two types of elements: cubic
elements with edges of 0.2 m long, and rectangularparallelepiped elements of 0.1 m×0.2 m×0.2 m.
1.2m(3) Input ground motion
The Building and Housing Research Center in Iranrecorded the unique main shock of the Bam earthquake [12]. The acceleration of the horizontal component of the earthquake is oriented in the in-plate direction (Fig. 8).

0.2m

1.2m

z

Figure 8. Input ground motion


y z
x

Element Ax(a) Straight joint wall

Element B

4.3 Results
The relative displacements of the center of gravity ofelement A and B (Fig. 7) are shown in Figs. 9 and 10.
The relative displacements from 0.0 to 1.5 s are used
to compare the results of the refined DEM and the
proposed method. In addition, the seismic behaviors
of structures obtained by the two methods are compared. The results of the refined DEM are referred toas “RIGID” and the results of the proposed method
are referred as “DEFORM”.

  (b) English bond wall

Figure 7. Analytical model(1) History of the relative displacement
Figures 9 and 10 show the relative displacementsfrom 0.0 to 1.5 s. Before initializing ground motion, aself-weight was applied to the structures. The groundmotion was instigated after the structures reached astable condition under their own weight.

(2) Material properties
Table 3 shows the material properties of the bricks
and mortar. The strengths of the mortar are based onfield experimental results obtained during a survey ofthe damage resulting from the 2003 Bam earthquake[11].
Table 3. Material properties
Brick
Density p (kg/m3)
1.8 × 103
Young’s modulus E (N/m2) 9.8 × 106
Poisson’s ratio v
0.25
Tensile strength (N/m2)
2
Shear strength (N/m )Friction angle
Compressive strength (N/m2)-a) Straight joint wall model
Figures 9 (a), (b), and (c) illustrate the relative displacements in the x-, y-, and z-directions, respectively.
From Fig. 9 (a), the amplitude obtained using the refined DEM (RIGID) is nearly 0.0 m and the amplitude obtained using the proposed method (DEFORM)
is about-0.0002 m from 0.0 to 0.3 s. These analyses were conducted after initializing the self-weight, so that theinitial difference between both methods is due toPoisson’s effect under this self-weight. In Fig. 9 (b),Mortar
1.8 × 103
9.8 × 106
0.25
4.6 × 103
2.9 × 10332°
4.9 × 10550

SEISMIC BEHAVIOR ANALYSIS OF MASONRY STRUCTURES USING THE DISTINCT ELEMENT METHOD CONSIDERING ELEMENT DEFORMABILITY

the displacement is observed to be nearly 0.0 s in
RIGID because the refined DEM does not consider
Poisson’s effect. In comparison, Poisson’s effect isobserved using the DEFORM method. The differences in Fig. 9 (a), (b), and (c) are due to the deformation
caused by Poisson’s effect.

mation of elements is ignored. As a result, the modeled structure collapses easily. On the other hand, in
DEFORM, compressive forces do act between elements because of Poisson’s effect. These compressive
forces work as a prestress, so that the structure isstronger and collapses later than in RIGID.

b) English bond wall model
Figure 10 shows the results of an English bond wallmodel. Like the straight joint wall model, the difference from 0.0 to 0.3 s is due to self-weighted Poisson’s effect. As shown in Fig. 10 (b), the displacement in RIGID is 0.0 m, whereas Poisson’s effect isobserved in DEFORM. Therefore, by investigating
the differences in the two methods, the proposed
method is confirmed to express Poisson’s effect byconsidering the deformation of the element.

b) English bond wall model
As shown in Fig. 12, a crack is generated with a diagonal orientation but the structure does not collapse in
RIGID. This is because the joints of an English bondwall are complicated and they enhance the interlocking between elements. In DEFORM, a large crack isgenerated at (a) and the upper part collapses at (b).
The elements deform in DEFORM, so that the interlocking of joints is weaker than in RIGID.
Compared with other models, the English bondwall model is strongest in RIGID. The straight joint
wall model collapses easily due to the straight joint.
However, the English bond wall model has complicated joints and does not readily collapse. In DEFORM, both models collapse, but the number of fallen bricks is smaller in the English bond wall model
than in the straight joint wall model.

(2) Seismic behavior
The seismic behavior of the structure at 4.0, 4.5, 5.0,
5.5, and 6.0 s is compared in this section.
a) Straight joint wall model
As shown in Fig. 11, failure occurs at (a) and thestructure collapses at (c) in RIGID. In DEFORM, thestructure collapses later than in RIGID.
This model has straight joints, so the failure
spreads easily. In RIGID, compressive force due toPoisson’s effect by the self-weight of horizontally adjacent elements cannot be applied because the defor-(a) x-direction

5.Conclusions
In this study, we propose a distinct element methodthat considers element deformability. The deformationof an element can be expressed by dividing the element into a deformation part and a contact part.

(b) y-direction(c) z-direction

Figure 9. History of relative displacement (straight joint wall model)(a) x-direction(b) y-direction(c) z-direction

Figure 10. History of relative displacement (English bond wall model)
51

A. FURUKAWA, R. HORIKAWA, J. KIYONO, K. TOKI(a) 4.0 s(b) 4.5 s (d) 5.5 s(c) 5.0 s(e) 6.0 s

Figure 11. Deformation form (straight joint wall model)(a) 4.0 s(b) 4.5 s(d) 5.5 s(c) 5.0 s(e) 6.0 s

Figure 12. Deformation form (English bond wall model)
Horizontal displacement against the vertical load
is calculated using the proposed method. The results
obtained using the proposed method agree well withthe results of the FEM. The proposed method is confirmed to correctly express Poisson’s effect.
The seismic behaviors of simple masonry walls
were simulated using the proposed method and the refined DEM. The effect of considering deformation
was investigated by comparing the results of bothmethods. Ground motion was input in the in-plane direction. In the refined DEM, compression force due toself-weight between horizontally adjacent elements
cannot be expressed, which means that a straight joint
wall model collapses easily. In the proposed method,compressive forces are generated between horizontally adjacent elements because the method can consider
Poisson’s effect. The structure collapsed later, which
is modeled better in the proposed method than in therefined DEM owing to these compressive forces. An
English bond wall model collapses less easily than astraight joint model because of the interlocking effectdue to the complexity of the joints. In the refined

DEM, the effect of interlocking is clearly seen; however, in the proposed method, the effect is weaker as aresult of element deformation.
When ground motion was input in the out-ofplane direction, the results of the two methods were
almost the same, because there is only one element in
the out-of-plane direction and the influence of Poisson’s effect is small. In real masonry buildings, thecritical failure pattern is out-of-plane failure. The
seismic behavior of masonry buildings, however, isvery complicated since if one wall vibrates in the outof-plane direction, the orthogonal wall vibrates in thein-plane direction. Therefore, the authors would liketo investigate the influence of element deformability
by analysis through the simulation of buildings in afurther study.
Since the validity of the proposed method has
not been confirmed by the experiment, the authors
would like to compare the analytical results with anexperiment to examine the validity of the proposed
method in a future study.

52

SEISMIC BEHAVIOR ANALYSIS OF MASONRY STRUCTURES USING THE DISTINCT ELEMENT METHOD CONSIDERING ELEMENT DEFORMABILITY

Acknowledgement

tin of Earthquake Engineering, Vol. 10, pp. 269-283,
2012.

This work was supported by JSPS KAKENHI
Grant Number 00380585.

[7] A. Furukawa, J. Kiyono, and K. Toki, Proposal of a numerical simulation method for elastic, failure and collapse behaviors of structures and its application toReferences

seismic response analysis of masonry walls, Journal ofDisaster Research, Vol. 6, No. 1, pp. 51-68, 2011.

[1] OCHA (Office for the Coordination of Humanitarian Af-[8] P. B. Lourenco, Analysis of masonry structures with in-fairs), http://www.unocha.org/ [Accessed: Dec, 2012].

terface elements, theory and applications, Delft Uni-

[2] O. C. Zienkiewicz and R. L. Taylor, The finite element

versity of Technology, Faculty of Civil Engineering,th

TU-DELFT report no. 03-21-22-0-01, 1994.

method, 5 edition, Vol. 1, 2, 3, Butterworth Heine-[9] U. M. Ascher and L. R. Petzold, Computer methods formenn: Oxford, U.K., 2000.
[3] K. Oguni, M. Hori, and K. Sakaguchi, Proposal of new

ordinary differential equations and differential-algebra-

FEM for analysis of failure phenomena, Journal of theic equations, Society for Industrial Mathematics, 1998.
[10] M. A. Ghannad, A. Bakhshi, S. E. Mousavi Eshkiki, A.

JSCE, No. 766/I-68, pp. 203-217, 2004.
[4] P. A. Cundall and O. D. L. Strack, A discrete numerical

Khosravifar, Y. Bozorgnia, and A. A. Taheri Behbah-model for granular assemblies, Geotechnique, 29, pp.

ani, A study on seismic vulnerability of rural houses in

47-65, 1979.

Iran, Proc. of the First European Conf. on Earthquake
Engineering and Seismology, Paper No. 680, 2006.

[5] K. Meguro and M. Hakuno, Fracture analysis of concrete structures by the modified distinct element meth-[11] J. Kiyono and A. Kalantari, Collapse mechanism ofod, Proc. of JSCE, Structural Eng./Earthquake Eng.,

adobe and masonry structures during the 2003 IranVol. 6, No. 2, pp. 283s-294s,1989.

Bam earthquake, Bulletin of Earthquake Research Institute, University of Tokyo, 79, pp. 157-161, 2004.

[6] T. Nakagawa, T. Narafu, H. Imai, T. Hanazato, Q. Ali,and C. Minowa, Collapse behavior of a brick masonry[12] BHRC, The very urgent preliminary report on Bamhouse using a shaking table and numerical simulation

earthquake of Dec. 26, 2003. http://www.bhrc.gov.ir,

based on the extended distinct element method, Bulle-

2003 [Accessed: Dec. 2012].

53