Stochastic strength prediction of masonry structures : a methodological approach or a way forward ?

Today, there are several computational models to predict the mechanical behaviour of masonry structures subjected to external loading. Such models require the input of material parameters to describe the mechanical behaviour and strength of masonry constructions. Although such masonry material parameters are characterised by stochastic‐probabilistic nature, engineers are assigning the same material properties throughout the structure to be analysed. The aim of this paper is to propose a methodology which considers material spatial variability and stochastic strength prediction for masonry structures. The methodology  is  illustrated on a case study covering the  in‐plane behaviour of a  low bond strength masonry wall panel containing an opening. A 2D non‐linear computational model based on the Discrete Element Method (DEM) is used. The computational results are compared against those obtained from the experimental findings in terms of failure mode and structural capacity. It is shown that computational models which consider the spatial variability of masonry material properties better predict the load carrying capacity, stiffness and failure mode of the masonry structures. These observations provide new insights into structural behaviour of masonry constructions and lead to suggestions for improving assessment techniques for masonry structures.


Introduction
Masonry is heterogeneous material in nature which consists of units (i.e. clay/concrete bricks, stones etc.) and mortar joints. Even in the same structure, the mechanical properties of masonry units and mortar vary. Such variability affect the mechanical response of masonry constructions [1,2].
Variability in masonry materials have been studied by several researchers in the past. For example, Melbourne et al. [3] discussed that bricks fired at higher temperatures in the kiln present better quality and durability characteristics to those fired at lower temperatures. Moreover, stone inclusions found in the raw clay materials can also vary the mechanical properties of clay masonry bricks, in some cases increasing their compressive strength. Also, variation in ageing and deterioration of bricks even in the same structure is another factor which influence their condition and physical characteristics.
Variations in mortar could arise due to: a) the composition and quality of the materials used during the mortar making process; b) the interaction of the mortar with the adjacent masonry units; c) the orientation of the mortar joints in the structure; and d) material degradation as a result of the ageing effects. According to Brocken et al. [4], mortar between masonry units with high moisture content will cure differently to those between units with lower moisture content. In addition, the characteristics of bed and head/perpend mortar joints could be different. According to Dialer [5], the strength of the head joints is usually lower than the strength of the bed joints. This is a result of the greater degree of mortar shrinkage in the perpend joints. Also, perpend mortar joints are often not filled fully with mortar. According to Mann and Muller [6], changes of the strength and stiffness of the perpend joints can generate a discontinuous or non-uniform stress distribution in masonry. Moreover, workmanship effects such as incomplete filling of mortar joints could introduce variability in unit-to-mortar bond strength [7][8][9][10]. Such defects could be presented in varying degrees and the global behaviour of the brickwork will reflect their combined effect. However, assessing the overall effect of workmanship on the behaviour and strength of brickwork is not a straightforward issue [1,7].
Today, there are several computational tools to predict the inservice and load carrying capacity of masonry. These range from the classical plastic solution methods [8] to the advanced non-linear computational formulations (e.g. finite and discrete element methods of analysis). Computational models require the input of material properties to describe the mechanical behaviour and strength of masonry constructions. These must be selected carefully by the modeller/engineer.
Although masonry units and mortar are characterised by variability in strength, engineers are typically use non-spatial models (e.g. same material properties throughout the masonry domain) for estimating the mechanical response of masonry constructions [9]. It is only recently that some research work has been done to address material variability in numerical models for the analysis of masonry constructions.
For example, Fyfe et al. [10] investigated, with the use of Monte-Carlo simulations, the effect of bond tensile strength and bed joint thickness variation on the mechanical response of laterally loaded masonry walls. Fyfe et al. used a 3D macroelement software based on the FEM. A mean failure load and a safety factor were obtained. A few years later, Li et al. [11][12][13] using a 3D FEM model and Monte-Carlo simulations, predicted the lateral resistance of masonry walls subjected to one way horizontal bending taking into account the unit-tounit spatial variability. Results demonstrated that models with spatial variability of material properties better represent the experimental behaviour of masonry walls. Similar findings were also obtained from Moradabadi et al. [14], who used a macro-model developed in ABAQUS software to predict the maximum compressive strength of masonry prisms and unreinforced masonry walls in two way bending [12]. Zhu et al. [15] studied the mechanical behaviour of compressed masonry wallets made from hollow concrete elements. The compressive and tensile strength of masonry obtained using a Monte-Carlo type simulations and a Latin Hypercube sampling approach.
From the above, although some work has been done in the past to address material variability in numerical models, such studies are mainly focusing on the use of macro-elements based on the finite element method combined with Monte Carlo simulations. In addition, there is still no systematic procedure to consider the spatial variability of masonry materials (i.e. bricks and mortar joints) for the analysis of low bond strength masonry constructions. This paper presents a methodological approach able to incorporate the spatial variability of mortar strength properties in low bond strength masonry constructions. A numerical model based on the Discrete Element Method (DEM) of analysis has been used. The model was initially developed to model rock sliding in which failure occurs along the joint. This has similarities with the behaviour of low bond strength masonry in which the predominant failure mode is at the brick to mortar interface. The methodology proposed herein is applied on a case study covering the in-plane strength of a low bond strength masonry wall panel with opening [16][17][18]. Computational results are compared against experimental ones and the findings are discussed.
The remainder of this manuscript is structured as follows: Section 2 discusses the methodology of the probabilistic analysis. Section 3 presents an overview of the discrete element method for modelling masonry. Section 4 describes the experimental procedure and results from testing the masonry wall panels with opening, and Section 5 presents the details of the development simplified micro-model. The stochastic analysis of the spatial and no-spatial model are presented in Section 6 and compared against the experimental findings. Conclusions and recommendations for future works are outlined in Section 7.

