Journal of Natural Disaster Science, Volume 33, Number 1, 2012, pp11-36

Numerical Simulation of the Failure Propagation ofMasonry Buildings during an Earthquake
Aiko FURUKAWA*
Junji KIYONO**
Kenzo TOKI****Assistant Professor, Civil and Earth Resources Engineering, Kyoto University**Professor, Department of Urban Management, Kyoto University***Professor, Disaster Mitigation of Urban Cultural Heritage, Ritsumeikan University(Received December 2, 2010 Accepted October 4, 2011)
ABSTRACTThe process of failure propagation of masonry buildings during earthquakes is simulatedusing a refined version of the distinct element method that simulates three-dimensionalelastic, failure, and collapse behaviors of structures. Models with a flat roof and models
with a vault roof are considered, and their failure propagation mechanisms are examined.
The influence of the direction of the input ground motion on failure propagation is also
investigated. Moreover, the effectiveness of three reinforcement measures is compared.
One measure is increasing the mortar strength, the second is increasing the thickness ofthe bearing walls, and the third is introducing wooden columns and beams. Among thesemeasures, increasing the mortar strength is found to be the most effective. Increasing the
thickness of the bearing walls and introducing wooden columns and beams are found to beeffective only if the roof has sufficient integrity since the collapse of the roof depends on
the integrity of the bearing walls and the roof itself.
Keyword: failure propagation process; masonry building; earthquake; distinct element
method

1. Introduction

posed by the European Seismological Commission in
1998, which was modified from the MSK scale (1964)
to be applicable to various modern structures [3].
Similar to the MSK scale, EMS98 categorizes the vulnerability classes of buildings from A to F, and adobebuildings are classified as vulnerability class A, whichis the weakest type of building. EMS98 also classifies
damage to buildings into five grades from G1 to G5(G1: negligible to slight damage, G2: moderate damage, G3: substantial to heavy damage, G4: very heavy
damage, G5: destruction). The relationship between
the damage grade of adobe buildings and the EMS98

The principal cause of death due to earthquakes is thecollapse of buildings, which has accounted for 75% ofearthquake fatalities over the last century [1]. In particular, the collapse of masonry buildings has killed
a large number of victims [2]. An adobe building is
a masonry structure made of sun-dried bricks (adobe
bricks) and mortar and has particularly low earthquake resistance. Such structures collapse even at low
intensities of ground motion and collapse rapidly athigh intensities. EMS98 is a macroseismic scale pro11

A. FURUKAWA, J. KIYONO, K. TOKI

intensity is defined as follows. At intensity V, a few(0-20%) buildings have grade G1 damage. At intensity VI, a few adobe buildings have grade G2 damage
and many (20-60%) have grade G1 damage. At intensity VII, a few adobe buildings have grade G4 damage
and many have grade G3 damage. At intensity VIII, a
few adobe buildings have grade G5 damage and many
have grade G4 damage. At intensity IX, the damage
to many adobe buildings is grade G5, and at intensity
X, the damage to most (60-100%) adobe buildings isgrade G5.

A method based on discontinuum modeling is moresuitable for analyzing failure and collapse phenomena. Examples of numerical methods for a discontinuum are the distinct element method (DEM) [10],discontinuous deformation analysis (DDA) [11], andthe rigid-body spring method (RBSM) [12]. However,these methods have disadvantages. In the DEM, a
way of determining the spring constant from the material properties has not been established, and the values need to be quantified experimentally. Therefore,
the reliability of the results is not high. For example,
Alexandris et al. [13] adopted the DEM to investigate
the crack patterns and collapse mechanism of stone
masonry structures subjected to severe earthquake excitations. They calibrated the analytical parameters bycomparing the numerical prediction of the responsewith experimental results. The calibrated parameters
depend on the scale of the structures, so they do notalways give a reliable solution if the calibrated parameters are applied to models of different scale. The
disadvantage of DDA is that the formulation is complicated since inequality equations need to be solved.
Therefore, application to a three-dimensional problem
is considered difficult even though it is possible intheory. The RBSM is not suitable for large displacement problems since it does not consider the recontact
of an element with neighboring elements other thanthe element with which contact is initially set.

Therefore, it is necessary to improve the earthquake
resistance of these primarily weak masonry buildings to reduce the number of casualties. However,knowledge is still limited as to how the failure begins
and proceeds, how buildings collapse, and how earthquake resistance can effectively be improved. With
this background, a series of seismic behaviors—fromelastic to failure to collapse behaviors—is simulatedusing a refined version of the distinct element method[4].
Among numerical simulation methods, the finite element method (FEM) is the most common for the analysis of a continuum [5]. It can deal with both elastic
and plastic behaviors, but it has difficulty in solving
failure and collapse phenomena since it is based on
the mechanics of the continuum and uses a continuous
shape function. Ghannad et al. [6] conducted seismic
simulation of typical Iranian masonry buildings using
the FEM, and assessed the initiation of cracks fromstress distributions. The crack pattern obtained by the
analysis was in good agreement with actual crack occurrence during the Bam earthquake. Cao et al. [7]computed the earthquake response of adobe structureswith a joint element that expresses the separation andsliding between bricks. Gambarotta and Lagomarsino[8] and Lourenco et al. [9] introduced joint elements
at the interface as well as inside each brick, and simulated cracks inside the bricks. They confirmed thevalidity of their simulation by comparing the simulation results with the experimental results. In this way,
the FEM can simulate discontinuous phenomena suchas cracking, separation, and sliding through the use of
joint elements. It is still difficult, however, to simulate
total collapse behavior employing the FEM.

As an alternative, the present paper uses a refinedversion of distinct element method, which simulates
three-dimensional elastic, failure, and collapse behaviors [4]. The method models the structure as an assembly of rigid elements and cannot handle Poisson’s
effect like the DEM, but unlike the DEM, interaction
between elements is modeled by multiple springs andmultiple dashpots attached to the surfaces of the elements. The spring constant of each spring is theoretically determinable.
Using this method, this study simulates the seismic
behavior of simplified adobe buildings in Iran. Models
with a flat roof and models with a vault roof are considered, and their failure propagation mechanisms are
examined and compared. The influence of the directionof the input ground motion on the process of failure
propagation is also investigated. The outcomes of thesimulation are compared with examples from past
12

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

earthquakes to confirm the reliability of the simulation.
Moreover, this study investigates the effectiveness ofreinforcement measures. Through a damage survey, it
is pointed out that poor-quality mortar, weak wall-towall connection, and lack of out-of-plane resistance ofwalls are important issues in considering the seismic
vulnerability of rural houses in Iran [14]. Therefore,
this study considers three reinforcement measures.
One is increasing the mortar strength, which increases
brick-to-brick connections and wall-to-wall connections as a result. The second is increasing the thicknessof the bearing walls, which increases their out-of-plane
resistance. The third is introducing wooden columns
and beams, which increases the out-of-plane resistanceof walls and wall-to-wall connections. The effects ofroof type, direction of ground motion, and reinforcement measures on failure propagation are investigatedin terms of tensile cracking, building deformation, andscattering of bricks determined by simulation.

The surface of an element is divided into small segments as shown in Fig. 1 (c). The segment in the figure is rectangular, but the method does not limit thegeometry of the segment. The black points indicate
the representative point of each segment, and the relative displacement or contact displacement between
elements is computed for these points. Such points
are referred to as contact points or master points inthis study. One restoring spring and one combinationof contact spring and dashpot are attached to one segment (Fig. 1 (d) ) at each of the representative points
in Fig. 1 (c). The spring constant for each segment isderived on the basis of the stress-strain relationship ofthe material and the segment area.
Forces acting on each element are obtained by summing the restoring force, contact force, and other external forces such as the gravitational force and inertial
force of an earthquake. The behavior of an element
consists of the translational behavior of the center ofgravity and the rotational behavior around the centerof gravity. The translational and rotational behaviors ofeach element are computed explicitly by solving Newton’s law of motion and Euler’s equation of motion.

