Effect of Precompression and Material Uncertainty on the In-Plane Behavior of URM Pier–Spandrel Systems

Theoretical and experimental studies on loadbearing masonry walls have shown the significant influence of the axial load level (i.e., precompression) and wall aspect ratio on in-plane lateral resistance. Nonetheless, the impact of the precompression and spatial variability of the material properties needs to be further investigated at the scale of walls with openings. This study presents a stochastic analysis of unreinforced (URM) pier–spandrel systems subjected to both axial loads on piers and lateral loads, considering the spatial variation in material properties. A discontinuum-based computational model was utilized to assess the force–displacement behavior of a benchmark pier–spandrel structure under different vertical precompression levels on piers. A total of 750 simulations were carried out to propagate material uncertainties in lateral load analysis. The proposed modeling strategy, based on the discrete element method, explicitly represents joint openings, sliding, and crushing phenomena at the contact points defined between the adjacent discrete rigid blocks. According to the validated computational modeling strategy, meaningful inferences were made regarding the effect of the precompression level on the maximum displacement and ultimate lateral load-carrying capacity of the benchmark URM pier–spandrel system. The results showed that vertical pressure on piers had considerable influence on the displacement ductility of the system while yielding less variation in the displacement capacity. Furthermore, the appealing feature of the spatial probabilistic analysis is noted in the variation in the lateral load-carrying capacity of the structural system.

Keywords: discrete element method; masonry; material uncertainty; pier–spandrel system; probabilistic analysis; Monte Carlo simulation; spatial variability

The vast majority of residential buildings and architectural heritage were built using unreinforced masonry (URM), a composite material consisting of units (stone, clay brick, earth blocks, etc.) and often a binding agent such as mortar. Each constituent of masonry shows significant uncertainty in terms of its mechanical properties with a highly nonlinear material behavior [[1]]. Therefore, a probabilistic approach considering predefined variations in material properties varying spatially can effectively predict the behavior and capacity of URM structures. Recently, the probabilistic analysis of masonry structures has been studied by various researchers [[2], [4], [6], [8], [10]]. In most of the existing works on the probabilistic analysis of masonry structures, the mechanical properties of masonry are assumed to be uniformly distributed within the structure, while only very few papers consider the spatial variation in those properties [[3], [9], [11], [13]]. As evidenced by those studies, the spatial variation could change the damage pattern and variation in response parameters. Thus, more research is necessary to understand and quantify such effects.

A typical combination of masonry composite consists of strong units and weak mortar joints, where mortar joints delineate weak planes prone to failure. This phenomenon yields localized cracks in the form of joint openings and sliding along the head and bed joints in URM walls. Conversely, weak unit–strong joint action causes opposite behavior, in which the cracks predominantly pass through the masonry units [[15], [17]]. In each combination of masonry constituents, progressive damage can be noticed in URM walls subjected to lateral loading due to successive failures at the masonry units and mortar joints. As witnessed in the past and recent earthquakes, URM buildings are vulnerable to seismic loads (unless reinforced), particularly when no proper connection is implemented between the structural elements and orthogonal walls, among many other factors [[18], [20]]. The in- and out-of-plane failure of masonry walls are the two common collapse modes that can develop independently or interact with each other based on the construction quality, wall configuration, and openings. However, when façades are properly connected to orthogonal walls and floor systems are significantly stiffer than URM walls, in-plane mechanisms become more likely to develop under earthquake forces that are distributed to the in-plane walls as proportional to their own lateral stiffness.

Different computational modeling techniques are available to simulate the in-plane behavior of URM walls [[21]]. The first approach considers the composite nature of masonry as a continuum (denoted as macromodeling), where the entire URM structure is modeled via an equivalent homogeneous medium with averaged material properties corresponding to the main features of masonry (e.g., [[22], [24], [26]]). The main drawback of this technique is that it may not accurately capture the damage-induced anisotropy and failure mechanism governed by the morphological features of the masonry [[28]] unless a computational homogenization technique is adopted to capture failure localizations that stem from the mesoscale damage as demonstrated in the literature (see, e.g., [[29], [31]]). Thus, macromodeling is commonly used for the large-scale analysis of masonry buildings. The second technique, so-called micromodeling, is proposed to overcome these challenges by considering the masonry constituents (i.e., masonry units, mortar, and unit–mortar interfaces) in the mesoscale (e.g., [[32], [34]]). However, since this approach is computationally demanding with a large number of input parameters, it is only applicable to small-scale problems; hence, it is impractical for real-life large-scale buildings. The alternative and the most common approach, simplified micromodeling, is the midway between those techniques where the masonry units are explicitly represented as expanding up to the half thickness of the mortar joints, and mortar joints are indicated via zero-thickness interfaces (e.g., [[36], [38], [40], [42], [44]]). This approach reduces the required number of input parameters and provides less computational demand without compromising the detailed and explicit representation of essential failure modes in the URM walls (i.e., sliding and joint opening). Specifically, cracking in masonry units, joint opening and sliding at the unit mortar interfaces, and masonry crushing can be addressed in the simplified micromodeling, which will be discussed in the next section in detail. The well-known numerical modeling strategies used in micromodeling are the discrete element method (DEM) and discontinuous finite element analysis. Throughout this research, the DEM approach was utilized to simulate the in-plane behavior of a URM pier–spandrel system.