Methodology
The approach commonly used by engineers/modellers to ignore the spatial variability of masonry in their computational models is found to be problematic. An alternative method which reflects the heterogeneity of masonry is proposed here. According to the method, the spatial variability of material parameters is considered as input for the numerical model (i.e. material parameters with a known probability distribution). Distribution properties of masonry material parameters can be determined by performing a series of in-situ small scale tests on masonry constituents e.g. direct tension, direct shear, compressive strength of bricks and mortar etc [19]. Computational models with stochastic spatial material parameters can then be carried out. Realizations of stochastic material variables can be randomly assigned to the model (hereafter called "spatial models"). The output of the simulation (e.g. the ultimate load or the load at first crack) will be in a form of a probabilistic distribution. The proposed methodology is shown in Fig. 1.

Overview of the Discrete Element Method and UDEC software
In this study, the commercial 2D software UDEC (Universal Distinct Element Code) based on the Discrete Element Method (DEM) has been used. Within UDEC, masonry units are internally discretized into triangular finite elements, which assumed to behave according to a linear elastic stressstrain law; since the predominant failure mechanism in low bond strength masonry is in mortar joints and not masonry units. Masonry units are separated by zero-thickness interfaces to represented mortar joints. At the interfaces, the bricks are connected to each other by sets of point contacts. These contact points are located at the outside perimeter of the masonry units and are created at the edges or corners of the bricks and at the finite element nodes. As the number of point contacts in the model increases, so does the accuracy of the stress distribution obtained. The mechanical behaviour of the contacts is shown in Fig. 2. The contact model follows the soft-contact approach: Joint stiffnesses are defined in the normal (kn) and in the shear (ks) direction. Small interpenetration between neighbouring blocks in contact are allowed. The contact stiffnesses represents the deformability of the mortar. In this study, the normal stiffness has been taken according to Equation 1: where k n is the normal stiffness of the joint, while E m the elastic modulus of the mortar and t m is the thickness of the mortar joints. The maximum allowable tensile stress in the mortar is limited by the tensile strength parameter (f t ). If the normal stress exceeds the tensile strength, the normal stress is set to zero. The shear behaviour follows the Coulomb-law. If the shear stress exceeds the maximum shear stress, the internal friction angle (φ) and the cohesion (c) is set to their residual values (φ res , c res ). Non-associative flow rule was applied; hence the dilation angle (ψ) was set to zero.
In the model, the Newtonian equations of motion were solved with explicit time integration, using the central difference method. To reach static equilibrium artificial damping forces was introduced in the system. These viscous forces were decreased as the system reached its equilibrium state. For further information about the DEM, the reader is kindly directed to [20].

Experimental testing of the masonry wall panels with opening
In the laboratory, four masonry wall panels (S1 to S4) were tested. Each panel contained approximately 2 m opening (Fig.   3). All panels had a soldier course immediately above the opening. This is a typical architectural feature of façades in the UK. The dimensions of the bricks were 215 mm × 102.5 mm × 65 mm. All mortar joints were approximately 10 mm thick. Mortar was made from 1:12 (OPC:sand) weigh-batched mortar. The bricks and mortar were selected to produce brickwork with a low bond strength characteristic. A single vertical point load applied at mid-span and at the top of the panel. Load applied incrementally and mid-span deflections were recorded at each load increment. Cracks of approximately greater than 0.15 mm in width (i.e. hairline cracks) were monitored and recorded.