2. Analytical Method
2.1 Basic conceptThis study employs a refined version of the DEM tosimulate a series of structural dynamic behaviors fromelastic to failure to collapse phenomena. A structure ismodeled as an assembly of rigid elements, and interaction between the elements is modeled with multiple
springs and multiple dashpots that are attached to thesurfaces of the elements. The elements are rigid, but
the method allows the simulation of structural deformation by permitting penetration between elements.

(a) Restoring spring

Figure 1 (a) shows a spring for computing the restoring
force (restoring spring), which models the elasticity ofelements. The restoring spring is set between continuous elements. Structural failure is modeled as breakage of the restoring spring, at which time the restoring
spring is replaced with a contact spring and a contact
dashpot (Fig. 1 (b) ). Figure 1 (b) shows the spring anddashpot for computing the contact force (contact springand dashpot) and modeling the contact, separation, andrecontact between elements. The dashpots are introduced to express energy dissipation due to the contact.
Structural collapse behavior is obtained using these
springs and dashpots. The elements shown in Figs. 1 (a)and (b) are rectangular parallelepipeds, but the method
does not limit the geometry of the elements.

(b) Contact spring and dashpot(c) Segments and contact points(d) Multiple springs and multiple dashpots
Fig. 1 Basic concept of the analysis method
13

A. FURUKAWA, J. KIYONO, K. TOKI

2.2 Spring constant of each element
As described in Section 2.1, there are two types ofsprings, namely restoring and contact springs. It is assumed that the spring constants of the restoring spring
and those of the contact springs are the same. It isconsidered that each segment has its own spring.

kn =E A /(1 −ν A2 )ks =E
(1 −ν )"
2

ks =

,

kn =

E
,
2(1 + ν )"

ks =

(1)
where E is Young’s modulus, v is Poisson’s ratio, and" is the distance from the surface at which the springis connected to the center of gravity.

E B /(1 −ν B2 )

1
.
"A
"B
+
E A / 2(1 + ν A ) E B / 2(1 + ν B )1," A − tM / 2tM" B − tM / 2++
E A /(1 − ν A2 ) E M /(1 − ν M2 ) E B /(1 − ν B2 )1,
" A − tM / 2" B − tM / 2tM++
E A / 2(1 + ν A ) E M / 2(1 + ν M ) E B / 2(1 + ν B )(3)
where tM is the mortar thickness, EM is Young’s modulus, and vM is Poisson’s ratio of the mortar. The normal direction of forces is the direction perpendicular
to the surface of the master point of element A.

2.3 Modeling of elastic behavior
It is assumed that two elements, A and B, are continuous, and that a contact point of element A is continuous with element B as shown in Fig. 2 (a). The contact point of element A is the “master point” and thecontact point of element B is the “slave point.” Theslave point does not have to be the center of the segment. The combinations of continuous elements foreach segment are determined in the initial stage of thecomputation.

Let σ and τ be the normal and shear stresses acting
at the contact point, and let un and us be the relative
displacements between the adjacent master and slave
points in the normal and shear directions. The relation
between traction (σ, τ) and relative displacement (un,us) is then written as

σ = knun ,τ = ksus .

(4)
It is noted that the method cannot handle Poisson’s
effect since it considers the contact between two elements.

Let GA and GB be the centers of gravity of elements A
and B, respectively. Let " A be the distance from GA
to the surface of element A in contact. Let " B be thedistance from GB to the surface of element B in contact. Let EA and EB be Young’s moduli and vA and vB
be Poisson’s ratios of elements A and B.

Element A
GA

contact point(master point)A

It is assumed that both elements have their own spring
constants in the normal and shear directions based on
their material properties as shown in Fig. 1 (b). The
spring constants per area for each element are obtained from Eq. (1). Assuming that these springs areconnected in series, the spring constants between elements per area, k and k , are shown in Fig. 2 (b).
n

+

,

"B

(2)
The spring constant between elements connected bymortar is also obtained in a similar manner. For example, in masonry structures, bricks are often connectedwith mortar. In this case, the spring constant per area
between elements (bricks) is obtained as

Springs are set for both the normal and shear (tangential) directions of the surface. Let us denote thearea of the segment as dA and the relative (contact)displacement at the surface segment as un and us. The
subscripts n an s indicate the values in the normal andshear directions, respectively. The spring constants
per area in the normal and shear directions, kn and ks,
are obtained as follows [4].

kn =

1
"A

contact point(slave point)

B

GB
Element B

Element A
GA

A

contact point(master point)

contact point(slave point)
B

GB
Element B(a) Continuous elements (b) Restoring spring between elements
Fig. 2 Modeling of the restoring spring between elementss

14

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

interface between bricks, and we do not discriminate
failure at the mortar from failure at the interface between brick and mortar. The restoring spring between
elements is linear until any failure occurs. In tension,
linear tensile reaction force acts until it reaches thetensile strength. After tensile failure, no tensile reaction force acts between these elements until they arein contact again. In compression, the linear compressive reaction force acts until it reaches the compression strength, and the compressive reaction force
then obeys the Mohr-Coulomb cap model. During theunloading process after the failure, the compressive
strength decreases linearly with the initial stiffness,and it becomes 0 when elements are separated. The
shear reaction force is also linear until it reaches theshear strength, and it then obeys the Mohr-Coulomb
cap model. During the unloading process after the
failure, the shear strength decreases linearly with theinitial stiffness and obeys the limit of Eq. (6). Since
the face of the elements is divided into many segments and each segment has its own spring, these
springs gradually break during loading.

2.4 Modeling of failure phenomena
The elastic behavior of structures is demonstrated by
the linear multiple restoring springs between continuous elements until the restoring force of a spring
reaches its elastic limit. The elastic limits are modeledusing the criteria of tension, shear, and compression
failure. When a spring reaches one of these limits, it
is judged that failure has occurred at that segment of
the spring. After the failure, the restoring spring is
replaced with a contact spring and dashpot at this segment. The method can trace the expansion of failure
between elements. The three failure modes—namely,
tension, shear, and compression failure modes—aredefined based on the Mohr-Coulomb cap model asfollows.
(1) Tension failure mode
For the tension failure mode, the parameter considered is tensile strength ft. When the normal stress of
spring σ exceeds the tensile strength, the restoring
spring is assumed to be broken by the tension failure.
The yield function has the following form (Fig. 3).
f1 (σ ) = σ − f t .

(5)
The normal restoring stress cannot exceed this limit.
(2) Shear failure mode
For the shear failure mode, the Coulomb friction envelope is used. The parameters considered are bondstrength c and friction angle ϕ. The yield function has
the following form (Fig. 3).
f 2 (σ ) =| τ | +σ tan φ − c .

Shear stress τ
Shearμ=WDQφ
failureφφCompression failure
IP
fm

(6)

Fc

Tensileσ failureσ

IWft Normal stressFig. 3 Modeling of failure phenomena

The shear restoring stress cannot exceed this limit.
(3) Compression failure mode
For the compression mode, an ellipsoid cap model isused. The yield function has the following form (Fig. 3).
(7)
where fm is the compressive strength and Cs is thematerial model parameter. Cs = 9 is adopted on thebasis of past research [15]. When the restoring stressexceeds this limit, both the normal and shear restoring
stresses are reduced in the same proportion to meetthis limit.

2.5 Modeling of contact and recontact between
elements
If a segment of an element is in contact with anotherelement with which the segment is not continuous via
the restoring spring, the contact spring and dashpot
generate contact force between the elements. Contact
between a segment and the surface of another element
is detected at each time step for all segments that are
not continuous with other elements via a restoring
spring.

In real structures, failure occurs not only at the mortar, but also inside bricks and at the interface between
brick and mortar. However, in this study, it is assumed
that failure does not occur inside bricks but at the

Figure 4 (a) illustrates the situation where a contact
point of element A is in contact with discontinuous element B, and Fig. 4 (b) shows the contact spring anddashpots between these elements.

f 3 (σ ) = σ 2 + C sτ 2 − f m ,215

A. FURUKAWA, J. KIYONO, K. TOKI