The lateral load-carrying capacity and failure mechanism of URM walls are primarily affected by material properties, boundary conditions, precompression stresses, geometrical properties, and wall cross section morphologies [[45], [47]]. In the case of URM pier–spandrel systems, the coupling effect between the pier and spandrel provides additional complexity to the problem along with the other parameters. Recently, pier–spandrel systems' behavior has been investigated experimentally and numerically to better understand the coupling effect, load-carrying capacity, and the associated collapse mechanisms [[9], [49], [51], [53]]. In this context, this research aims to strengthen the existing knowledge in the literature, offering a detailed look at the lateral load-bearing capacity of URM pier–spandrel structures under different precompression forces, including material uncertainty and its spatial variation throughout the structure. A simplified micromodeling approach, denoted as discrete rigid block analysis (D-RBA), was utilized to analyze the reference URM pier–spandrel system. The validated modeling technique was further utilized to conduct a comprehensive parametric analysis to quantify the influence of vertical precompression forces on the lateral load-bearing capacity and displacement behavior of the structural system. Three levels of precompression stresses, i.e., low, medium, and high stress (approximately corresponding to 0.2 MPa, 0.38 MPa, and 0.75 MPa, respectively), are considered. The variation in the response was quantified, and the results were compared to highlight the effect of the vertical precompression stress.

In this study, a DEM-based approach was utilized to simulate the structural behavior of URM pier–spandrel systems. The adopted modeling strategy was first developed by Cundall [[55]] to analyze rock mechanic problems and simulate the progressive collapse mechanism of jointed rock masses. In the late 1990s, DEM was used to assess the seismic behavior of unreinforced masonry arches and walls, where the stone skeleton of the structures was represented as a group of rigid blocks [[56]]. In the last several decades, DEM has become a strong alternative to the conventional finite element models and has been employed excessively in the structural assessment of masonry structures for static and time-history analysis [[48], [57], [59], [61], [63], [65]].

In DEM, the composite and discontinuous nature of masonry is represented by a system of rigid blocks that can interact with each other along their contact surfaces. It is worth noting that, within the DEM framework, the distinct block may be rigid and/or deformable depending on the desired outcomes from the computational model. Throughout this study, we used rigid blocks; hence, deformations were lumped at the contact surfaces. As mentioned earlier, the proposed DEM-based approach is denoted as a simplified micromodeling technique in the literature [[67]], in which the masonry units are expanded up to half thickness of the mortar joints, and mortar joints are indicated as zero-thickness interfaces as shown in Figure 1. The mechanical interaction between the discrete blocks was predicted using the point contact hypothesis, where each contact point had three orthogonal springs with limited tensile strength and cohesive shear behavior [[68]].

Overlapping is allowed among adjacent blocks that are controlled by the predefined contact stiffness in the normal ( $k_n$ ) and shear ( $k_s$ ) directions at each contact point (see Figure 2). Typically, the Coulomb–Slip joint model is utilized in the shear direction, which requires a cohesion and friction angle ( $\tau =c+tan\varphi \sigma$ , where $c$ is cohesion, $\varphi$ is the friction angle, and $\sigma$ is normal stress), whereas tensile stress is limited by the predefined strength value ( $f_T$ ). Throughout this study, fracture-energy-based contact constitutive models were utilized to better capture the postpeak response of the material in tension, shear, and compression. A linear softening behavior was adopted for all three cases, where tension–shear damage was coupled using a single damage parameter, as explained in [[69]]. In the adopted modeling strategy, each masonry unit is represented via two blocks (with a potential crack surface in between) to capture tensile failure (cracking) in the units. The readers are referred to the most recent study of the authors for a detailed explanation of the proposed approach denoted as discrete rigid block analysis (D-RBA) [[10]].

