Discrete Rigid Block Analysis to Assess Settlement Induced Damage in Unreinforced Masonry Façades

Abstract

In this study, a system of discontinuous rigid blocks is employed to simulate the possible damage mechanisms in unreinforced masonry (URM) façades and load-bearing frame systems subjected to settlement using the discrete element method (DEM). First, the employed modeling strategy is validated utilizing the available experimental results presented in the literature. Once there is a good agreement between the computational models and experimental findings, a sensitivity analysis is performed to quantify the influence of the input parameters defined in the DEM-based numerical model. Finally, the proposed modeling strategy is further utilized to assess the damage pattern that may develop in a URM façade due to uniform and non-uniform settlement profiles. The results of this study clearly show that the discrete rigid block analysis (D-RBA) provides robust numerical solutions that can be employed to visualize and assess the possible damage patterns and related collapse mechanisms of URM masonry systems as an alternative modeling strategy to standard continuum-based solutions.

1. Introduction

Unreinforced masonry (URM) construction constitutes the vast majority of the world’s building stock, including historical structures worldwide. Although mechanical and thermal characteristics of masonry provide desirable construction features, it has a relatively weak tensile strength at mortar joints. This salient feature of masonry may yield local and global structural problems in masonry buildings, such as cracking, joint dislocation, and sliding failure in the load-bearing component of masonry structures (e.g., walls, arches, domes, etc.) throughout their service life. Additionally, extreme natural events (e.g., earthquakes and tornados), high-temperature fluctuations, and severe changes in boundary conditions (e.g., soil settlements and landslides) may trigger tension failure at the mortar joints and result in a partial or total collapse of masonry structures.

In this paper, settlement-induced damage progression (SIDP) in URM facades is studied using the discrete element method (DEM). This study demonstrates the possible crack patterns and failure modes of unreinforced masonry structures subjected to various forms of settlement.

The accuracy of the DEM models is ensured by comparing them to prior experimental results, resulting in an ability to test a large variety of scenarios. This, in turn, allows studying, in detail, the contribution of multiple variables (e.g., joint tensile strength and joint frictional resistance). As a result, the paper provides a better understanding of the cracking behavior and damage progression in masonry structures due to soil settlement.

Since the early 1980s, both experimental and computational works have been presented in the literature regarding settlement-induced failures in masonry structures. Boscardin et al. [1] used analytic models and field data to develop a procedure that helps evaluate the tolerance of a brick-bearing wall for settlement caused by tunneling. A homogeneous isotropic material behavior was assumed, and the brick-bearing wall was modeled as a linear-elastic deep beam. The ground deflections were imposed on the structure, and the critical tensile strain concept was used to evaluate the damage. The results showed that the tolerance to differential settlement decreases when a structure is subjected to increasing lateral strain, indicating buildings sited next to excavation are less tolerant than similar structures settling under their self-weight because of the lateral strain caused by responses from excavation. Burland and Wroth [2] discussed the damage in load-bearing brickwork walls due to settlement. Relative deflections were estimated based on the concept of critical tensile strain. Son et al. [3] performed damage assessment of brick-bearing, brick-infilled, and open-frame structures subjected to excavation-induced ground settlements. A two-dimensional (2D) discrete element method was utilized to simulate cracking in masonry units, and the emphasis was given to four-story structures on a shallow foundation. It was aimed to evaluate the elastic and non-linear response of the different structural types with different soil conditions (i.e., soft and stiff) considering progressive ground settlement. The results indicated that deriving conclusions by only looking at the elastic behavior of the structure can be misleading; thus, crack-induced damage progression and the soil-structure interaction phenomena should also be taken into account to obtain a better understanding of the relationship between structure and soil settlement. Recently, Camós et al. [4] analyzed the structural response of one-story masonry buildings affected by ground movements from L9 Metro tunnel construction in Barcelona. Two computational modeling strategies were utilized: the equivalent beam method and nonlinear finite element analysis (FEA). It was found that the largest plastic strain development occurs near the corners of doors and windows, where a high-stress concentration is expected. Furthermore, real tunnel-induced cracks obtained from the in situ damage survey were closely captured by evaluating the maximum crack openings determined by nonlinear FEA, demonstrated as diagonal and vertical cracks in windows and doors. Giardina et al. [5] performed a laboratory test carried out on a scaled masonry façade subjected to tunneling-induced settlements, where particular attention was given to the soil-structure interaction phenomena. The experiment reproduced a circular tunnel excavation in soft soil under a masonry structure. A pre-defined settlement profile was applied to a scaled building model by continuously monitoring the deflection and crack progression in the structure. The experiment showed that soil-structure interaction plays an important role in overall structural response because the stress redistribution within the façade depends on the relative stiffness between the façade and soil interface. In another recent research, interpreted crack patterns of masonry buildings subjected to landslide-induced settlements is analyzed using the load path method (LPM) by Palmisano [6]. The LPM searches for the load path that contains the lowest value of total strain energy instead of the exact path. Using macro-models obtained by the LPM, the structural behavior of the masonry walls subjected to landslides was evaluated.