and contact forces acting at each point are obtained byconsidering the area of the master point. Assuming dA
is the area of the master point, the effective area canbe half the segment area dA to avoid double countingof the contact force, since the contact points of bothelements A and B can be the master points. Therefore,
assuming that dA/2 is the effective area of the masterpoint, the spring constants and damping coefficientsare

The spring constant and the contact forces in the normal and shear directions are calculated in the samemanner as for the restoring force. The differencesfrom the case for the restoring force are that the contact force is generated only while the compression
force acts and that the shear force is bounded by thefriction limit.

τ = σ tan φ ,(8)where ϕ is the friction angle.
The 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 ,c s = 2hs m ave k s ,GAcontact point(master point)

GB

B
Element B

d s = C s Δu s / Δt ,2.7 Structural damping
In this study, structural damping that expresses energy
dissipation in structural vibration is considered. Thisdiffers from energy dissipation of the contact expressed by damping constants cn and cs.
Damping constant c is introduced as mass-proportional damping as
c = αm ,(13)
where α is a variable defining the damping constant. α
can be calculated as
α = 2hi ω i ,(14)
where ωi is the angular frequency and hi is the corresponding damping constant.
Mass-proportional damping is adopted since the method solves the equation of motion for each element
and mass-proportional damping is easily applicable.
The dashpot attached between contact elements is a
stiffness-proportional damping. We plan to investigate
an appropriate model of damping in a future study.

Element A
GA
A

contact point(slave point)A

es = K su s ,(12)where Δun and Δus are the increments of relative displacement and Δt is the time interval. The variables
written in bold face are vectors.

(10)
are the mass densities of ele-

contact point(slave point)

C s = c s dA / 2 .

d n = C n Δu n / Δt ,The damping constants should be evaluated according
to the properties of the elements, but this study adopts
critical damping (hn = hs = 1.0) by considering the factthat most structural components tend not to bounce
greatly and that their oscillation tends to quickly disappear when they collide with each other.

Element A

C n = c n dA / 2 ,e n = K nu n ,m ave = ρ A " A + ρ B " B ,ρB

K s = k s dA / 2 ,(11)
In consideration of this, the size of the segments
should be small and the same for elements A and B.
Finally, the spring forces (en, es) and damping forces(dn, ds) in the normal and shear directions are writtenas(9)
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. In this study, mave iscalculated aswhere ρ A andments A and B.

K n = k n dA / 2 ,contact point(master point)
B

GB
Element B(a) Discontinuous ele- (b) Contact spring andments in contact
dashpot between elements
Fig. 4 Modeling of the contact spring and dashpot between elements
2.6 Restoring and contact forces acting at eachelement
The restoring and contact forces are calculated at all
segments for both elements A and B. The restoring

2.8 Equations of motion
Equations of motion can be constructed using the restoring and contact forces and other external forces.
16

ω

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

the solution is conditionally stable [16]. The following inequality is used to determine the time interval
for the computation.

The motion of each element is obtained by solving thetwo equations of motion. One is the equation for thetranslational motion of the center of gravity, and theother is the equation for the rotational motion aroundthe center of gravity.
(1) Translational motion of the center of gravity
The forces acting on an element are the sum of external forces, such as the gravitational force and inertial
force due to an earthquake, and the restoring and contact forces between elements. The equation of motion
for the translational motion is§E /(1 − ν 2 ) ·¸.
Δt ≤ min¨ " /
¨
¸
ρ
©
¹

(19)
The term on the right hand side of Eq. (19) is obtained
for all segments of all elements, and the time interval
is determined. If the solution is unstable using thistime step, a smaller value is used until stability is attained.

mx g (t ) + cx g (t ) = mg − mz(t ) + ¦ F (t ) ,3. Analysis Outline(15)
where xg(t) is the displacement vector of the center ofgravity of an element at time t, m is the mass of theelement, c is the damping constant of the element, g is
the gravitational acceleration vector, z t is the groundacceleration vector at time t, and ΣF(t) is the sumof the restoring and contact force vectors at time t.
(2) Rotational motion around the center of gravity
First, theω angular velocity vector ω(t) is obtained bysolving theω following Euler equation of motion.

3.1 Modeling of bricks and mortar
In this study, individual components of the masonry
structure shown in Fig. 5 (a) (i.e., brick and mortar
joints) are modeled in a simple manner as shown in
Fig. 5 (b). The bricks are modeled with rigid elements
and the mortar joint between elements is modeledwith multiple springs and multiple dashpots. The sizeof one element is the sum of the brick size and thethickness of mortar together. The multiple springs andmultiple dashpots interact with the surfaces of adjacent elements.

 (t ) + Ȧ (t ) × IȦ (t ) = ¦ R (t )r (t ) × R (t ) F (t ) .
IȦ(16)Here, I is the tensor of the moment of inertia, r(t) isthe vector between the center of gravity and the point
where force F(t) is applied. R(t) is the matrix representing the transformation from the absolute coordinate system to the inertial frame of reference.

The modeling in this study is three-dimensional,
and the elements are modeled with rigid rectangularparallelepipeds and hexahedrons. Faces surrounding
the elements are divided into segments. The interval
between contact points of neighboring segments is set
to a quarter of the shortest edge length according toa past study [4]. The method was verified through a
comparison with the experimental results of monotonic loading and the analytical results of free vibrationby the FEM [4]. Verification by comparison with cyclic loading and shaking table tests is very importantand comprises one of our future plans. The dynamic
failure patterns and scattering of elements, which will
be shown later, are numerical solutions based on theassumptions shown in the previous section. Moreover,
adobe brick is a very brittle material compared to other masonry elements such as burnt brick, stone, andconcrete block, and the adobe bricks themselves fracture during earthquakes. However, due to the limitation of the analytical method, we do not model failure
inside such bricks. Therefore, the solution focuses on
structures where mortar failure is dominant and brick

The vector from center of gravity xg(t) to arbitrary
point xp(t) in an element, xgp(t), is obtained by solving
the following differential equation using the angularvelocity vector.

x gp (t ) = (R T (t )Ȧ(t )) × x gp (t ) ,(17)
where ω(t) is the angular velocity vector in the inertial frame of reference obtained by solving Eq. (16).
The coordinates of point xp(t) are then obtained as follows.

x p (t ) = x g (t ) + x gp (t ) .

(18)
In this study, Eqs. (15) – (18) are solved explicitly using the central difference scheme.
2.9 Critical time interval
Since the equations of motion are solved explicitly,17

A. FURUKAWA, J. KIYONO, K. TOKI

failure is negligible. In this respect, the method may
not be appropriate for structures where failure inside
the bricks occurs with high probability. Modelingof failure of the brick itself is also one of our futureplans.

Mortar joint
Brick(a) Masonry structureelements of 0.2-m length to express their deformationand separation.
The roof consists of 27 beams that run in the y direction. The dimensions of the beams are 0.2 m × 0.2 m
× 3.4 m (length). These roof beams are also divided
into elements of 0.2-m length to express their deformation and separation. The roof beams rest on thewalls and there are mortar joints between the beams
and walls. Wooden beams are also modeled withrigid elements in the same manner as for the bricks.
The conjunction between wooden elements has thestrength of wood.

Spring and dashpot
Element(b) Analytical modelFig. 5 Analytical modeling of masonry structures(2) Vault roof model
The vault roof model has two bearing walls, two nonbearing walls, and a vault roof. Modeling lower than
2.4 m is the same as that for the flat roof model.

3.2 Analytical models
Figure 6 illustrates the two masonry models analyzed
in this study. One has a flat roof made of wooden
beams, and the other has a vault roof made of adobe
bricks. The directions (x, y, z) are also shown in Fig.
6. In each model, the walls facing the y direction support the roof; therefore, these walls are referred to as
bearing walls and the other two walls are referred toas nonbearing walls. The basement is modeled by afixed element (not shown in the figure) and the modelis placed on the basement.
(1) Flat roof model
The flat roof model has two bearing walls, two nonbearing walls, and a flat roof. The external width of
the building is 5.4 m × 3.4 m, and the internal widthis 5.0 m × 3.0 m. The height of the walls is 2.4 m. The
thickness of walls is 0.2 m. The walls are composedof bricks with dimensions of 0.2 m × 0.1 cm × 0.1
m. There are mortar joints with a thickness of 0.01
m between the bricks. Bricks are stacked according
to the English bond model. The connections between
walls and their perpendicular walls are constructed bymortar and bricks stacked with a shift of half the bricklength in alternating layers.