In DEM, the rigid block movements are computed by integrating translational and rotational equations of motion using the central difference method for each block to obtain new translational ( ${\dot{u}}_i^{t+}$ ) and rotational velocities ( ${\omega }_i^{t+}$ ), given in Equation (1) (velocities are evaluated at midtime intervals, i.e., $t^{+}=t+\Delta t/2$ , $t^{-}=t-\Delta t/2$ , $\Delta t$ : time step):

(1) $$\begin{gathered} {\dot{u}}_i^{t+}={\dot{u}}_i^{t-}+\frac{\Delta t}{m}\Sigma F_i^t-\lambda \Sigma F_i^{t}sgn{\dot{u}}^{t-} \\ {\omega }_i^{t+}={\omega }_i^{t-}+\frac{\Delta t}{I}\Sigma M_i^t-\lambda \Sigma M_i^{t}sgn{\mathsf{\omega }}^{t-} \end{gathered}$$

where $\Sigma F_i^t$ , $\Sigma M_i^t$ , $I$ , and $m$ represent the unbalanced force vector (i.e., the sum of the contact forces, self-weight, and external forces), the unbalanced moment vector (i.e., the sum of moments produced by contact and external forces), the moment of inertia, and the block mass, respectively. The quasistatic solutions were obtained from the equations of motions using Cundall's local damping formulation, using damping force (i.e., $F_d=\lambda \Sigma F_i^{t}sgn{\dot{u}}^{t-}$ ) and moment (i.e., $M_d=\lambda \Sigma M_i^{t}sgn{\mathsf{\omega }}^{t-}$ ). It is noted that the damping force and moment were calculated based on the magnitude of $\Sigma F_i^t$ and $\Sigma M_i^t$ , respectively, and applied as opposed to the motion ( $\lambda$ : a nondimensional force-based damping constant (default value is 0.8)) [[71]]. After calculating the new velocities ( ${\dot{u}}_i^{t+},{\omega }_i^{t+}$ ), the block positions were updated, and the corresponding action/reaction forces were obtained based on the relative block displacement increments in the normal ( $\mathsf{\Delta }{u}_n$ ) and shear ( $\mathsf{\Delta }{u}_s$ ) directions. The elastic contact stress increments were computed (i.e., the elastic normal stress increment $\mathsf{\Delta }\sigma =k_n\mathsf{\Delta }{u}_n$ and elastic shear stress increment $\mathsf{\Delta }\tau =k_s\mathsf{\Delta }{u}_s$ , with $k_n$ and $k_s$ being the normal and shear stiffnesses, respectively) and added to the old stresses (obtained in the previous step). If the new stress state at the contact point violated the implemented failure criteria, it was corrected and utilized in the next step. Finally, the new contact stresses were multiplied with the associated contact area to obtain the contact forces, later utilized in the equations of motion (Equation (1)). The explained computational procedure was performed in an explicit time-marching scheme until reaching the quasistatic equilibrium from the numerical solution. All DEM-based simulations were run using a three-dimensional discrete element code, 3DEC, developed by ITASCA Consulting Group, Minneapolis, CO, USA [[72]]. Numerical stability was ensured during the analysis by utilizing sufficiently small time steps ( $\Delta t<\Delta t_{cr}$ ). The critical time step ( $\Delta t_{cr}$ ) was automatically computed by 3DEC based on the minimum mass and maximum contact stiffness defined in the numerical model ( $\Delta t_{cr}=0.2\sqrt{m_{min}/k_{n,max}}$ ).

This section presents the proposed data preparation strategy and the reference study utilized as a benchmark. The authors already presented the validation of the computational model (see [[10]]); hence, only an overview of the test setup and analysis model are discussed in this paper.

In the benchmark study, a URM pier–spandrel structure built with tuff stones and lime mortar was subjected to in-plane pushover testing with 200 kN precompression loads kept constant on each pier and displacement-controlled lateral loading until reaching approximately 1% drift [[73]]. Each pier was 1.70 m long and 3.62 m high, whereas the opening between them was 1.70 m long and 2.30 m high. The spandrel that connected the piers to each other was 1.00 m high, with a timber lintel over the opening to support the masonry above. Identical boundary conditions and geometrical features were used in the numerical model where the lateral loads were applied to the cantilevered piers, as shown in Figure 3a. Lateral loads were prescribed as constant velocity boundary conditions to mimic the displacement control testing while the base shear was recorded during the analysis. As shown in Figure 3b,c, according to the analysis results, a similar collapse mechanism, including the rocking behavior of piers, flexural cracks at the pier–spandrel connections, and diagonal tension cracks at the spandrel, was obtained that was in line with the experimental findings. It is worth noting that the material and contact properties were taken from the reference study without conducting any calibration procedure. The readers are referred to the early publication of the authors for further details [[10]].