In contrast to the existing literature summarized above, this study examines the cracking phenomenon in masonry walls via a discontinuum-based modeling approach where masonry units are represented as discrete (or distinct) rigid blocks, and mortar joints are replicated as zero-thickness interfaces. The described DEM-based analysis of masonry constructions has been used to assess various structures as discussed in a seminal article by Lemos [7].

First, the proposed modeling strategy is validated based on the small and large-scale experimental results, and then it is applied to assess a typical existing URM facade in downtown Omaha, Nebraska (United States). Once a validated model is established, a series of sensitivity analyses on the joint tensile strength and friction angle is performed. A brief background of DEM, validation of the computational approach, parametric analyses, and an application on real URM façade are presented in the following sections, respectively. Finally, the results are discussed, and necessary inferences are made regarding the damage progression and input material properties influencing the displacement response of masonry walls subjected to various forms of support settlement.

2. Computational Modeling Approach: Discrete Rigid Block Analysis (D-RBA)

Masonry is a highly non-linear, orthotropic, and heterogeneous material that makes it difficult to analyze accurately via simple methods. Moreover, frequently, the cracks are localized at the mortar joints, following a distinct crack propagation path. As a result of this, morphological features of masonry construction have an important role in the structural behavior and capacity of URM walls and other types of structural systems [8,9,10]. Although continuum-based solutions (e.g., standard FEA) provide valuable information about the capacity and smeared damage representation, more sophisticated algorithms or interface elements embedded into the FE model are required to obtain localized failures in FEA, as discussed in other studies [11,12,13]. Here, we propose an alternative approach, called discrete rigid block analysis (D-RBA), to simulate masonry façades based on the DEM. It should be noted that since the lack of tensile strength governs the examined collapse mechanism, no compression failure is implemented at the joints or within the blocks, as proven to be a valuable approach in previous studies [14,15,16].

The employed computational method is first developed by Cundall [17] to simulate jointed rock masses. In the last several decades, DEM has been used to analyze various masonry structures with different levels of complexity, considering deformable and rigid blocks in the material and structural level [18,19,20,21,22]. The main motivation of DEM is to simulate the interaction between regular or irregular-shaped blocks that go under large deformations using automatic contact detection algorithms to capture sliding, joint opening, and total contact loss. The numerical procedure of D-RBA relies on integrating translational and rotational equations of motions of rigid block center of mass. Note that a 3D rigid block has six degrees of freedom, three translational and three rotational. Simply, new velocities, evaluated at the end-points of the time step, are computed explicitly using the central difference method as given in Equations (1) and (2). The time step (∆t) is divided into mid-intervals and denoted as $t^{+}=t+\Delta t/2$ and $t^{-}=t-\Delta t/2$. The quasi-static solutions are obtained by applying artificial damping to prevent any unnecessary noise in the solutions.