The roof has a semicircular profile with a radius of
1.5 m. The dimensions of the bricks for the walls arebasically the same as those of the flat roof model, but
if the brick is adjoining the roof, it is then cut to makethe wall semicircular. The vault roof has a thicknessof 0.2 m and is divided into bricks having a section
area of 0.2 m × 0.2 m in the radial direction and alength of close to 0.2 m in the tangential direction.
Element 2
Element 1

5.4m

Roof thickness
0.2m

3.4m
0.2m
2.4m

z
y

Wooden
Beam

1.5m
3m

Door 1.2m×2.0m
Window 1.0×1.0mx(a) Flat roof model(b) Vault roof modelFig. 6 Masonry models

3.3 Analytical cases(1) Mortar strength
Kiyono et al. conducted a field experiment to measure
the mortar strength of a damaged adobe building during a reconnaissance survey of the damage resulting
from the Bam earthquake in Iran [17]. The results areshown in Table 1 and are referred to as the results forweak mortar in this study.

One of the bearing walls has a door opening and awindow opening. The size of the door is 1.2 m (width)
× 2.0 m (height). The size of the window is 1.0 m × 1.0
m, and the height of the lower side of the window is 1.0
m. There are wooden beams above the openings that
support the elements above the openings. The dimensions of the beams above the openings are 0.2 m ×
0.1 m × 1.4 m (length). These beams are divided into

Ghannand et al. conducted a laboratory experiment tomeasure the strength of the adobe bricks used in rural
regions in Iran [6]. The results are shown in Table 1
18

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

Table 1 Strength of mortar and wood
VariableWeakmortar [17]Normalmortar [6]Strongmortar [18]WoodTensile strength ft (N/m2)

4.6 × 103

1.27 × 105

1.42 × 106

1.1 × 107

Shear strength c (N/m2)

2.9 × 103

3.04 × 104

1.19 × 106

9.0 × 106

Friction angle ϕ

32°

No data(32° used)

36°

0°

Compressive strength (N/m2)

No data(4.9 × 105 used)

4.9 × 105

2.66 × 107

4.5 × 107

and are referred to as the results for normal mortar inthis study.

(2) Analytical cases
Five models are considered for flat and vault roof
models as presented by Fig. 7 and Table 2.
● Flat20W and Vault20W are models with weakmortar and bearing walls 0.2 m thick.
● Flat20N and Vault20N are models with normal
mortar and bearing walls 0.2 m thick.
● Flat20S and Vault20S are models with strong
mortar and bearing walls 0.2 m thick.
● Flat40N and Vault40N are models with normal
mortar and bearing walls 0.4 m thick.
● Flat20R and Vault20R are models with normal
mortar, bearing walls 0.2 m thick, six wooden
columns supporting each bearing wall, twowooden columns supporting each nonbearing
wall, and four wooden beams around the base of
the roof. It is assumed that the bottom elementsof the wooden columns are fixed to the ground,and that all the wooden members are connected
to other members with mortar.

Since the compression strength was not measured forweak mortar [17], the value for normal mortar is used.
In addition, since the friction angle was not given for
normal mortar [6], the value for weak mortar is used.
The outcome of the simulation does not depend onthese values since all models only suffered tensile
failure in this study as explained later.
For comparison, the mortar strength of bricks used inEurope measured by Marzahn is given in Table 1 andis referred to as the strength of strong mortar in thisstudy [18].
The results show that the mortar used in Europe ismuch stronger than that used in Iran, and that buildings damaged during the Bam earthquake had especially weak mortar strength. It is unrealistic to apply
strong mortar to weak adobe bricks, but this case ischosen to investigate the effect of mortar strength.

Weakmortar

The variation in the mortar strength of adobe buildings may be large, but this study conducted a deterministic analysis without considering variation inproperties. We discuss the effect of mortar strengthby simply comparing the results for three different
mortar strengths. Failure of the brick itself is ignored
considering that the most failure occurs at the mortar
conjunction.

Thin bearing
wall of 20cm

Wooden columns and beams are divided into many
elements to express deformation and breakage.
The breakage is modeled for the mortar listed in Table 1.

Strongmortar(Flat20S)(Flat20W)

Normalmortar(Flat20N)

Thick bearing
wall of 40cm
(Flat40N)

Wooden columns
and beams(Flat20R)(a) Flat roof model

19

A. FURUKAWA, J. KIYONO, K. TOKI

listed in Table 3.
Weakmortar

Structural damping is assumed to be 2% mass-proportional damping. The value of α, which defines the
mass-proportional damping, is obtained so that thedamping constant for the first natural frequency in they direction is 2%. The computation of the natural frequency is described in Section 3.5.

(Vault20W)

NormalmortarThin bearing
wall of 20cm
(Vault 20N)

Strongmortar(Vault20S)

Thick bearing
wall of 40cm
(Vault40N)

3.5 Natural frequency
The history of impact acceleration in the x and y
directions is input separately for investigating thenatural vibration characteristics in each direction. The
impact acceleration is first input into the undamped
masonry models, and the displacement histories of
elements 1 and 2 (Fig. 6) are computed. Element 1 islocated at the center of the nonbearing wall at a heightof 2.0 m, and element 2 is at the center of the bearing
wall at a height of 2.0 m. The time interval of the impact acceleration is 0.01 s, the acceleration at the firsttime step is set as 1 m/s2, and the acceleration for the
rest of the time is set as 0 m/s2.

Wooden columns
and beams(Vault20R)(b) Vault roof modelFig. 7 Analytical cases
3.4 Material properties
Table 3 presents the material properties. The massdensity and Young’s modulus of adobe brick arechosen on the basis of experimental results for adobe
bricks obtained by Ghannand et al. [6]. The mass density and Young’s modulus of the mortar are assumed
to be the same as those of the adobe brick. The material properties for the wooden columns and beams are

Table 4 lists the first natural frequencies in each direction obtained by Fourier transformation of the free-

Table 2 Analytical cases
Model nameRoof type

Thickness ofbearing wall

Mortar strength

Reinforcement measure

Weak

-

20 cm

Normal-StrongIncreasing mortar strength

Flat20W
Flat20N
Flat20S

Flat

Flat40N

40 cm

NormalIncreasing thickness of bearing wall

Flat20R

20 cm

NormalIntroducing wooden columns and beams

Weak

-

Normal-StrongIncreasing mortar strength

Vault20W
Vault20N
Vault20S

Vault

20 cm

Vault40N

40 cm

NormalIncreasing thickness of bearing wallVault20R

20 cm

NormalIntroducing wooden columns and beams

Table 3 Material properties
Variable
Mass density (kg/m3)
Young’s modulus (N/m2)
Poisson’s ratio
Structural damping

Adobe brick
1.8 × 103
9.8 × 107
0.25

MortarWood
1.8 × 103
7.0 × 102
9.8 × 107
6.3 × 1080.25
0.3
mass-proportional damping(2% damping at the first-order natural frequency)
20

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

Table 4 First-order natural frequencies
Model nameRoof typex direction

y direction

Flat20S/Flat20N/Flat20S
Vault20S/Vault20N/Vault20S

Flat

6.64 Hz6.05 Hz

Flat40N
Vault40N
Flat20R
Vault20R

Vault

3.91 Hz3.13 HzFlat

7.03 Hz6.08 Hz

Vault

4.49 Hz3.32 HzFlat

6.84 Hz6.45 Hz

Vault

4.29 Hz3.32 Hzvibration responses. The value in the x direction is
obtained from the response of the center of gravity ofelement 1, and the value in the y direction is obtained
from the response of the center of gravity of element 2.