The adopted probabilistic analysis considered the variation in the joint and unit properties for masonry as random and dependent parameters. The random parameters were determined based on their effect on the strength and displacement capacity of the pier–spandrel system. These parameters are the compressive strength of masonry ( $f_{c,m}$ ), the joint (interface) tensile strength ( $f_{t,j}$ ), the tensile strength of units ( $f_{t,u}$ ), and the joint friction angle ( $\varphi$ ), as shown in Table 1. The experimental results [[73], [75]] were used to obtain the mean value and coefficient of variation (CoV) of the random parameters.

The dependent variables were presumed to be fully correlated with the random variables, such that their mean and CoV were the multiples of the statistical parameters of the random variables. As such, the cohesion of the units and joints was defined by the following relations: $c_{unit}=1.5f_{t,unit}$ and $c_{bond}=1.5f_{t,bond}$ , respectively. Similarly, the fracture energies were determined following the suggested ductility index values in [[76]]. Additionally, it is noted that the elastic stiffness of the lintel was considered constant during the analyses (adopted as $E_{lintel}=0.3$ GPa). A normal distribution was used for compressive strength and friction angle parameters, whereas a lognormal distribution was assigned to the tensile strength of units and joints, complying with the earlier work of the authors and the pertinent literature [[5], [7], [16], [77]]. The variables and their statistical distributions are presented in Table 1.

Once the statistical distributions and their parameters were defined, the Latin Hypercube Sampling (LHS) method [[79]] was used to obtain the sample values for Monte Carlo simulations. For each vertical pressure level, 250 simulations were run, i.e., adding up to 750 simulations in total. No correlation between the random variables was considered. The spatial variability of the parameters was set up in the following way: Each joint (or contact plane) in the structural system was defined with an identification number to assign the parameter values later. The identification numbers were defined separately for 367 units/blocks and 1211 unit–mortar interfaces to assign the mechanical properties effectively. In each of the 250 simulations, 367 values of unit tensile strength and 1211 values of joint tensile strength, friction angle, and masonry compressive strength were sampled using the statistical distributions and parameters given in Table 1. Once the contact properties were determined for each contact surface, either between or within the blocks, they were assigned to the associated contact points, providing homogeneous material property distribution through the contact plane. It is underlined that no spatial correlation was considered for the mechanical properties among the adjacent joints since it did not necessarily improve the understanding/quality of the probabilistic analyses, as discussed in [[41]]. The sampling process was repeated for 250 simulations. Once the parameters were generated for 250 simulations, the same set was used for all three vertical pressure levels. Three vertical stresses (or precompression loads), i.e., low, medium, and high, were determined to correspond approximately to 0.2 MPa, 0.38 MPa, and 0.75 MPa, respectively. Therefore, a total of 750 pushover analyses were performed.

The effect of the precompression stress on both the force and displacement capacities of the URM pier–spandrel system was investigated using the stochastic discontinuum-based model. The lateral force was applied up to a lateral displacement of 70 mm, corresponding to a 2.4% drift ratio. However, the force–displacement curves were truncated at 80% of the maximum lateral force once passing the ultimate load. The force–displacement behavior of the pier–spandrel system is presented as the relationship of the normalized displacement with the lateral load and with a seismic coefficient in the double vertical-axis plots. The normalized displacement represents the drift value, the ultimate displacement normalized by the height of the application point of the lateral loading (2.91 m from the ground level). On the other hand, the seismic coefficient is the lateral force normalized by the weight of the pier–spandrel system plus the vertical forces applied to the piers.

As expected and depicted in Figure 4a–c, the increasing precompression stress level increased the lateral load-carrying capacity. On the other hand, the seismic coefficient was between 0.42–0.48 for the low pressure, whereas it decreased to 0.34–0.37 for the medium pressure and 0.25–0.27 for the high pressure. The decrease stems from the increased vertical forces on the piers. The increase in the load capacity brought about a decrease in the lateral displacement capacity of the system, which is noted in Table 2.