Development of the computational model based on the Discrete Element Method
A geometric model of the masonry wall panel tested in the laboratory (Fig. 4) was created in the discrete element model. Within the model, each brick was represented by a deformable block separated by a zero-thickness interface at each mortar bed and perpend joint. The mortar joints were modelled using elastic-perfectly plastic coulomb slip-joint area contact option [20]. In the present work, normal distribution of the material parameters for the masonry constitutive model has been assumed. The mean values and the coefficient of variance (CoV) were obtained from the literature [1,19]. The material properties used in the computational model are shown in Table 2. The joint normal stiffness and shear stiffness were related to each other according to the following expression [16]: The tensile and cohesive strength of the mortar joints were also related to each other using the following equation [19]: A series of computational models with non-spatial and spatial variability were developed. In the former case, a realization for a material parameter was created and assigned to every brick and mortar joint in the masonry wall panel. In latter case, models with spatial variability were designed to capture the range of strength properties of mortar joints. In particular, a realization of the material parameter was independently created for every brick and mortar joint. In this study, no attempt was made to correlate the material parameters of the neighbouring bricks and mortar joints.
In the computational model, the bottom edges of the wall panel were modelled as rigid supports whilst the vertical edges of the wall panel were left free. Initially the model was brought into a state of equilibrium under its own self-weight. Then, a constant vertical velocity was applied at the load spreader plate on the top of the wall panel. Convergence tests were also carried out on the magnitude of the velocity to be applied to the spreader plate to make sure that a quasi-static loading condition was achieved. The velocity was converted to a vertical displacement and the force acting on the spreader plate for each load increment was calculated with the use of a FISH function in UDEC. Hence, a load versus midspan displacement relationship for each material realisation in the masonry wall panel was determined. Results from these were then compared against experimental findings.
To investigate the influence of the different material properties on the global behaviour of the masonry wall panel, four sets (Table 3) of simulations were carried out by considering only a subset of the properties as stochastic variables, while the other ones were handled as deterministic variables. Mean and CoV values were assigned to the deterministic parameters. Each set of analysis contained 300 spatial and 300 non-spatial numerical simulations. The number of simulations was set high enough to have confident results from a statistical point of view. The total computational time for a single set of simulations (300 simulations) was approximately 48 hours. The difference between the defined mean value of a masonry material property and the average of the first n realizations can be seen in Fig. 5. Fig. 6 shows the difference between the prescribed standard deviation and the deviation of the first n realizations. The difference from the prescribed average and standard deviation value after 300 simulations is less than 1.7% and 1.8% respectively.