of 2%. Figure 9 (a) is the response spectrum forthe ground motion from 0 to 3.5 s. Of the horizontal components, it is found that the L component isdominant from 0 to 3.5 s. Figure 9 (b) is the response
spectrum for the ground motion from 0 to 10 s. The L
component is less than the T component for a natural
frequency of less than 6 Hz and is greater than the Tcomponent for a natural frequency greater than 6 Hz.
All vault roof models have natural frequencies of lessthan 6.0 Hz (Table 4), so the T component is dominant for the vault models from 3.5 to 10 s. Although
all flat roof models have natural frequencies higher
than 6.0 Hz, the natural frequencies decrease as thefailure proceeds, so it can be concluded that the Tcomponent is dominant in the failure process after 3.5
s except for the flat models, such as Flat20S, which
do not show remarkable failure behavior. From thesefindings, T is referred to as the latter dominant component and L, the initial dominant component in thisstudy. The gravitational acceleration is considered to

The input acceleration is low enough to prevent
failure of the mortar, so difference due to the mortar
strength is not presented. Models with a flat roof havehigher natural frequencies than those with a vault
roof since wooden beams are stiffer and lighter than
adobe bricks. Models with thick bearing walls havethe highest natural frequencies in the x direction andsecond-highest natural frequency in the y direction.
The models with wooden columns and beams havethe second-highest natural frequency in the x direction
and the highest natural frequency in the y direction.
3.6 Time interval
Substituting the material properties listed in Table 3
and the dimension of each element into Eq. (19), Δt ≤5.03 × 10–5 s is obtained for all cases. Therefore, thetime interval of Δt = 5.0 × 10–5 s is used.

Tbe g = {0.0 0.0 − 9.8} .

Two cases for the input ground motion direction are
considered as presented in Table 5. In the case of
TLV, the T, L, and V components are input in the x, y,
and z directions, respectively. In the case of LTV, the L,T, and V components are input in the x, y, and z directions, respectively.

3.7 Input ground motion
The Building and Housing Research Center in Iran recorded the unique main shock of the Bam earthquake[19]. The maximum accelerations of the two horizontal (T, L) and vertical (V) components of the earthquake after correction were 6.2344 m/s2 (T), 7.7828
m/s2 (L), and 9.7995 m/s2 (V) (Figs. 8 (a) – (c) ). The
direction of the accelerometer for the L component is
N278E. The vertical component is very large. For thehorizontal components, the peak ground accelerationof the L component is greater than that of the T component.

In the analysis, the gravitational force is input first.
After structural stability under gravitational force isconfirmed, the ground motion is input. We confirmedthat all models stand stably under gravitational force
without failure.

Figure 9 presents the displacement response spectraof the three components for a damping coefficient21

A. FURUKAWA, J. KIYONO, K. TOKI(a) T component(b) L component(c) V component

Fig. 8 Time histories of input ground motion

(a) 0–3.5 s(b) 0–10 sFig. 9 Displacement response spectra for a damping coefficient of 2%
Table 5 Cases of the input ground motion direction
Case

x direction

z direction

TLV

T component

L component

V component

LTV

L component

T component

V component

4. Resultsa) TLV (Fig. 10)
Since the mortar is weak, there was tensile failure
throughout the walls. Tensile failure occurred first in
the bearing wall at 1.20 s as shown in the picture andthen in the nonbearing wall at 2.00 s as shown in the
picture. The roof beams on both edges then separated
at 2.80 s owing to the different vibration characteristics of the two walls. There were diagonal cracks
around the openings in the bearing wall at 4.80 s, andthe nonbearing walls then failed from the top andbricks of the nonbearing walls fell. The bearing walls
did not fall, so the roof did not collapse.

4.1 Failure propagation of flat roof models(1) Flat20W
The seismic behavior of the flat roof model with weakmortar (Flat20W) is shown in Figs. 10 and 11 for the
TLV and LTV cases, respectively. Failure propagation is expressed in the following manner. The faceof an element where no failure occurs is shown inlight gray. The face of an element where the mortar
connection is broken by tensile failure is shown in
red and surrounded by black lines. Neither shear nor
compression failure occurred.

Ly direction

T

Fig. 10 Seismic behavior of Flat20W, TLV
22

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

T

L

Fig. 11 Seismic behavior of Flat20W, LTVb) LTV (Fig. 11)
Similar to the case of TLV, there was tensile failure
throughout the walls, first in the bearing wall andthen in the nonbearing wall. Separation of roof beams
at both edges is also seen. Deformation of the bearing wall in the out-of-plane direction is seen at 4.4s, which is not seen in the case of TLV. There werehorizontal cracks in the bearing wall at 4.80 s, andthe roof then collapsed with the failure of the bearing
walls.
c) Effect of the input ground motion direction
In the case of LTV, since the latter dominant component, T, is input in the out-of-plane direction of the
bearing wall, the bearing wall deforms more readily,leading to the total collapse of the roof. The increase in
out-of-plane deformation increases the P–∆ effect of theroof weight, where roof weight P acting on the bearing
wall with deformation ∆ enforces a moment of P∆ on
the bearing wall, and thus, the building collapses.

(2) Flat20N
The seismic behavior of the flat roof model with normal mortar (Flat20N) is shown in Figs. 12 and 13 forthe TLV and LTV cases, respectively. Figures 12 (a)and 13 (a) are the front views and Figs. 12 (b) and 13(b) are the rear views. Elements 1, 2, and 3 are elements whose histories of displacement are presented
later in the paper. Failure propagation is expressed ina manner similar to that for Flat20W. Neither shear
nor compression failure occurred.
a) TLV (Fig. 12)
As seen in Fig. 12 (a), there was no tensile failureof the walls. First, diagonal cracks formed from thecorner of the openings in the bearing wall at 2.00 s asshown in the picture. Horizontal cracks are also seen
between the window and door openings. The numberof diagonal and horizontal cracks then increased in
the bearing wall at 2.80 s. There were vertical cracks
at the edge of the bearing wall with openings at 3.60 s.
These vertical cracks allowed the nonbearing wall tomove more easily, and horizontal cracks formed in thenonbearing wall at 4.80 s. As seen in Fig. 12 (b), there
were no diagonal cracks in the bearing wall without
openings. Horizontal cracks were generated by outof-plane deformation.
b) LTV (Fig. 13)
Similar to the case of TLV, there were diagonal cracks
and horizontal cracks in the bearing wall with openings at 2.00 s as seen in Fig. 13 (a). Vertical cracks inthe bearing wall generated at 2.80 s, earlier than in thecase of TLV. The bearing wall then deformed in theout-of-plane direction at 6.00 s and fell in this direction. The roof collapsed with the failure of the bearing
wall with openings. As seen in Fig. 13 (b), the opposite bearing wall without openings deformed at 6.00 s
because the roof pushed the wall outward, and it lostits balance and fell with the roof.

In the case of TLV, the latter dominant component, T,
is input in the in-plane direction of the bearing wall,and clear diagonal cracks are seen in the bearing wall.
However, deformation of the bearing walls in the outof-plane direction is smaller and the walls do not fall,
which prevents the total collapse of the roof.
The initial dominant component, L, affects the initial
failure propagation, but the time in which the L component is dominant is short. Therefore, the collapse
process is dominated by the latter and longer dominant component, T.
It is thus found that the total collapse of a building
depends on the integrity of the bearing walls against
input acceleration.

23

A. FURUKAWA, J. KIYONO, K. TOKI

L

T(a) Front view

Element 3(b) Rear viewFig. 12 Seismic behavior of Flat20N, TLVTL(a) Front viewTL(a) Front view(b) Rear viewFig. 13 Seismic behavior of Flat20N, LTV

24

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKEc) Effect of the input ground motion direction
The building collapsed in the case of LTV but not in
the case of TLV. In the LTV case, the initial dominantcomponent, L, deformed the load-bearing wall in thein-plane direction and generated diagonal and horizontal cracks. The latter dominant component, T, thendeformed the load-bearing wall so that it collapsed in
the out-of-plane direction. The direction of the input
ground motion is important as it controls the occurrence of structural collapse.