$${\dot{u}}_i^{t+}=\left[\left(1-\frac{c\Delta t}{2}\right){\dot{u}}_i^{t-}+\left(\frac{\Sigma F_i^t}{m}+g_i\right)\Delta t\right]\frac{1}{1+\frac{c\Delta t}{2}}$$

(1)

$${\omega }_i^{t+}=\left[\left(1-\frac{c\Delta t}{2}\right){\omega }_i^{t-}+\frac{\Sigma M_i^t}{I}\Delta t\right]\frac{1}{1+\frac{c\Delta t}{2}}$$

(2)

where $u$, $F$, ω, and $M$ are the displacement, force, angular velocity, and moment vectors. Moreover, $g$, $I$, and $c$ represent gravitational acceleration vector, the approximate moment of inertia, and the mass-propositional damping constant.

Once the new velocities are obtained, the block positions are updated based on the calculated translational and rotation increments ($\Delta u_i={\dot{u}}_i^{t+}\Delta t,\Delta {\theta }_i={\omega }_i^{t+}\Delta t$). Accordingly, relative point-contact displacements are determined among the adjacent blocks, where each point-contact orthogonal springs are defined, illustrated in Figure 1. In the normal direction, tensile strength is limited via pre-defined value ($f_T$), whereas linearly elastic behavior is assumed in compression, as mentioned earlier. Coulomb slip-joint model is used in the shear direction, which requires contact friction angle and cohesion, denoted as $c$ and $\varphi$, respectively (Figure 1), to calculate the shear strength depending on the vertical stress ($\sigma )$ at the joint (i.e., $\tau =c+\sigma tan\varphi$). It is worth noting that the assigned normal ($k_n$) and shear ($k_s$) spring stiffnesses control elastic response at the joints, as well as the amount of overlapping between blocks (in compression). Note that the contact stress increments are computed utilizing the relative displacements between blocks (i.e., $\Delta \sigma =k_n\Delta {\delta }_n$,$\Delta \tau =k_s\Delta {\delta }_n$, see Figure 1) [23]. Then, the new contact stresses are obtained (${\sigma }^{new}={\sigma }^{old}+\Delta \sigma ,{\tau }^{new}={\tau }^{old}+\Delta \tau$) and updated according to the considered contact constitutive model.

Finally, the contact stresses are multiplied with the assigned contact area to compute the forces employed in the equations of motions. Since the executed explicit computational procedure is conditionally stable, $\Delta t$ should be equal to or smaller than the critical time step ($\Delta t_{cr}$) to keep the numerical stability during calculations. Here, the time step is determined as $0.1\sqrt{m_{min}/k_{n,max}}$, considering a similar analogy to a simple degree of freedom system, suggested in [24]. Throughout this research, a commercial three-dimension discrete element code (3DEC-7.0), developed by ITASCA, is used [25]. In the next section, validation of the applied discontinuum-based modeling approach is presented.

3. Validation Studies

Two experimental studies, presented by Portioli and Cascini [26] and Giardina et al. [5], are utilized to validate our modeling strategy. A simple URM pier and spandrel system was tested in the former test, while in the latter, a scaled URM façade was built and exposed to non-uniform support settlement.

Validation Study—1: The small-scale dry-joint masonry portal frame, made up of tuff stones, was built and tested using vertically moving support by Portioli and Cascini [25]. It aimed to better understand the ultimate displacement and collapse mechanism of stone masonry portal frames exposed to support settlement. In the reference study, different boundary conditions (i.e., fixed and roller) were considered during the experimental campaign; however, only one of them (in which the moving support is fixed in the lateral direction) is employed in this study for validation purposes. The modeled portal frame comprises 12-courses in height and nine bricks in length with a three brick-length opening (brick dimension: 100 mm × 100 mm × 50 mm). Note that vertical displacement was imposed to the left masonry pier, where the supporting plate was restricted along the lateral direction, and supported reactions were monitored during the experiment. For further details regarding the test setup, the readers are referred to the reference study [26].