Results and discussion
In this section, the influence of material variability on the mechanical behaviour of the masonry wall panel were investigated. Monte-Carlo type of analyses were performed to determine the effect of variation of masonry material properties in the masonry wall panel with respect to the load at which first visible crack (crack width 0.15 mm) develops and ultimate load carrying capacity (i.e. maximum load that the wall can carry). The failure mode of the masonry wall panel has also investigated and compared with the experimental findings.
Both experimental and computational results showed that the behaviour of the masonry wall panels is characterised by an initial cracking in the soldier course of the panel followed by flexural cracks in the bed joint at the supports. As the load in the panel increased, diagonal stepped cracks appeared from the corner of the opening towards the mid-span, until the panel failed due to shear (or excessive diagonal tension). Images of a typical masonry wall panel tested in the laboratory can be seen in Figure 7. It should be noted that in the case of a computational model representing a symmetric masonry structure with a symmetrically applied load and nonspatial variability in material properties, a symmetric failure mode would occur. However, in the experiment, this is not the case due to inherent variability of masonry material parameters, imperfections in construction of the wall, application of loading etc.   Fig.8, all spatial models failed with an asymmetric failure mechanism. Failure mode #2 is characterized by the cracks arising from the bending of the soldier course. Above it, diagonal cracks can be seen from the application point of the external load to the corner of the opening. Tensile crack of the bed joints near to the opposite can be seen. Failure mode #1 and #3 include shear failure of the soldier course at the corner and at the middle part of the opening. Failure mode #4 contains significant diagonal cracks from both corners to the external load, while shear failure of the soldier course occurs typically at one of the corners.
Failure mode #5 is characterized by the shear failure at or near above the base of the model. In case of a non-spatial symmetric models the failure mode is expected to be symmetric as well. The distribution of the different types of failure mechanisms as obtained from the 300 simulations for the material variability as per Set #1 (See Table 3) are shown in Fig. 9.   Fig. 11 compares the numerical against the experimental load-displacement relationships in the masonry wall panel in case of spatial and non-spatial analysis, respectively. The numerical load-displacements curves were obtained from the results of Set 1 (see Table 3), when all of the material properties were taken as stochastic parameters. The experimental curves stop at a point before ultimate load has been reached. The little peaks in the curves shown in Fig. 10 represent relaxation of the loading and moment re-distribution in the panel due to the formation of a new crack. When a crack propagates, there is an abrupt loss of stiffness in the panel. From Fig. 10 and Fig. 11, it is evident that the spatial models give a better approximation of experimental load-displacement curves both in terms of the ultimate load, deflection at failure and stiffness.   In the numerical model, the load at first visible crack was recorded when the normal displacement at a mortar joint (i.e. interface) was equal with 0.15 mm. Fig. 12 shows the distribution in the form of a histogram of the load at first visible crack as predicted by the numerical model for the Set 1 material variables. From Fig. 12, the deviation of spatial and non-spatial models in case of the load at first crack is very similar. However, the distribution of the load at first crack obtained from the computational models was lower than the one obtained from the experimental study. This was due to the fact that the experimental load at first crack captured when the crack width was approximately 0.15 mm, and not exactly 0.15 mm. Fig. 13 shows the histogram of the ultimate load bearing capacity of the masonry wall panel as predicted by the DEM model and using the material paramters from Set #1. From  Fig. 13, the predicted mean value of the ultimate load (4.73 kN for the non-spatial, and 4.44 kN for the spatial analysis) of the masonry brickwork wall panel compares quite well with the range of values obtained from the experiments (see Table  1, mean 4.76 kN). The ultimate load shows significantly higher deviation in case of non-spatial analysis. The smaller standard deviation of the spatial model can be explained with the averaging of the neighbouring material properties in the computational model. In the case of non-spatial simulations, if a realization of a material property gives an extremely small/high value, this value is assigned to every joint/brick of the model causing extremities in the value of ultimate load as well.  Table 3). Figure 13. Histogram of the ultimate load that the panel can carry as predicted by the numerical model (for the Set #1 material variables, Table 3). Fig. 14 shows the correlation between the load at which first crack occurs and the ultimate load that the masonry wall panel can carry with respect to material properties with spatial and non-spatial material variability. The computational models with and without spatial material variability do not show correlation between the load at which first crack occurs and the load carrying capacity.  Table 4 shows the mean and coefficient of variance (CoV) of the ultimate load bearing capacity that the masonry wall panel can carry for the different set of materials parameters ( Table 3). The sensitivity of the model was analysed with different sets of simulation, where some of the input parameters were handled as deterministic ones, while others were considered as stochastic parameters. From Fig. 15 and Fig. 16, the ultimate load bearing capacity is significantly influenced from the mortar strength parameters (i.e. tensile strength and cohesion of mortar), while the effect of brick and mortar elastic parameters has less significant effect on it. Moreover, for the investigated range of the strength parameters, the relationship between the ultimate load and the cohesive and tensile strength appears to be linear. When the strength parameters were handled as stochastic (Set 3, Table 3 and 4), the standard deviation of the ultimate loads were significantly higher in case of computational models with spatial variability in material properties. On the contrary, when the elastic parameters were handled as stochastic (Set 2 and Set 4, Table 3 and 4) the spatial models showed higher values of the coefficient of variance (CoV).

Conclusions
Masonry constructions are characterised by inherent variability. Such variability arises from: a) the masonry units; b) the mortar joints; c) the workmanship effects; and d) weathering effects in the structure. Although it is well known that masonry material properties are characterised by spatial variability even in the same structure, engineers are often assigning the same material properties throughout the structure to be analysed. The aim of this paper is to highlight the material parameter identification problem and proposes a methodology to incorporate the variability of masonry unit and mortar joint properties into a simplified micro-model based on the discrete element method of analysis.   A 2D non-linear discrete element model of a masonry wall panel with opening has been developed. The effect of spatial variability of material properties on the mechanical behaviour of the masonry wall panel investigated. Numerical results of the model which included spatial and non-special variability of material properties were compared against experimental findings. From the results analyses, it was shown that the model was able to capture the ultimate load and failure mode of the masonry wall panel. Also, the model with the spatial variability of material properties gives a better representation of the behaviour of the real structure when compared to the case of non-spatial analyses. The elastic properties of the brick do not influence the load carrying capacity of the masonry wall panel. On the other hand, the strength parameters of the mortar joints are driving the behaviour and load carrying capacity of the wall. Cohesive and tensile strength are the predominant parameters which influence the load carrying capacity of the wall. In the case of non-spatial analysis, high deviation of the load carrying capacity of the wall panel observed. However, for spatial analysis, a decreased deviation was observed. This was a result of the averaging effect of the stronger and weaker joints in the masonry wall panel. In the future, further analysis will be carried out to investigate the material variability problem for bricks and mortar joints in different masonry constructions.