cracks with peak ground acceleration in both cases.
The latter dominant component, T, did not generatecracks.
d) Effect of mortar strength
From a comparison of Figs. 10–15, the effect of mortar is clear. Increasing mortar strength is effective inpreventing the collapse of structures.
(4) Flat40N
The seismic behavior of the flat roof model with thickbearing walls and normal mortar (Flat40N) is shownin Figs. 16 and 17 for the TLV and LTV cases, respectively. There was tensile failure but no shear or compression failure.

(3) Flat20S
The seismic behavior of the flat roof model withstrong mortar (Flat20S) is shown in Figs. 14 and 15
for the TLV and LTV cases, respectively. There wastensile failure but no shear or compression failure.

a) TLV (Fig. 16)
Compared with the Flat20N model (Fig. 12), there
were few cracks in the bearing wall. However, thecrack patterns in the walls were similar to those of
Flat20N.
b) LTV (Fig. 17)
Compared with the Flat20N model (Fig. 13), thecrack patterns in the walls were similar, but the number of cracks and the out-of-plane deformation of the
bearing wall were smaller and the structure did notcollapse.
c) Effect of the thickness of the bearing wall
By increasing the thickness of the bearing wall, thenumber of cracks decreased and the structure did
not collapse in the case of LTV. Thus, increasing the
thickness of bearing walls is also effective in increasing the seismic resistance of a structure.

a) TLV (Fig. 14)
In the case of TLV, there was only tensile failure between the floor and the bearing wall. There was no
failure of the walls or roof. Failure occurred at 3.40 s,
at which time the L component was dominant. The L
component made the bearing wall move in the out-ofplane direction and generated cracks at the base of thewall.
b) LTV (Fig. 15)
In the case of LTV, diagonal and horizontal cracks
formed in the bearing wall at 3.40 s, at which time
the L component was dominant. The L component deformed the bearing wall in the in-plane direction andgenerated these cracks.
c) Effect of the input ground motion direction
For the flat roof models with strong mortar, the initial
dominant component, L, generated a small number ofL

T
Fig. 14 Seismic behavior of Flat20S, TLV

T

L
Fig. 15 Seismic behavior of Flat20S, LTV
25

A. FURUKAWA, J. KIYONO, K. TOKI

L

T
Fig. 16 Seismic behavior of Flat40N, TLV

L

T
Fig. 17 Seismic behavior of Flat40N, LTV(5) Flat20R
The seismic behavior of the flat roof model withwooden columns and beams and normal mortar(Flat20R) is shown in Figs. 18 and 19 for the TLV and
LTV cases, respectively. There was tensile failure but
no shear or compression failure.

ber of cracks, especially diagonal cracks, decreased,
and the structure did not collapse in the case of LTV.
Reinforcement using wooden columns and beams isthus also found to be effective.
(6) Comparison of reinforcement measures
From a comparison of Flat20W, Flat20N, and Flat20S,
it is successfully simulated that increasing the mortar
strength effectively increases the seismic strength ofthe model, which is natural since failure is only introduced in the mortar. We can also see how the failure
behaviors differ depending on the mortar strength.

a) TLV (Fig. 18)
Compared with the Flat20N model (Fig. 12), there
were fewer cracks in the bearing wall. The crack patterns in the walls also differed and there were fewer
diagonal cracks.
b) LTV (Fig. 19)
Compared with the Flat20N model (Fig. 13), thecrack patterns in the walls differed. In addition, there
were fewer cracks and less out-of-plane deformationof the bearing wall, and the bearing wall did not collapse.
c) Effect of wooden columns and beams
By introducing wooden columns and beams, the num-L

In addition, from a comparison of Flat20N, Flat40N,
and Flat20R, it is found that increasing the thicknessof the bearing walls and reinforcement with wooden
columns and beams also increases the seismic
strength of the model even though many cracks are
generated in the walls.

T
Fig. 18 Seismic behavior of Flat20R, TLV

T

L
Fig. 19 Seismic behavior of Flat20R, LTV
26

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKEThe displacement histories of Flat20S, Flat40N, and
Flat20R are compared with those of Flat20N to investigate the effectiveness of increasing the mortar
strength, increasing the thickness of the walls, andintroducing wooden columns and beams.

thickness of the nonbearing wall is the same. On theother hand, the displacements of elements 2 and 3 on
the bearing walls differ between Flat20N and Flat40N
since the bearing walls of Flat40N are thicker. By increasing the thickness of the bearing wall, the out-ofplane deformation decreases, which prevents collapse
in the case of LTV.

Figures 20 and 21 show the histories of displacementof elements 1, 2, and 3. Element 1 is at the center ofthe nonbearing wall, so its out-of-plane direction, x,
is investigated. On the other hand, elements 2 and 3
are at the centers of the bearing walls, so their out-ofplane direction, y, is investigated. The height of theseelements is 2.0 m. Figures 20 and 21 show the results
for the TLV and LTV cases, respectively. Since theFlat20N model collapsed in the case of LTV, the history of Flat20N is off the scale in Figs. 21 (b) and (c).

Compared with the displacements of Flat20N and
Flat40N, those of Flat20R and Flat20S are muchsmaller for the three elements in both directions. The
effectiveness of increasing the mortar strength and reinforcement with wooden columns beams is clear andis superior to increasing the thickness of the bearing
walls. Flat20S has less displacement, so increasing
the mortar strength is found to be the most effective
approach. Flat20R is the second-best model in termsof reducing the displacement response.

By comparing Fig. 20 (a) with Fig. 21 (a), it can beseen that the LTV case (Fig. 21 (a) ) has larger amplitudes in the x direction at 0–3.5 s because the larger L
component results in larger displacement. Similarly,
the TLV case (Figs. 20 (b) and (c) ) has larger amplitudes in the y direction at 0–3.5 s. Later, there is atendency for the amplitude in the other direction to belarger.

(7) Summary of the failure propagation process
The failure propagation process of the flat roof models is summarized as follows.
1) Diagonal and horizontal cracks form in walls fromthe corners of openings.
2) These cracks increase the out-of-plane deformation of the roof.
3) The out-of-plane deformation induces verticalcracks in the corners.

Comparing Flat20N and Flat40N, the displacement ofelement 1 on the nonbearing wall is similar since the(a) Element 1, x(b) Element 2, y(c) Element 3, yFig. 20 History of displacement, TLV(a) Element 1, x(b) Element 2, yFig. 21 History of displacement, LTV
27(c) Element 3, y

A. FURUKAWA, J. KIYONO, K. TOKI

4) The vertical cracks separate wall-to-wall connections.
5) The out-of-plane deformation of walls increases
since walls can move more independently.
6) The increase in out-of-plane deformation increases
the P–∆ effect of the roof weight.
7) The wall falls with the roof.

direction owing to the input ground motion. No clear
difference due to the direction of the input groundmotion is seen since the structure collapsed mainly
under its own weight owing to the instability of theroof.
(2) Vault20N
The seismic behavior of the vault roof model with normal mortar (Vault20N) is shown in Figs. 24 and 25 for
the TLV and LTV cases, respectively. There was only
tensile failure and no shear or compression failure.

4.2 Failure propagation of vault roof models(1) Vault20W
The seismic behavior of the vault roof model withweak mortar (Vault20W) is shown in Figs. 22 and 23
for the TLV and LTV cases, respectively. The process of failure propagation is expressed in the samemanner as for the flat roof models. There was neither
shear nor compression failure.

It is confirmed that Vault20N was not damaged if
there was no ground motion. However, small initial
ground motion triggered movement of the roof, andthis unbalanced the weight distribution and the roof
collapsed. Since the vault roof pushed the bearing
walls outward, the bearing wall with openings easily
cracked with small input ground motion at 1.20 s as
seen in the picture. The difference in the failure process due to the direction of the input ground motion
was small. In both cases, the bearing walls started
collapsing with the small initial ground motion, andthe nonbearing walls collapsed later with larger amplitudes of ground motion.