The applicable material properties (i.e., density and contact surface friction angle) are directly taken from the cited study, while the contact stiffness ($k_n$) is predicted as 1 GPa/m, and the normal to shear stiffness ratio is taken as 0.4. As shown in Figure 2, the applied discontinuum-based rigid block approach provides good predictions in terms of collapse mechanism (Figure 2a) and force-reaction, recorded at the base of the settled masonry pier (see Figure 2b-left). Additionally, a sensitivity analysis is performed by considering similar discontinuum-models in the literature (e.g., [27,28]) to better observe the influence of contact stiffness. The results show parallel base reactions together with no considerable structural behavior difference. A smoother base reaction vs. vertical displacement curve is obtained for a lower stiffness ($k_n=0.1$ GPa/m), whereas relatively sudden reaction change (in other words, more brittle macro behavior) is observed for a higher stiffness ($k_n=25$ GPa/m), given in Figure 2b-right.

Validation Study—2: Once the D-RBA is validated on a small-scale pier-spandrel system, a more complex structural configuration is analyzed (Figure 3a). A recent experimental work presented by Giardina et al. [5] is utilized to validate our computational approach where the tunneling-induced settlement effects on masonry façades were investigated. The experiment was set up to replicate the event, where a circular tunnel being excavated in soft soil under a masonry structure on a wooden pile foundation. The structure was located in the hogging zone of the tunnel-induced settlement profile (downward concavity). Note that a 1/10th scaled masonry façade was constructed with 1010 scaled solid bricks (dimensions of 25 mm × 40 mm × 50 mm) and 2 mm thick layers of low strength lime mortar. The base of the scaled masonry façade was connected to a flexible steel profile to be able to control the settlement beneath the structure. The interface stiffness was taken into account by simulating soil-structure interaction by inserting rubber between the masonry façade and the steel profile. Several monitoring points on the façade were monitored during the experiment to trace the lateral and vertical deformations, as illustrated in Figure 3.

First, dead loads (self-weight + floor loads) were applied (readers referred to reference study for further details). Then, the steel profile was adjusted in 23 steps of 0.5 mm, reaching a final displacement of 11.5 mm on the scaled model. The settlement profile is represented by the gaussian error function presented by Peck [29] given in Equation (3).

$$S_v\left(\mathit{x}\right)=S_{v,max}{e}^{-\frac{x^2}{2i_x^2}}$$

(3)

where $S_{v,max}$ is the maximum settlement measured above the tunnel axis, $\mathit{x}$ is the horizontal distance from the tunnel axis and $i_x$ denotes the horizontal distance between the vertical tunnel axis and the point of inflection.

Here, $i_x$ is taken as 8 m, determined by Giardina et al. [5]. The same loading procedure (i.e., dead load + support settlement) is followed through the discrete rigid block analysis. It is worth noting that in the discontinuum-based model, the continuous settlement profile is replaced with a series of discrete vertical displacements, which are defined for each support block as demonstrated in Figure 3b. Once the support settlement profile is determined, the vertical velocity is prescribed to the corresponding support blocks with a very low displacement rate to obtain quasi-static solutions from the motion equations, as explained in the previous section. The contact properties are determined in line with the reference study, where joint tensile strength, cohesion, and friction angle are defined as 0.1 × 10⁶ Pa, 0.15 × 10⁶ Pa, and 35 degrees, respectively, whereas the contact stiffness is assumed 11 GPa/m, determined via sensitivity analysis as performed in the previous validation study with the same normal to shear stiffness ratio.

The discrete rigid block models show similar trends compared to the experimental results given in Figure 4. The outcomes of this validation study show a strong similarity between the horizontal and vertical displacement histories, recorded from monitoring points ($M{P}_1$ and $M{P}_2$), see Figure 4a. Moreover, observed crack pattern and global structural mechanism obtained from D-RBA exhibit parallel behavior compared to experimental findings presented in the reference study [5].

4. Parametric Analysis

