Experimental, numerical and analytical investigations of masonry corners: Influence of the horizontal pseudo-static load orientation

Construction and Building Materials 344 (2022) 127969 0950-0618/© 2022 Published by Elsevier Ltd. Experimental, numerical and analytical investigations of masonry corners: Influence of the horizontal pseudo-static load orientation Carla Colombo a , Nathana¨ el Savalle a , Anjali Mehrotra b , * , Marco Francesco Funari a , c , Paulo B. Lourenço a a University of Minho, ISISE, Guimar˜ aes, Portugal b TU Delft, Delft, the Netherlands c Department of Civil and Environmental Engineering, University of Surrey, Guildford GU2 7XH, UK ARTICLE INFO Keywords: Limit analysis Flexural mechanism Rocking-sliding mechanism Seismic assessment Tilting tests Dry-joint stiffness Discrete element modelling ABSTRACT The present work aims to expand the knowledge of the behaviour of masonry corners, which are capital to obtain an integral seismic response in masonry buildings. In particular, the influence of the seismic load orientation (from π /4 to π /2) is investigated experimentally, numerically and analytically. Both units and interfaces have been subjected to a material characterisation process, following which pseudo-static 1:4 scaled experiments on a tilting table have been conducted on a symmetric dry-joint masonry corner. The experimental results have also been simulated using a discrete element model. Finally, a new analytical limit analysis model has been developed, which considers both experimental and numerical observations and accounts for rocking-sliding and flexural mechanisms. In general, a good agreement is found between the three approaches, both in terms of collapse mechanism and load multiplier. 1. Introduction Simulation of the seismic behaviour of historical masonry structures is a challenging task. It depends on several factors, such as the mechanical properties of the constitutive materials (units and mortar), the geometrical configuration of the structure, the connections between load-bearing components, the stiffness of the horizontal diaphragms, etc. When an earthquake occurs, localised out-of-plane (OOP) failure mechanisms are prone to form in historical masonry structures, particularly in the absence of box-like or integral structural behaviour [1,2]. In this framework, following post-seismic damage surveys carried out after the Irpinia and Syracuse earthquakes in Italy, in [3], an abacus of local failure mechanisms was compiled that may be assessed through simple analytical formulations mainly based on the theorems of the limit analysis (LA). Once a mechanism is selected, a set of equilibrated generalised forces and a set of compatible generalised virtual displacements are determined. Then, the work done by the generalised forces in equilibrium with the internal stresses for the given set of generalised virtual displacements is computed. The application of the static theorem leads to a lower-bound or safe solution based on equilibrium equations, while the application of the kinematic theorem provides an upper-bound multiplier of the collapse load factor [4]. Thus, the solution that satisfies the hypotheses of the pre-mentioned theorems and equilibrium, compatibility and material conditions is the correct solution and provides the theoretical actual collapse load multiplier for the specific problem. In the last decade, following Heyman’s pioneering work [4], to investigate the most reasonable collapse mechanisms, different algorithms were implemented into user-defined routines of analysis, adopting the lower or upper bound theorem of LA [5–8]. Among the damage mechanisms found in masonry structures, failure of masonry corners is commonly observed during earthquakes [9,10]. Certain geometrical features or structural details usually increase the probability of corner failure, among them: i) the lack of confinement of the building by adjacent structures, ii) the presence of openings, iii) the presence of chimneys, and iv) the lack of rigid horizontal diaphragms. However, limited investigations have been conducted concerning the seismic vulnerability of masonry corners. To that end, in [5], an analytical formula based on LA was proposed to determine masonry corners’ associated horizontal load multiplier. A similar methodology * Corresponding author. E-mail addresses: [email protected] (C. Colombo), [email protected] (N. Savalle), [email protected] (A. Mehrotra), m.funari@surrey. ac.uk (M.F. Funari), [email protected] (P.B. Lourenço). Contents lists available at ScienceDirect Construction and Building Materials journal homepage: www.elsevier.com/locate/conbuildmat https://doi.org/10.1016/j.conbuildmat.2022.127969 Received 15 February 2022; Received in revised form 16 May 2022; Accepted 24 May 2022 Construction and Building Materials 344 (2022) 127969 2 was also applied in [11,12], which included the actual value of the frictional resistance defined as a weighted value dependent on the crack line angle. In accordance with post-earthquake surveys, these models also account for openings, which are a significant factor of weakness for this mechanism [10]. Recently, LA of masonry corners was also used as the first step of a workflow to identify the most likely collapsing macroblock’s geometry [13]. The second step consisted of a simplified analytical dynamic rocking analysis of the identified macroblock [13]. Such an approach leads to an overall efficient and engineering-oriented tool, similar to what has also been proposed for masonry façades [14]. However, while most studies focused on the rocking-sliding mechanism of masonry corners, the experimental tests presented by [15] evidenced through tilting tests that a flexural mechanism can also develop. Consequently, a separate LA formulation for this flexural mechanism was proposed and calibrated against experimental results [15]. Nevertheless, the failure of corners was only investigated when loaded symmetrically with respect to both walls (i.e. π /4), thereby revealing a gap in the literature. In general, LA-based tools neglect the structure’s global behaviour and only focus on assessing pre-defined local failure mechanisms [3,5,7]. Nonetheless, in many cases, there is a need to investigate the non-linear global behaviour of masonry structures in more detail, and this may be achieved by adopting sophisticated numerical methodologies, such as the Finite Element Method (FEM) [16–19] or the Discrete Element Method (DEM) [20–24]. Such approaches model masonry using different representation scales: equivalent continuum, macroblocks, or discrete representations (i.e. micro-modelling). DEM, in particular, enables the explicit representation of the units and joints [21,25], thereby making it well-suited for modelling the response of blocky masonry structures. DEM has been applied with success to the analysis of both unreinforced and retrofitted, dry-joint and mortared, masonry wall assemblies [22,24,26–33]. Despite the accuracy, the computational efficiency of the available numerical methods is rarely compatible with the need to have a rigorous real-time post-earthquake assessment [34]. Thus, researchers are committed to developing alternative modelling approaches and practical tools to decrease the computational cost without compromising the accuracy [35–38]. As far as masonry built cultural heritage is concerned, joint behaviour can broadly be akin to dry-joint behaviour. Indeed, if present, it is common that mortar deteriorated through time losing its bonding properties. Therefore, this paper focuses on dry-joint masonry specimens. Single-leaf and scaled (1:4) masonry corners tested on a tilting table are adopted as they provide a relatively simple validation framework for numerical and analytical models. The objectives of the present work are summarised below: •Provide experimental data to validate numerical and analytical methodologies to study the stability of masonry corners. •Provide an illustrative example of the complete characterisation of dry-joint masonry specimens at low normal stress, including dryjoint stiffness. •Quantify the effect of the seismic orientation on the failure mechanism and associated collapse load multiplier of masonry corners. •Improve existing analytical LA models to account for the orientation of the seismic loading in an integrated way. •Validate LA predictions of the collapse load multiplier and failure mechanism against experimental and numerical DEM results. The experiments are reproduced numerically using the Discrete Element software 3DEC, while a novel analytical LA formulation for masonry corners is also developed. Such a general formulation allows to consider, jointly or separately, rocking-sliding and flexural mechanisms noted by [15] and accounts for the orientation of the seismic loading. In short, this work investigates the influence of the orientation of the seismic load (adopting static lateral loading) on the capacity and failure mechanism of dry-joint masonry corners. The paper is organised as follows. Section 2 describes in detail the experimental campaign that concerns the mechanical characterisation of both the blocks and interfaces and the tilting tests of the masonry corner prototypes. Section 3 is devoted to the DEM model and its validation against the experimental results. Experimental and numerical observations are subsequently adopted to formulate a new analytical formulation based on LA in Section 4. Finally, relevant conclusions are drawn in Section 5. 2. Experimental campaign In this section, the experiments involving single-leaf dry-joint masonry corners are presented. The campaign studies the influence of the seismic loading direction on masonry connections and constitutes a reference for validating analytical and numerical results. Furthermore, the investigation describes the characterisation of the material, including unit and interface properties. 2.1. Material characterisation The masonry corners were built of dry calcium silicate (CS) blocks, Fig. 1. Experimental setup sketches: (a) normal joint stiffness evaluation by means of joint closure tests and (b) tangential joint stiffness evaluation through direct shear box tests. C. Colombo et al. Construction and Building Materials 344 (2022) 127969 3 provided by the Xella Company [39], which were sawn to dimensions t×l×h=57.3±0.6(1.1%)×114.9±0.5(0.5%)×70.3±0.2(0.3%)mm3 and density equal to ρ =1876±17(0.9%)kg⋅m-3. Here, the standard deviation and coefficient of variation (CoV) are also given. The compressive strength fcr of the units was evaluated, following the ASTM D 2938-95 [40] standard, through five uniaxial displacement-controlled compression tests on prismatic bricks. The tests provided a value equal to fcr = 36.5±2.4(6.6%)MPa. The Young’s modulus was subsequently obtained through ten compression tests along the three orthogonal material directions [41]. Besides a small anisotropy due to the industrial vibrocompaction production process, a value of 7.0±1.16(16.5%)GPa was found to characterise the material. Joint normal kn and shear ks stiffnesses describe the variation of the normal (or shear) stress with the normal (or shear) displacement. They constitute critical numerical parameters in discrete element modelling approaches, though no clear guidelines already exist to support their identification. Herein, the normal joint stiffness kn was characterised through classical joint closure tests [42–46], which consisted of uniaxial compression tests on two-stacked CS blocks (Fig. 1a). The displacement reading was obtained using four monitoring LVDTs placed around the blocks in the vicinity of the dry-joint, supported by two aluminium rings bounding the extremity. The tests were displacement-controlled through an additional LVDT measuring the actuator movement and consisted of five loading-unloading cycles. At the beginning of each loading cycle, the joint was completely unloaded to obtain stress-displacement relationships for very low stresses, here up to 0.01 MPa. Fig. 2a shows an example of the obtained experimental curves, in which the displacement due to block elastic deformation has already been removed based on the previously evaluated Young’s modulus. After an initial phase of joint closure, where relatively low stresses trigger large displacements, the slope (i.e. the joint stiffness) dramatically increases for larger stress levels. Each experimental curve has been fitted using an exponential function (Fig. 2a), which expresses the stress as a function of the joint closure dj and two empirical constants A and B as follows [42]: σ =A×eB×dj(1) Differentiating both sides of Eq. (1) with respect to the normal joint displacement dj, the normal joint stiffness reads: kn=d σ ddj=B× σ (2) From the 25 tests conducted, two outliers have been disregarded, and the mean value of the empirical constant B was found equal to B= 195.6±49.4(25.3%)mm-1. It is noted that the coefficient of variation of this joint property is much larger than the ones found for the mechanical properties of the unit (6.6% for compressive strength and 16.5% for Young’s modulus). This difference is explained by the larger variability of the physical properties of the joint (e.g. surface’s roughness, geometrical tolerance) which directly influence the joint stiffness compared to volumetric properties. Direct shear tests were also conducted to evaluate the tangential stiffness k s of the dry-joints. The specimens consisted of two CS blocks of 55.82 ±1.17 (2.1%)mm length, 57.65 ±0.27 (0.5%)mm width and 27.15 ±0.81 (3.0%)mm height. The test setup was composed of i) a bottom block located within the box base and pushed by a horizontal force FH and ii) a fixed top block subjected to a constant normal force FV (Fig. 1b). The data acquisition system comprised two load cells measuring the vertical and horizontal forces FV and FH, while the relative normal and shear displacement between the bottom and the top blocks was monitored by a Digital Image Correlation (DIC) system, which has already been employed in the literature for extracting the interface stiffness, proving its efficiency against other types of acquisitions [47]. This non-contact full-field measurement method solely identifies the relative displacement between the two blocks, whereas the shear box displacement acquisition system, i.e. vertical and horizontal LVDTs, also measures the external compliances owing to the shear box rig. Fig. 2b provides an example of the obtained shear stressdisplacement curve assuming a pre-compression of 10 kPa. The first window (I) corresponds to the so-called “elastic-stick” regime, in which the increase of shear stress leads to very small relative displacements [47,48]. In the second window (II), micro-slips progressively develop: at the micro-scale, some asperities are in adhesion (i.e. stuck), while other asperities experience relative tangential displacements. Such a phase is a transition from window I to window III, which corresponds to the gross slip phase, where all the asperities experience relative tangential displacement: the value of shear stress reaches the shear strength of the joint. Ten tests were carried out and no dilatancy has been noticed, similar to [49]. The tangential stiffness ks was estimated by considering the slope of the elastic-stick phase (I), where the stress grows with a linear trend with respect to the displacement. Since the vertical stresses involved in the Fig. 2. (a) Fitting process to evaluate the stress dependency of joint normal stiffness and (b) example of the stress-displacement curve obtained from direct shear tests. C. Colombo et al. Construction and Building Materials 344 (2022) 127969 4 experiments lie between 0 kPa and 10 kPa, the shear tests assumed 5 kPa and 10 kPa for the vertical stresses. Five tests for each precompression level were performed, obtaining a mean value of ks= 0.33 ±0.07 (22.8%)GPa⋅m−1 at 5 kPa and ks=0.49 ± 0.25 (51.1%)GPa⋅m−1 at 10 kPa. Again, a high coefficient of variation is obtained for the interface stiffness. Moreover, it is noted that the ratio between shear and normal stiffness for these two levels, using the previous expression for the normal stiffness, is 0.34 and 0.25, which is also in line with results found by others [50,51]. Finally, the shear (III) strength was extracted from the shear box tests considering 11 normal loading conditions, i.e. from 3 kPa to 50 kPa [52]. Assuming a Coulomb behaviour for the joint, the friction coefficient was estimated through linear regression of all the experimental outcomes and results in μ =0.67 (R2=0.97). The material characterisation process was completed by identifying the friction coefficient also through pseudo-static tilting tests, which was then compared with the results from the direct shear tests. Two bottom blocks were fixed to the tilting table, while a top unit was left free to slide. The table was then progressively inclined to obtain the block-block friction coefficient, while a digital inclinometer attached to the table allowed the reading of the sliding angle. A total of 50 tests were performed, leading to a mean value of μ =0.66 ±0.06 (9.5%), with a range between 0.50 and 0.78. This value is close to the one identified through direct shear tests, especially when considering the experimental scatter. The obtained variability can be related to the surface’s dustiness and also to some geometrical effects due to the presence of two interfaces instead of one. The mean value μ =0.66 is hereafter assumed globally representative of the dry-joint behaviour of the CS blocks, given the agreement with the shear box tests and the fact that this testing condition (i.e. layout of the joint) is physically closer to the corner’s test setup. 2.2. Test setup Tilting table tests have often been adopted as a tool to study collapse mechanisms of masonry structures [15,53,54]. A simple apparatus approximates the seismic forces through the progressive tilting of one extremity of the table: seismic forces are experimentally represented as static horizontal actions proportional to the specimen mass. Despite their simplicity, tilting tests estimate the horizontal collapse load multipliers of given structures and make it possible to observe the associated failure mechanisms [54]. The present experimental campaign focuses on studying the capacity of masonry corners to withstand horizontal forces with respect to seismic load orientation. Experimentally, the building-up process varied the direction of the corner bisector plan towards the rotational axis of the hinge (Fig. 3). More in detail, five load orientations equally spaced by π /16 were analysed, i.e. π /4,5 π /16,3 π /8,7 π /16 and π /2. The specimen consisted of a single-leaf dry-joint masonry corner, which comprised two walls of, ideally, 804 mm in length and 633 mm in height (Fig. 3 and Fig. 4). The courses had seven full blocks for each wall, staggering the joints between the layers and reaching nine horizontal courses (Fig. 4a). A total of 25 tests (5 per configuration) were performed at the laboratory of the University of Minho in Guimar˜ aes (Portugal). The adopted testing rig consisted of a 1.5m ×1.5 m steel table welded on a double T frame. The table was connected to an actuator through a steel wire that allows the uplifting of one extremity, while the other end rotated around a hinge [56]. Tilting speed was set to 0.1◦⋅s−1 before an angle of 9◦and 0.02◦⋅s−1 afterwards, in order to ensure quasi-static conditions close to the specimen collapse. Sliding between the first course and the steel table was prevented by fixing this course to the table with a double-faced tape, while the uprightness and repeatability of the building up process were ensured using a timber guiding system. The free ends of the walls were clamped with an L-shaped profile and rubber pads to simulate ideal boundary conditions (Fig. 4b). A digital inclinometer installed on the table captured the rotation angle, while the collapse mechanism was recorded by two digital camFig. 3. Top view sketch of the corner position by varying the bisector direction. Fig. 4. (a) Corner 5 π /16 test setup and (b) side view of the clamping system. C. Colombo et al. Construction and Building Materials 344 (2022) 127969 5 eras placed orthogonally to the plane of each sidewall. Finally, three lights were installed in the vicinity of the corner to allow a recording speed of 10 Hz with a shutter speed of 4 ms. 2.3. Experimental results In this section, the outcomes obtained from the experimental campaign are presented. The results in terms of load multiplier (i.e. ratio between horizontal and vertical loading) λ α are listed in Table 1, while full details of the experimental dataset are publicly available in [55]. As stated, the repeatability and reliability of the results were ensured by including five tests for each direction (the CoV is very small, demonstrating the low scatter). Additionally, similar collapse mechanisms have been found for tests belonging to the same group. As Fig. 5 illustrates, the collapse load multiplier decreases following a non-linear trend from π /4 to π /2. It can be observed that the highest capacity of the system is obtained when the load direction is coincident with the corner bisector, i.e. π /4, resulting in an increase in the capacity by 21% compared to the π /2 case. Indeed, when rotating the corner towards the π /4 orientation, the out-of-plane component acting on each wall is reduced. Additionally, the different types of failure mechanisms, namely rocking, sliding and horizontal flexural, depend on the orientation of the corner. The collapse mechanism for the loading conditions π /4,5 π /16,3 π / 8,7 π /16 and π /2 are depicted in Fig. 6. The configurations for which the load is applied in the direction of the bisector, i.e. π /4 (Fig. 6a), develop a mixed mechanism that comprises both in-plane and out-ofplane components in each wall. Both walls show a predominant outof-plane component, appearing at first with an out-of-plane sliding of the units. This initial phase is then followed by an uplift of the blocks that contributes to the formation of three cylindrical hinges in both walls: i) the first one is located in the vicinity of the boundary conditions, i.e. line O’A’ and line O’’A’’, and appears fully open, ii) the second one is fully formed and shows up in a central wall area corresponding to the inner wall side, i.e. line O’B’ and line O’’B’’, iii) the third one develops in an area close to the wall intersection and it is only triggered in the left wall, line O’C’, and completely open in the right wall, line O’’C’’. Fig. 6b displays the collapse mechanism of the 5 π /16 case. In this configuration, the left wall shows a similar mechanism to the π /4 case, while the right wall presents a higher displacement in-plane, followed by an out-of-plane mechanism. Similar collapse mechanisms have been registered for 3 π /8 (Fig. 6c), 7 π /16 (Fig. 6d) cases and the weakest condition, namely π /2 (Fig. 6e). The right wall works in its plane, showing a clear initial sliding phase of the units followed by a rotation around O’’ (Fig. 6c, d, e-right). At the same time, the orthogonal wall experiences a predominant out-of-plane mechanism, which results in the creation of two macroblocks, namely O’A’B’ and O’B’C’ (Fig. 6c, d, eleft). It is worth noting that the macroblock O’A’B’ is formed thanks to two cylindrical hinges fully formed, i.e. O’A’ and O’B’, while the onset of the O’C’ is triggered but not fully formed. In particular, the cylindrical hinge O’C’ dissipates a minimal amount of energy, resulting in pseudo orthogonality between the two walls in the kinematic collapse configuration. 3. Numerical (DEM) modelling Numerical analyses of the corner specimens have been conducted using DEM in 3DEC [57]. In this section, the model setup is first briefly presented, following which its predictions are validated against the experimental results. 3.1. Model setup Using DEM, the individual units making up a structure can be modelled as either rigid or deformable blocks. In the latter case, the blocks are discretised using a finite element mesh. In the present work, the specimens are modelled as an assembly of rigid blocks, as these are more suitable for quasi-static analyses [57]. All the system deformability is therefore concentrated in the joints and each block is characterised by the experimental density of 1876kg⋅m−3 and 6 degrees of freedom. Interaction between the blocks is modelled using zero-thickness nonlinear springs (point contacts), which are automatically created when two block faces are determined to be in contact. The point contacts approach allows the modelling of large displacements and complete separation between blocks, which makes it particularly suitable for problems of stability and/or collapse. In the case of rigid blocks (as used herein), point contacts (sub-contacts) are generated through triangulation of the block faces and are usually created at each block face’s vertices [57]. By default, rectangular block faces are divided into two triangles, resulting in four-point contacts (Fig. 7); however, more contact points are usually required to capture the contact stress distribution accurately [21]. In this study, the triangulation density is increased so that each block comprises 10, 7 and 6 contact points along its length, height, and thickness, respectively. Using the contact points approach, the resultant force increment ΔF at each contact point can be determined through the relative incremental displacement Δu between the blocks at that point, the contact spring stiffness k, as well as the subcontact area Ac. The latter is typically equal to 1/3 of the areas of the triangles containing the contact point (Fig. 7). As a soft contact approach is used here, the extent of interpenetration between blocks is controlled by the stiffness of the springs. Though Section 2.1 evidenced the dependency of stiffness on axial stress, in this work, stiffness is considered homogeneous and is taken equal to the average value in the specimens, as given next. From the geometry and density of the specimens, the axial stress at mid-height of the specimen is σ = (0.633×1.876 ×9.81)/2=5.82 kPa. According to Eq. (2), the Table 1 Experimental results, load multiplierλ α . Load orientation α Load Multiplier λ α Mean Value λ CoV [%] π /4 0.287 2.0 5 π /16 0.270 3.9 3 π /8 0.248 3.0 7 π /16 0.239 1.3 π /2 0.236 3.6 Fig. 5. Experimental load collapse multiplier λ α obtained for all the orientations tested. C. Colombo et al. Construction and Building Materials 344 (2022) 127969 6 Fig. 6. Failure mechanisms for (a) Corner π /4, (b) Corner 5 π /16, (c) Corner 3 π /8, (d) Corner 7 π /16 and (e) Corner π /2. C. Colombo et al. Construction and Building Materials 344 (2022) 127969 7 average normal stiffness of the dry-joints therefore equals kn=195.6× σ =1.14 GPa⋅m-1. Similarly, the tangential stiffness is determined to be ks=0.33 GPa⋅m-1. The non-linear response of the point contacts is controlled via the tensile strength ft in the normal direction and the Coulomb-slip joint model in the shear direction (Fig. 8), which consequently depends on the values specified for the cohesion c and friction coefficient μ (indicated by μ =tanφ in Fig. 8). Herein, null values are selected for the tensile strength and the cohesion of the dry-joints, while the dilatancy angle ψ is set to zero, according to the direct shear tests (Section 2.1) and literature [49]. On the other hand, the friction coefficient is set to the value of μ =0.66 (φ=34◦)obtained through the tilting tests. With respect to the boundary conditions, the blocks at the end of the walls are fixed to reproduce the clamped conditions of the experiment. Similarly, the bottom course is fixed in order to model the effect of the double-faced tape (see Fig. 9, shown here for the π /4 corner test setup). In terms of loading, to simulate the tilting table, the vertical component of gravity is kept constant at10 m⋅s-2, while the horizontal component is progressively increased in steps of 0.05 m⋅s-2. For each load increment, the analysis is run until the maximum velocity in the system is less than 0.001 m⋅s-1, and only once the convergence criterion is satisfied is the loading increased. Failure is considered to occur with a lack of convergence of the maximum velocity, and the simulations are continued until the maximum displacement of the system exceeds 0.7 m, to obtain the complete collapse of the structure. The equations of motion of the system are solved by means of an explicit time-stepping scheme, which makes use of a central-difference algorithm [21]. At each timestep, the relative velocity across a subcontact (contact point) is first determined and then converted into a relative incremental displacement Δu. This displacement is resolved into its normal and shear components, Δun and Δus, and is consequently used to determine the force increments. The total normal and shear forces for the sub-contacts are subsequently updated and adjusted based on the contact constitutive relations (in this case, the Coulomb-slip joint model). The sub-contact forces are added to the forces/moments acting on the relevant block centroids, and the corresponding equations of motion are then integrated to obtain the new block velocities and positions, which serve as input for the next calculation cycle. Note that to ensure numerical stability, a suitably small timestep is required for the analyses, which is automatically calculated by the program and is directly proportional to the smallest nodal mass and inversely proportional to the maximum contact stiffness of the system. Consequently, a decrease in nodal mass tends to increase the overall solve time (e.g. more contact points, resulting from an increase in the triangulation density). For static and pseudo-static solutions, the equations of motion are damped in such a way that they reach a force equilibrium state as rapidly as possible. Conceptually similar to dynamic relaxation, a damping force, proportional to the velocity of the blocks, is artificially added into the equations of motion to avoid oscillations around the equilibrium position [57]. Here, local damping with the default value of 0.8 is used, whereby an additional/artificial damping force is applied on each node. This force is proportional to the magnitude of the unbalanced force at that same node, in a direction such that energy is always dissipated [57]. This form of damping is usually recommended for Fig. 7. Mechanical representation of contact between blocks. Fig. 8. Coulomb friction joint (a) normal stress versus normal displacement and (b) shear stress versus shear displacement behaviour (adapted from [57]). Fig. 9. Final 3DEC model. C. Colombo et al. Construction and Building Materials 344 (2022) 127969 8 static/pseudo-static analyses and can minimise oscillations arising from abrupt failure of the model [22,57]. 3.2. Experimental validation The results of the numerical models have been evaluated in terms of the load multipliers λ α and associated collapse mechanisms obtained for each load orientation angle α . In terms of the collapse load multipliers, as Table 2 and Fig. 10 indicate, the predictions of the numerical model generally tend to follow the same trend as the experiments, with the load multiplier globally decreasing in a non-linear manner from π /4 to π /2. However, the numerical model appears to consistently overestimate λ α , with an average overprediction of about 10%. This slight overestimation is likely to be partly due to the fixed boundary conditions assumed in the numerical model, which do not permit any rotations and consequently increase the stiffness of the walls in the out-of-plane (OOP) direction. The effect of this assumption increases as the OOP load component acting on the individual walls grows and explains why, in the case of the π /2 orientation, where one of the walls is loaded entirely in the OOP direction, a higher load multiplier is required to cause a local collapse. In addition to the effect of the boundary conditions mentioned above, the overestimation of the load multiplier is also attributed to other assumptions of the numerical model, such as perfect (zero-thickness) joints between blocks and uniform cross-sections within the blocks. This first assumption is most likely not the case in reality for the vertical (head) joints, where small gaps between the blocks could decrease the in-plane and flexural capacities of the walls. Additionally, given that the CS blocks used in the experiments have been manually sawn, it is unlikely that they would have perfectly uniform crosssections. Given that the walls are already quite slender in the out-of-plane (OOP) direction (H/t=10), small differences in the thickness of the blocks (decrease of a few mm) may increase the overall slenderness and therefore decrease the capacity of the wall in the out-of-plane direction as well. Thus, two sets of additional numerical simulations have been run with a reduced thickness (increased slenderness). The first uses a 5% reduction in thickness (and therefore an approximately 5% increase in slenderness), while the second uses a 10% reduction in thickness (or conversely a 10% increase in slenderness). Note that a thickness reduction higher than 10% (i.e. 5.7 mm) was not considered reasonable given the near-perfect edges of the CS blocks. The results of these analyses can be found in Fig. 11a and Table 3, where it can be seen that although an increase in slenderness does cause a slight decrease in the collapse load multiplier, the numerical simulations still overestimate λ α , with an average overprediction of 8% and 6% for the 5% and 10% thickness reduction cases respectively. Thus, the non-uniform cross-sections of the CS blocks are unlikely to be a significant cause of the discrepancy between the numerical and experimental results. Table 2 Comparison of the experimental and numerical load multipliers λ α , with Difference (%) = (λNUM α −λEXP α )/λEXP α ×100. Load orientation α Load Multiplier λ α Difference (%) Experimental Numerical π /4 0.287 0.315 9.8 5 π /16 0.270 0.290 7.4 3 π /8 0.248 0.270 8.9 7 π /16 0.239 0.265 10.9 π /2 0.236 0.270 14.4 Fig. 10. Variation of collapse load multiplier with load orientation for both the experimental tests and numerical models. Fig. 11. Variation of collapse load multiplier with load orientation for both the experimental tests and numerical models, for (a) different slenderness (H/t)ratios and (b) different stiffnesses (k n,ks). C. Colombo et al. Construction and Building Materials 344 (2022) 127969 9 In fact, another source of discrepancy could be attributed to the values adopted for the joint normal stiffness kn, which was herein evaluated through classical joint closure tests. In [44], the authors have demonstrated that the local normal stiffness of joints belonging to a staggered masonry wall is approximately 2 to 3 times lower than that evaluated through joint closure tests. This is due to the existence of misalignment and dimension differences between blocks in a wall, where contact is never only between two blocks (case of the joint closure tests). To that end, two additional sets of simulations have been conducted – the first one reducing the normal stiffness by a factor of 2 (or 50%) and the second reducing the normal stiffness by a factor of 3 (so to approximately 30% of the original stiffness). In the absence of similar information about the tangential stiffness ks, it was chosen to keep the initial ratio of ks/kn constant for these analyses. The results of these additional analyses can be found in Fig. 11b and Table 4. As Fig. 11b illustrates, the numerical models demonstrate a sensitivity to the joint stiffness. Decreasing this parameter by a factor of 2 causes a substantial decrease in the collapse load multiplier, with the numerical predictions now lying within the range of the experimental scatter, with an average overprediction of 4% (with respect to the mean experimental values). Similarly, reducing the joint stiffness by a factor of 3 further reduces the collapse load multiplier, with the numerical predictions now providing a lower bound to the experimental results and, on average, underpredicting λ α by 3%. These analyses thus demonstrate the importance of evaluating the joint stiffness on the whole masonry wall, and not just through joint closure tests. Nevertheless, it is worth underlining that the sensitivity analyses investigate the variation of a single parameter individually (i.e. slenderness or stiffness), while most likely the global outcome includes the influence of combined factors. For this reason, it is understood that the adopted range of stiffness variation may be reduced if concurrently considering other sources of uncertainties. With respect to the collapse mechanisms, a near-perfect agreement has been found between the numerical models and the experiments in terms of the blocks involved in each mechanism (see Fig. 12). Note that as similar mechanisms were obtained from the four sets of parametric analyses, only the mechanisms from the original numerical analyses are presented here. A good agreement was also observed in terms of the shape of the failure mechanism and the associated hinge locations, as illustrated by Fig. 13, with the formation of horizontal flexural mechanisms in the out-of-plane loaded walls and rocking-sliding mechanisms in the in-plane walls. The two walls have also been observed to retain their orthogonality relative to each other. The main difference observed is in terms of the location of hinge B’ on the OOP-loaded wall, which occurred one block closer to the corner in the numerical models, resulting in a slightly larger macroblock O’A’B’ than the experimental tests for all considered load orientations. 4. Analytical model The failure mechanism of the corner through an analytical model has been recently investigated by [15]. A three-dimensional macroblock model with frictional joints to analyse failure mechanisms was developed, with crack patterns and load factors derived through the kinematic approach of Limit Analysis (LA). This formulation considers either the rocking–sliding mechanism or the horizontal flexure phenomena Table 3 Comparison of experimental and numerical load multipliers λ α for different slenderness (H/t)ratios, with Difference (%) = (λNUM α −λEXP α )/λEXP α ×100. Load Multiplier λ α Load orientation α Experimental Numerical (H/t)⋅1.05 Difference (%) Numerical (H/t)⋅1.10 Difference (%) π /4 0.287 0.310 8.0 0.305 6.3 5 π /16 0.270 0.285 5.6 0.280 3.7 3 π /8 0.248 0.265 6.9 0.260 4.8 7 π /16 0.239 0.260 8.8 0.255 6.7 π /2 0.236 0.265 12.3 0.260 10.2 Table 4 Comparison of experimental and numerical load multipliers λ α for different stiffnesses kn and ks, with Difference (%) = (λNUM α −λEXP α )/λEXP α ×100. Load Multiplier λ α Load orientation α Experimental Numerical (k n , k s )⋅0.5 Difference (%) Numerical (k n , k s )⋅0.3 Difference (%) π /4 0.287 0.300 4.5 0.280 −2.4 5 π /16 0.270 0.270 0.0 0.255 −5.6 3 π /8 0.248 0.255 2.8 0.240 −3.2 7 π /16 0.239 0.250 4.6 0.235 −1.7 π /2 0.236 0.255 8.1 0.240 1.7 Fig. 12. Comparison in terms of blocks involved in the collapse mechanisms for the experimental and numerical tests, for the different load orientations. C. Colombo et al.