In both cases of input ground motion direction, thevault roof collapsed under its own weight. The vault
roof pushed the two supporting bearing walls outward, so the walls collapsed in the out-of-plane direction as soon as the simulation started. The nonbearing
walls did not collapse under their own weight because
they did not support the weight of the roof. Instead,the nonbearing walls collapsed in the out-of-plane

L

T

Fig. 22 Seismic behavior of Vault20W, TLV

T

L

Fig. 23 Seismic behavior of Vault20W, LTV
28

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

L

T

Fig. 24 Seismic behavior of Vault20N, TLV

T

L

Fig. 25 Seismic behavior of Vault20N, LTV(3) Vault20S
The seismic behavior of the vault roof model withstrong mortar (Vault20S) is shown in Figs. 26 and 27
for the TLV and LTV cases, respectively. There wasonly tensile failure and no shear or compression failure.

(4) Vault40N
The seismic behavior of the vault roof model withthick bearing walls and normal mortar (Vault40N) is
shown in Figs. 28 and 29 for the TLV and LTV cases,
respectively. There was only tensile failure and nocompression failure.

Vault20S has strong mortar strength, so the structure did not collapse even though the bearing wallwith openings was more cracked than in the case of
Flat20S.

L

Since the bearing walls are thicker, the roof collapsed
more slowly than Vault20N as shown in Figs. 24 and25. In the case of TLV (Fig. 28), the two bearing walls
did not collapse. However, in the case of LTV (Fig.
29), the bearing wall with openings fractured sinceT
Fig. 26 Seismic behavior of Vault20S, TLV

T

L
Fig. 27 Seismic behavior of Vault20S, LTV
29

A. FURUKAWA, J. KIYONO, K. TOKI

L

T

Fig. 28 Seismic behavior of Vault40N, TLV

T

L

Fig. 29 Seismic behavior of Vault40N, LTV
the dominant T component acted in the out-of-plane
direction of the bearing wall.

earthquake motion and collapsed during the period oflarge-amplitude ground motion. The wall was reinforced with wooden columns, but the bricks aroundthe openings scattered.

For the flat roof models, structural collapse depends
on the integrity of the bearing walls. Therefore, increasing the thickness of the bearing walls prevented
the collapse of the roof in the case of LTV. However,for the vault roof models, the thicker bearing walls
slow the collapse of the roof and they are themselves
less damaged, but this does not always prevent thecollapse of the roof. More specifically, if the roof
itself is weak or unstable, it collapses irrespective of
the strength of the bearing walls. Reinforcement of
the roof is thus necessary in the first place.

(6) Comparison of reinforcement measures
The displacement histories of Vault20S, Vault40N,
and Vault20R are compared with those of Vault20N to
investigate the effectiveness of increasing the mortar
strength, increasing the thickness of the bearing walls,
and introducing wooden columns and beams. Figures
32 and 33 show the histories of displacement of elements 1, 2, and 3.
Comparing Vault20N and Vault40N, it can be seenthat elements 1 and 2 fell in both cases but the elements of Vault40N fell later than those of Vault20N,
and that element 3 of Vault40N did not fall while thatof Vault20N did. These findings confirm the effectiveness of increasing the thickness of bearing walls.

(5) Vault20R
The seismic behavior of the reinforced vault roof
model with normal mortar (Vault20R) is shown in
Figs. 30 and 31 for the TLV and LTV cases, respectively. There was only tensile failure and no shear orcompression failure.

Comparing Vault40N and Vault20R, it can be seen thatelement 1 of Vault40N fell but that of Vault20R did
not since the nonbearing walls of Vault20R were reinforced and those of Vault40N were not. As for element
3, which did not fall in either model, the displacement

Since the four roofs are supported by wooden columns and surrounded by wooden beams, the roof collapsed more slowly than in the case of Vault40N (Figs.
28 and 29). However, the roof did not survive the30

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

L

T

Fig. 30 Seismic behavior of Vault20R, TLV

T

L

Fig. 31 Seismic behavior of Vault20R, LTV(a) Element 1, x(b) Element 2, y(c) Element 3, yFig. 32 History of displacement, TLV(a) Element 1, x(b) Element 2, y(c) Element 3, yFig. 33 History of displacement, LTV
of Vault20R was slightly larger than that of Vault40Nin the case of LTV since the great weight of the roof
pushed the bearing wall in the out-of-plane direction
and the thinner walls with wooden columns deformed
more easily than the thick walls without columns.

among all models. It can be seen that increasing the
mortar strength is the most effective approach.
(7) Summary of the failure propagation process
The failure propagation process of the vault roof models is different from that of the flat roof models, and issummarized as follows.

The displacement of Vault20S was the smallest
31

A. FURUKAWA, J. KIYONO, K. TOKI

1) The weight of the roof pushes the bearing walls in
the out-of-plane direction.
2) With small deformation of the buildings inducedby small initial ground motion, the heavy weightof the vaulted roof increases the out-of-plane deformation of the walls owing to the P–∆ effect.
3) The instability of the roof increases, and the roof
collapses.
4) If the bearing walls do not have sufficient strength,
they collapse due to the collapse of the vaulted
roof.

(2) Simulation results
For the flat roof models, the above collapse mechanism is confirmed by the simulation. For example,
Flat20N in the case of LTV collapsed because thebearing wall fractured, but Flat20N in the case of
TLV did not collapse because the bearing wall did
not. Moreover, by increasing the thickness of the
bearing wall of Flat20N, structural failure of Flat40N
was avoided even though the nonbearing wall was thesame.
For the vault roof models, thick bearing walls andwooden columns and beams had little effect. They
only delayed the failure, and collapse of the roof
could not be avoided. This is because the roof of theanalyzed model itself was weak. The roof requires
minimum levels of strength and integrity before reinforcements of bearing walls are effective.

4.3 Failure mechanism of past earthquake damage
and simulation results(1) Past earthquake damage
Photographs 1 (a) and (b) show two houses damagedby earthquakes. Photograph 1 (a) is of a house damaged during the 2003 Bam earthquake in Iran. Photograph 1 (b) is of a house damaged during the 1976
Friuli earthquake in Italy.
In Photo 1 (a), the nonbearing walls have been demolished, but the roof has not collapsed because thebearing walls are supported by the adjacent buildings.
In Photo 1 (b), the nonbearing walls are still standing, even though the roof and bearing walls havecollapsed. These two cases indicate that the stabilityof the roof depends on the two bearing walls and that
failure of the roof can be avoided if the bearing walls
are supported somehow, regardless of damage to nonbearing walls.

It is noted that the comparison between the analyticalresults for flat and vaulted roof models does not indicate that the vaulted roof model is far more vulnerable to earthquakes. The analyzed flat roof model has
much lighter roofs made only by wooden beams, but
in reality, mud and trees branches are used to cover
flat roofs depending on the area, and this makes theroof heavier and more vulnerable. Moreover, the vulnerability of the analyzed vault roof models changes
depending on the roof thickness, so we cannot simply
discuss the superiority of the flat and vaulted roof
models through these simulation results.

(a) Bam Earthquake, Iran (2003)(b) Friuli Earthquake, Italy (1976)

Photo 1 Buildings damaged by earthquakes
32

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

4.4 Scattering of bricks after earthquakes(1) Flat roof models
The areas of scattered fallen bricks are compared.
Overhead views from the basement of flat roof models are shown in Figs. 34 and 35.

to vibrate more greatly, and this may increase the areaof scattered fallen bricks. The effect of the directionof the input ground motion is clear from a comparisonof Figs. 34 and 35, since the flat roof models are damaged mainly by the earthquake ground motion.

In the case of LTV, Flat20W was severely damaged,
and the bricks in the nonbearing walls were mainly
scattered around these walls (Fig. 34 (a) ). In the caseof TLV, Flat20W and Flat20N collapsed, and the areaof scattered bricks for Flat20W was larger than thatfor Flat20N owing to the weakness of the mortar. It
can be seen that the mortar strength affects not only
the failure propagation process but also the area of
scattering. Walls connected with weaker mortar tend