In this section, a sensitivity analysis is performed to better understand the effect of two essential input parameters: joint tensile strength ($f_t$) and the friction angle ($\varphi$), used in the D-RBA. It is important to note that the ratio of cohesion to tensile strength kept constant as 1.5 during the parametric analysis, in line with the other studies [30,31].

Results show that tensile and frictional resistance affect the displacement at different magnitudes, as depicted in Figure 5. The change in friction angle from 20 degrees to 50 degrees only causes a 6% decrease in horizontal deflection, while changing tension from 0.01 MPa to 1 MPa causes a 29% percent decrease. Overall, it is important to note that horizontal deflection connected to change in tension is around 6% less than horizontal deflection connected to a change in frictional angle. Vertical displacement only decreased by 1% when changing the friction angle from 20 degrees to 50 degrees and 7% when changing tension from 0.01 MPa to 1 MPa. Hence, the sensitivity analysis results suggest that the joint tensile strength has a much larger influence on the results because the forces acting upon the mortar joints trigger a rapid tensile failure and yield different force-distribution paths at the joints. This finding also means that tensile failure becomes dominant over shear failure (sliding), which further underlines the importance of contact parameters on the structural response. Therefore, reliable data regarding the joint strength and quality are necessary as inputs for the computational models to obtain accurate predictions of structural behavior.

5. Applications of the D-RBA: Crack Patterns in URM Façade

In this section, the validated modeling strategy is applied to observe the possible crack patterns that can develop in a URM Façade with low-bond joint tensile strength. First, a hypothetical clay brick masonry façade is generated, inspired by a typical 19–20th century historic residential building located in the downtown area of Omaha, Nebraska, United States (Figure 6). The geometrical features are approximately measured based on the structural survey, where inner and outer column widths are determined as 0.65 and 0.75 m, respectively. The openings are regular, and the structure is symmetric along the vertical axis (further details regarding the geometrical properties can be found in Figure 6). Moreover, the material and contact properties, given in

Table 1. Material and contact properties.

$\mathit{\rho } \left(\frac{\mathit{k}\mathit{g}}{{\mathit{m}}^3}\right)$	${\mathit{k}}_{\mathit{n}}\left(\frac{\mathit{P}\mathit{a}}{\mathit{m}}\right)$	${\mathit{k}}_{\mathit{s}}\left(\frac{\mathit{P}\mathit{a}}{\mathit{m}}\right)$	${\mathit{f}}_{\mathit{t}}\left(\mathit{P}\mathit{a}\right)$	$\mathit{c}\left(\mathit{P}\mathit{a}\right)$	$\mathit{\varphi }(°)$
2000	10 × 109	4 × 109	0.1 × 106	$1.5f_t$	35

, are determined based on the related previous studies [15,16,32,33]. Note that only the self-weight of the structure is taken into account since the motivation of this research is to assess the failure mechanism due to support settlement, which is not affected by applied floor loads [34].

In total, six settlement-induced damage scenarios (from S-1 to S-6) are considered to better understand the corresponding damage pattern in the examined URM masonry façade. Specifically, the first three settlement profiles are defined as uniform vertical displacement, whereas, in the remaining scenarios non-uniform displacement field, along the length of the structure is considered. The maximum vertical displacement ($\Delta u_{max}$) is determined as 0.05 m, which gives the ultimate deflection ratio ($\Delta u_{max}/L_0$) of 0.00625, passing the very severe damage threshold according to the accepted damage criteria in the literature, which is discussed in detail by Serhal et al. [35]. The displacements are gradually applied by 0.005-m increments, and after each increment, it is ensured that the structure can stand under imposed settlement profile. As mentioned in the previous section, very low displacement rates (e.g., 1×10⁻³ m/s) are assigned to support blocks to avoid unnecessary structural vibrations and oscillation in the numerical solutions. Furthermore, from the computational cost point of view, each settlement analysis requires 2 h, using a computer system Intel^® Core^TM i7 CPU @ 2.7 GHz processor and 16 GB memory RAM.