The average maximum values of the peak lateral resisting force (Hmax) and their coefficient of variation are outlined in Table 2. It was observed that the variation in the maximum lateral load was much lower than the variability of the ultimate drift (θmax), complying with the existing literature and the authors' previous work on the same pier–spandrel system. For all cases considered, the variation in the maximum lateral load was less than 3%, indicating a highly reliable estimation of the load-carrying capacity. This phenomenon can also be observed in Figure 5, where the force–drift curves are given for both nonspatial and spatial analyses. Still, the ultimate drifts, considered at 0.8 Hmax on the postpeak descending branch of force–displacement curves, had significantly higher CoV, ranging between 23% and 54%, which highlights the noticeable variation in the displacement capacity. A different observation made in this study was the decrease in the variation in the average maximum lateral load when the spatial variation in the material properties was considered. This phenomenon is better observed in Figure 5a, where the response was visually more dispersed than the one in Figure 5b.

Another observation was the significant change in the drift response with the changing vertical precompression levels, see Table 2. The variability was much less for the high-pressure level, where presumably a more uniform failure pattern and force–drift behavior was observed, as seen in Figure 4c. On the contrary, the medium-level pressure results incorporated a highly variable response, also manifested in the CoV of maximum drifts.

Given that the failure modes were well separated, the statistical distributions of the load and drift capacities might be investigated and later used in a reliability assessment. As seen in Figure 6, for the high precompression loads where the behavior was more uniform, these distributions resembled beta and lognormal distributions. Further research is required to observe and quantify how the response parameters change together with their relationship with the failure mechanisms.

This study explored the lateral load-carrying capacity of URM pier–spandrel systems subjected to different precompression loads. The adopted computational model, validated earlier, explicitly represented the weak planes (i.e., unit–mortar interfaces) and masonry units as a system of rigid blocks and considered the uncertainty in the material properties. Moreover, the proposed modeling strategy offered a spatial variation in the strength properties of the mortar joints and units embedded in the contact points. There were 250 simulations performed for each precompression pressure, denoted as low, medium, and high. According to the analyses, the conclusions are summarized as follows.

Compared to the authors' previous study, where a nonspatial probabilistic analysis was performed, remarkably less variation in the average maximum lateral load was found considering the spatial variation in the material properties. The adopted spatial probabilistic D-RBA provided considerably less variation in the ultimate load (lower than 3%) as opposed to the ultimate drifts, which had a coefficient of variation ranging from 23% to 54%.

Furthermore, the effect of higher vertical pressure on the variation in ultimate drift (corresponding to 80% of the maximum lateral force) was noticed. Compared to medium and low precompression loads, more consistent force–displacement (or force vs. drift) responses were obtained with relatively less variation in the maximum displacement. This observation suggests that there should be less variation in the collapse mechanism in pier–spandrel systems when the vertical pressure gets higher. In future studies, the authors will further investigate this phenomenon by implementing automatic damage classification algorithms to exploit the type of collapse mechanism and the relationship between the applied vertical pressure. Overall, these conclusions strengthen the probabilistic analysis of URM structures and provide a detailed understanding of their structural behavior.

Graph: Figure 1 Illustration of the simplified micromodel of a regular masonry texture.

Graph: Figure 2 Representation of point contact hypothesis and spring-frictional element model.

Graph: Figure 3 Illustration of the loading pattern and results obtained from D-RBA: (a) representation of external loads, (b) collapse mechanism, and (c) failed contact points (tension combined with shear) along the unit–mortar interfaces (shown in red).

Graph: Figure 4 Force and seismic coefficient under varying drift for (a) low, (b) medium, and (c) high vertical pressures.

Graph: Figure 5 Force–drift curves obtained from (a) nonspatial analysis [[10]] and (b) spatial analysis of the present study (vertical force above each pier is 200 kN).

Graph: Figure 6 Ultimate force (Hmax) and drift (θmax) distributions for high precompression analyses, illustrated in (a) and (b), respectively.

Table 1 Variables, probability distributions, and statistical parameters.

Table 2 Average maximum lateral force and drift for each precompression level.

Software, writing—review and editing, investigation, and methodology, B.P.; writing—original draft preparation, methodology, and conceptualization, S.G.; conceptualization, methodology, investigation, and writing—review and editing, F.P. All authors have read and agreed to the published version of the manuscript.

Not applicable.

The authors declare no conflict of interest.