(2) Vault roof models
Overhead views from the basement of the vault roof
models are shown in Figs. 36 and 37. No effect ofthe direction of the input ground motion is seen in
the comparison of Figs. 36 and 37 because the vault
roof models collapsed owing to the weight imbalance
induced by small deformation due to the initial earthquake ground motion.

TL(a) Flat20W(b) Flat20N(c) Flat20S(d) Flat40N(e) Flat20R

Fig. 34 View from the basement of flat roof models, TLVLT(a) Flat20W(b) Flat20N(c) Flat20S(d) Flat40N

Fig. 35 View from the basement of flat roof models, LTV

33(e) Flat20R

g

A. FURUKAWA, J. KIYONO, K. TOKITL(a) Vault20W(b) Vault20N(c) Vault20S(d) Vault40N(e) Vault20R

Fig. 36 View from the basement of vault roof models, TLVLT(a) Vault20W(b) Vault20N(c) Vault20S(d) Vault40N(e) Vault20R

Fig. 37 View from the basement of vault roof models, LTVof the bearing walls is small, collapse of the roof can
be avoided. Therefore, a reinforcement that reduces
out-of-plane deformation is effective. Increasing the
thickness of the bearing walls and introducing wooden columns and beams are only partial reinforcements, but both are effective since they reduce outof-plane deformation. Increasing the mortar strength
is the most effective approach since it increases thewhole structural strength and reduces deformation on
the whole. It is also found that the failure propagationof the flat roof models is dependent on the directionsof the input ground motion since the out-of-plane deformation of the bearing walls depends on the characteristics of the ground motion.

5. Conclusions
In this paper, failure propagation of masonry buildings during earthquakes is simulated using a refined
DEM that can simulate three-dimensional elastic,
failure, and collapse behaviors of structures. The seismic behaviors of models with a flat roof or vault roof
are simulated and their failure mechanisms are discussed. In addition, the influence of the direction of
the input ground motion on the failure propagation isinvestigated. Moreover, the effectiveness of three reinforcement measures was compared. One measure isincreasing the mortar strength, the second is increasing the thickness of the bearing walls, and the third isintroducing wooden columns and beams, which havehigher tensile strength than the bricks.

For the vault roof models, the collapse of the roof is
induced by its heavy weight and instability triggeredby movement due to ground motion. Increasing the
mortar strength is the most effective approach since
it directly increases the strength of the roof. Increasing the thickness of the bearing walls and introducing
wooden columns and beams are effective in reducing
the instability of the roof since they reduce the out-

Even though the simulation is limited to one set ofearthquake records, the following findings are obtained through the simulations.
The stability of the roof of the flat roof models is dependent on the out-of-plane deformation of the bearing walls. Therefore, if the out-of-plane deformation
34

NUMERICAL SIMULATION OF THE FAILURE PROPAGATION OF MASONRY BUILDINGS DURING AN EARTHQUAKE

of-plane deformation of the bearing walls supporting
the roof. However, since these measures do not increase the strength of the roof, the roof collapses eventhough the bearing walls remain.

Luxembourg, 1998.
Furukawa A, Kiyono, J, Toki, K. Proposal of a
Numerical Simulation Method for Elastic, Failure and Collapse Behaviors of Structures andits Application to Seismic Response Analysis of
Masonry Walls. Journal of Disaster Research,
Vol. 6, No. 1, 2011 (to be published).
[5] Zienkiewicz, OC, Taylor, RL. The Finite Element Method. 5th ed. Oxford, UK: Butterworth
Heinemann vol. 1, 2, 3, 2000.
[6] Ghannad MA, Bakhshi A, Mousavi Eshkiki SE,
Khosravifar A, Bozorgnia Y, Taheri Behbahani
AA. A study on seismic vulnerability of rural
houses in Iran. Proceedings of the First European Conference on Earthquake Engineering andSeismology, paper no. 680, 2006.
[7] Cao Z, Watanabe H. Earthquake response prediction and retrofitting techniques of adobe
structures. Proceedings of the 13th World Conference on Earthquake Engineering, paper no.
2594, 2004.
[8] Gambarotta, L., Lagomarsino, S. Damage Models for the Seismic Response of Brick Masonry
Shear Walls. Part I: The Mortar Joint Model andits Applications. Earthquake Engineering andStructural Dynamics, Vol. 26, 423-429, 1997.
[9] Lourenco PB, Pots JG. Multisurface Interface
Model for Analysis of Masonry Structures. Journal of Engineering Mechanics, 123:660-668,
1997.
[10] Cundall PA, Strack ODL. A Discrete Numerical
Model for Granular Assemblies. Geotechnique,29:47-65, 1979.
[11] Shi GH. Block System Modeling by Discontinuous Deformation Analysis. Southampton, UK:
Computational Mechanics Publications, 1993.
[12] Kawai T. New Element Models in Discrete
Structural Analysis. Journal of the Society ofNaval Architects of Japan, 141:187-193, 1977.
[13] Alexandris A, Protopapa E, Psycharis I. Collapse mechanisms of masonry buildings derived
by the distinct element method. Proceedings ofthe 13th World Conference on Earthquake Engineering, paper no. 548, 2004.
[14] Bakhshi A, Bozorgnia Y, Ghannad MA, Khosravifar A, Eshkiki EM, Rofooei FR, Behbahani
AT. Seismic vulnerability of traditional housesin Iran. Seismo Adobe 2005 Seminar. Available[4]

It is found that the roof needs minimum levels ofstrength and integrity before the reinforcement ofbearing walls can become effective. If the roof has theminimum strength and integrity, the stability of theroof depends on the two bearing walls and failure of
the roof is avoided if the bearing walls are supported
somehow regardless of damage to nonbearing walls.
Increasing the thickness of the bearing walls and introducing wooden beams and columns are both effective in reducing the out-of-plane deformation of the
bearing walls. Increasing the mortar strength is themost effective approach since it increases the whole
structural integrity including the walls and the roof.
In a future study, we would like to develop a code for
evaluating the failure of the brick itself as well as an
appropriate damping model.
ACKNOWLEDGEMENTS
This work was supported by a Grant-in-Aid for the
Global COE Program “Global COE for Education,
Research and Development of Strategy for Disaster
Mitigation of Cultural Heritage and Historic Cities(Ritsumeikan University) ” from the Ministry of Education, Culture, Sports, Science and Technology ofJapan.
REFERENCES
[1]

[2]
[3]

Giardini D, Gruenthai G, Shedlock K, Zhang
P. The GSHAP Global Seismic Hazard Map.
In: Lee WH-K, Kanamori H, Jennings PC,
Kisslinger C, editors. International Handbook ofEarthquake and Engineering Seismology. London UK: Academic Press, 1233-1239 vol. 81B,
2003.
Coburn A, Spence R. Earthquake Protection.
2nd ed. Chichester: John Wiley and Sons, 2002.
Grunthal, G. (ed.) European Macroseismic Scale
1998 (EMS-98). Cahiers du Centre Européen
de Géodynamique et de Séismologie 15. Centre
Européen de Géodynamique et de Séismologie,
35

A. FURUKAWA, J. KIYONO, K. TOKI

at: http://www.3ei.org/images/peru05paper.pdf.
Accessed January 28, 2010.
[15] Lourenco PB. Analysis of Masonry Structureswith Interface Elements, Theory and Applications. Delft University of Technology, Facultyof Civil Engineering, TU-DELFT report no. 0321-22-0-01, 1994.
[16] Courant RL. Variational Methods for the Solution of Problems of Equilibrium and Vibration.
Bulletin of the American Mathematical Society,49:1-23, 1943.
[17] Kiyono J, Kalantari A. Collapse Mechanism of

Adobe and Masonry Structures during the 2003
Iran Bam Earthquake. Bulletin of Earthquake
Research Institute, University of Tokyo, 79:157161, 2004.
[18] Marzahn G. The Shear Strength of Dry-Stacked
Masonry Walls. Leipzig Annual Civil Engineering Report (LACER) No. 3. Germany: University of Leipzig, 1998.
[19] BHRC. The Very Urgent Preliminary Report on
Bam Earthquake of Dec 26, 2003. http://www.
bhrc.gov.ir, 2003.

36