The results of the D-RBA reveal that the uniform settlement profiles yield distinct localized diagonal cracks developing at the spandrels, starting from the early stages of the deflection. When the ultimate displacement is reached, the lateral deformation becomes evident, and an overturning behavior is noticed at the outer column, which can be seen in Figure 7a–c. Hence, even if we have a uniform vertical settlement, in case of slender columns (as in the analyzed URM façade), an overturning mechanism may be observed that can turn into a kinematic mechanism. In the case of non-uniform settlement scenarios (Figure 8a–c), almost doubled horizontal deflection with a lateral drift ratio (ξ) of 0.5–0.6% ($\xi =max. lateral displacement/H_3$) is obtained compared to uniform settlement for a given maximum vertical deflection. The results demonstrate a relatively different crack pattern compared to former analyses (S-1:S-3, see Figure 7a–c), and a less critical joint opening profile is noticed when the non-uniform support settlement is defined for columns, including $C_1$, $C_2$, $C_3$, and $C_4$, as shown in Figure 8c (S-6). This suggests that the internal force distribution among the column and spandrel system is shared more evenly compared to other uniform and non-uniform settlement scenarios (S-1:S-5), causing fewer joint openings in the URM façade.

Additionally, the effect of the in-plane slenderness ratio of the spandrel wall ($h/L_1$, see Figure 9a) is investigated to better understand its role on the cracking pattern and corresponding drift ratios. The experimental findings and computational models highlight that in-plane behavior of masonry walls subjected to lateral loading is a complex phenomenon and can be influenced by several factors, namely geometrical features (slenderness), material properties, boundary conditions, vertical loads, and morphology (or workmanship), as discussed in given references among many others [36,37,38,39,40]. However, compared to lateral loading, the settlement-induced failure mechanisms, and their relation to the macro-behavior of the URM Façades are less studied in the literature.

In this section, three slenderness ratio limits (i.e., $h/L_1$ < 1.0; 1.0 < $h/L_1$ < 1.5; $h/L_1$ > 1.5) are considered, and the second settlement scenario (S-2) is analyzed. Note that in the previous analyses, the slenderness ratio of spandrel walls (shown in Figure 7 and Figure 8) is less than 1.0 ($h/L_1$ ≈ 0.75). Although similar cracking patterns are observed in spandrel walls with different slenderness ratios, the considerable difference is noted in terms of drift ratios, calculated based on the lateral displacement of the monitoring point (MP, demonstrated in Figure 9a). The results indicate higher drift ratios when the slenderness ratio is larger than 1.5, as shown in Figure 9b.

6. Conclusions

This research presents a validated three-dimensional discontinuum-based modeling approach, called discrete rigid block analysis (D-RBA), to predict the structural behavior of unreinforced masonry facades subjected to various scenarios of soil settlement. Note that the proposed workflow can also be adopted in two-dimensional DEM-based rigid block analysis.

The applied modeling strategy is validated based on two experimental findings, where dry and mortar joint URM masonry systems are analyzed. In both cases, good agreement is found between the computational models and experimental results in terms of structural behavior and deformation histories. Furthermore, according to the sensitivity analyses performed using the validated numerical model, the joint tensile strength is found to be the most influential input parameter, significantly influencing the deformation behavior. Finally, the modeling approach is applied to the geometric representation of a historical URM building façade. Uniform and non-uniform settlement profiles are imposed on the support blocks, and corresponding cracking mechanisms are investigated. In all settlement-induced damage patterns, severely localized spandrel cracks are noticed when the ultimate deflection ratio is reached, which may threaten structural integrity. Furthermore, the influence of spandrel walls with different in-plane slenderness ratios are explored, and the lowest drift ratios are found when $h/L_1$ less than unity for the analyzed URM façade.

Hence, it is shown that the proposed D-RBA can be a robust solution that can be utilized in the damage assessment of URM façades suffering from settlement-induced cracks. In the future study, different façade openings and stiffness ratios regarding the soil-structure interaction will be considered, which will contribute to the current knowledge base on settlement-induced damage failures for unreinforced load-bearing masonry systems.