Numerical evaluation of test setups for determining the shear strength of masonry

The bond shear strength between masonry units and mortar is the weakest link in a masonry wall. Different material tests have been developed in order to characterize this bond behaviour. The objective of this study is to evaluate three common test setups through non-linear finite element analysis. The simulation method is based on our recent development of cohesive elements, which allows for the first time to fully capture the force-deformation characteristic of shear tests in 3D from the onset of loading until the residual shear strength and to retrieve typical shear failure modes observed in experiments. This study provides new insights into our understanding and interpretation of such shear tests: (1) elastic analysis, which has been widely used in the past, does not yield a stress distribution that is representative of the stress distribution at maximum resistance; (2) while friction coefficient is well estimated (the error is less than 10%), the local cohesion is underestimated by all three test setups of which the error lies between 13 and 32%; (3) the randomness of the material properties leads to a further underestimation of the mean value of the local cohesion; (4) differences in the material properties of the two joints of the triplet test units do not jeopardize the applicability of this test setup and estimations of the mean properties are obtained with similar reliability as for couplet tests.


Introduction
Under seismic loading, unreinforced masonry (URM) walls can fail due to in-plane or out-of-plane loading. For in-plane loading, walls can fail in shear, flexure or a combination of the two failure modes [1]. While the flexural strength is relatively independent of material properties and governed by static and kinematic boundary conditions as well as the geometry of the element, the shear strength is heavily dependent on material properties. The first shear cracks pass typically through the joints and form stair-stepped cracks (e.g. [2,3]). To predict the shear behavior of URM walls by means of analytical or numerical models, information on the joint shear strength including the cohesion and friction coefficient are required [4]. Different test procedures have been proposed by various researchers, indicating the difficulty of finding a consensus with regard to the best testing method. Interested readers can consult [5,6] for a review of the historical development of the various testing methods. This paper focuses on three frequently used test setups: The original version of this article was revised due to a retrospective Open Access order.
(1) A test setup applied by Van der Pluijm [7] (Fig. 1a). The results of this test series have been widely used for the calibration of numerical models [8]. (2) A test method proposed by Lourenço and his collaborators [9] (couplet test), which is similar to the classical shear box used in geomechanics (Fig. 1b). (3) The triplet test (Fig. 1c), which has been adopted by the European Committee for Normalization as the standard test for determining the joint shear properties [10]. For different experimental setups, cohesion is known to be the most sensitive parameter, while the friction coefficient is less sensitive to the setup [4].
In order to evaluate the performance of the various setups, previous studies compared the normal and shear stress distributions at the center line of the mortar joint that are obtained from elastic analysis. The following quality criteria were developed by Riddington [11] for the evaluation of the test setups and have been used since by others [5,12].
1. The shear and normal stress should be uniform along the joint length. 2. When failure is initiated at one point, the other parts of the joint should be close to failure too. 3. Tensile stresses in the joint should be avoided. 4. The failure should not be initiated at the edge. 5. The experiment should be easy to implement.
The objective of these criteria is to ensure that the material properties that are back-calculated from the test match as well as possible the local ones. To back calculate the cohesion and coefficient of friction from the experimental results, the maximum shear resistance needs to be determined at different normal stress levels. In the ideal case where criteria 1-4 are satisfied, the stress state at damage initiation is close to the stress state at maximum shear resistance and elastic analysis would lead to a representative shear stress distribution at maximum shear resistance. The estimated material properties obtained would be close to the local ones. However, as will be shown later, even the more complex setup by Van der Pluijm fails to fully meet the criteria 1, 2 and 3.
Early research on these setups, which was conducted in the 90s, concentrated on elastic analysis until damage initiation [6,11]. To better simulate the failure process and estimate more accurately the stress distribution for a small scale test like a shear test, a full representation of the geometry, i.e., representing brick, mortar and interfaces in the model, is necessary. This is known as detailed micro-modeling approach [5,8,[13][14][15]. Previous studies that used such modeling approaches have not been able to reproduce the full force-displacement relation for simple shear tests (e.g., [5,13,16]). This paper has five objectives: (1) evaluating to which extent conclusions on the performance of the test setups that were drawn from elastic analysis are representative when evaluating the maximum shear

Numerical formulation
The finite element analysis is conducted with the open source library Akantu [17]. In Sect. 2.1, the numerical framework is introduced. To fully capture the shear debonding process, a bilinear descending law is proposed in Sect. 2.2.

Explicit integration and extrinsic insertion
After spatial discretization of the weak form of equilibrium equation, the following well-known relationship is obtained: in which M is the mass matrix, € u is the acceleration vector, f int and f ext are the internal and external force vectors. The classical explicit second order central difference method is used here for time integration. The displacement, velocity and acceleration (u mþ1 , _ u mþ1 , € u mþ1 ) at time step m þ 1 are estimated by A constant time step Dt is used during simulation, which is confined by in which c l represents the longitudinal wave speed, l mine is the characteristic length of the minimum element, a is a safety factor, chosen to be 0.15 here.
Calculating f int requires a constitutive law for the bulk elements and a traction-separation law for cohesive element. For bulk elements, an isotropic linear elastic relation is assumed. Material nonlinearity comes from the cohesive elements [18,19].
Here the extrinsic approach [20] is used, for which cohesive elements are inserted dynamically during the simulation while the following criteria is met in which r c is the critical stress, r eff is the effective stress for the current state calculated by [20] for tension/shear t n ! 0 1 b Á ðjt s j À ljt n jÞ for compression/shear t n \0 in which t n ¼ r Á n, t s ¼ r Á s are the tractions in the normal n and tangential s directions of the facets, as illustrated in Fig. 2, b is the shear stress factor, l is the friction coefficient. Two situations, tension/shear t n ! 0 and compression/shear t n \0, are distinguished. For instance, if during the simulation inequality 6 is τ n i j satisfied on red facet, cohesive element will then be inserted. The nodes will then be splitted, e.g., nodes i, j were one node before insertion.
After insertion, the traction is determined by the following traction-separation law [21] T in which j ¼ G c;II =G c;I , b indicates the ratio between cohesion and tensile strength, T is a scalar value determined by the traction-separation law. In [20, 21], a linear descending law is used. While for quasi-brittle material, a bilinear descending law is preferable [22].

Bilinear descending law
The implemented bilinear descending law is illustrated in Fig. 3, which can be represented by the following equation where d c represents the effective separation upon which the cohesive element is totally damaged, d max is the current maximum effective separation, T max is the traction-separation at d ¼ d max , d h and r h indicate the position of the kink point which is determined by h and The effective separation is calculated in the same way as in [21] d ¼ For contact, the same node-to-node contact and friction algorithm as in [15] with an extension to 3D, in which the contact and friction forces, f m cont;ij and f m fric;ij , at time step m for node pair i, j is calculated by predicting the displacement and velocity at time step m þ 1.

Results and discussion
This section starts with a classical elastic analysis of Van der Pluijm's test setup (Sect. 3.1). In Sect. 3.2, a non-linear model is set up following the approach outlined in the previous section and the material properties are calibrated by fitting experimental results from Van der Pluijm's test [7]. A 3D analysis is further carried out to validate the 2D simulations. In Sect. 3.3, the influence of a random spatial variation of the material properties along the interface is studied. In Sect. 3.4, the simulation results are compared for different test setups with regard to their ability of estimating local strength parameters. This section ends by a discussion on the triplet test (Sect. 3.5). To facilitate comparison, all the test setups use the same brick dimensions (200 mm Â 100 mm Â 50 mm) and mortar thickness (15 mm) as in Van der Pluijm's test [7]. To further consider the fact that the mortar does not occupy the whole space between bricks, a setback distance of half the mortar thickness (7.5 mm) is left for all mortar layers ( Fig. 1) in our model.

Elastic analysis
As mentioned in the introduction, elastic analysis is only valid until damage initiation. In addition, previous studies evaluated the stress state at the center line of the mortar joint [6,11]. However, in shear tests, cracks develop commonly along the interface between unit and mortar, which is weaker than the mortar itself [23]. Due to the finite thickness of the mortar joint, the stresses along the interface can differ significantly from the stresses at the center line of the mortar joint, as will be shown in this section. This holds also for the Van der Pluijm test (Fig. 1a), which was especially designed to minimize the effect of the bending moment on the stress distribution within the mortar. The adopted elastic material properties are included in Table 1. As an example, for a normal stress of 0.5 MPa, the stress distribution is computed at the damage initiation state, i.e., at the instant when the criteria of Eq. 6 is first met at the boundary between unit and mortar and therefore the first cohesive element is inserted. For this state, the normal stress and shear stress distributions are shown in Fig. 4a. The distance from the insertion criteria, calculated by r c À r eff (in this subsection, to concentrate on the difference caused only by position, the same insertion criteria r c for mortar and for interface is temporarily assumed), is shown in Fig. 4b. As can be seen from this figure, the first cohesive element is inserted in the middle of the specimen. Note that the damage will normally be initiated only at the interface, because the mortar is generally much stronger than the interface, i.e., mortar has a higher r c as compared to the interface.
The test setup has been designed in such a way that the moment at the center line of the mortar joint is zero. This is achieved by applying the horizontal force F at that height. Since the joint has a certain thickness (15 mm here), the moment at the two interfaces between unit and mortar will not be zero. As a result, the normal stress distribution is no longer symmetric ( Fig. 4a) with regard to the y-axis (Fig. 1a). Due to the normal stresses that result from the bending moments at the heights of the interfaces, the left part of the interface 1 and the right part of the interface 2 have a lower shear resistance, as shown in Fig. 4b. This suggests that under uniform material properties, the cracks will propagate from the middle to the left along interface 1 and from the middle to the right along interface 2 at the same time.
This section shows that the additional moment that results from the thickness of the mortar joint is non  negligible. The stress distributions should be evaluated at the interfaces, instead of at the center line of the mortar joint, when accessing a shear test setup through criteria 1 to 4.

Calibration material properties with Van der Pluijm's test
To retrieve force-deformation characteristics for different compression levels, the fracture energy is assumed to be determined by the global normal stress using the following linear relation in which r N is the normal stress, a and b are constants. The parameters a, b and other material parameters are chosen by fitting the experimental curve in [7]. The calibrated material parameters are listed in Table 1.
More specifically: the brick elastic modulus and density are already given in [7], while the elastic modulus of mortar is selected to match the stiffness obtained from the experiment; a typical value of 0.15 is selected for the Poisson's ratio of brick and mortar; the inelastic parameters of cohesive elements are selected such that force-displacement characteristic matches experimental results, with the assumption that the mortar is three times stronger than the interface. Experimental information on mortar to interface strength ratio does not exist, but a parametric study can be found in [15]. A 2D analysis is firstly carried out. In the mesh, 1 0 564 second order plane strain elements are used, within which 444 9-node quadrilateral elements form a structured mesh for the mortar and brick, while 1 0 120 6-node second order triangular elements form an unstructured mesh for the loading shoe. The comparison with experimental data is presented in Fig. 5, which shows that the simulation results and experimental data match well. Further refining the mesh has limited influence on the maximum shear resistance, e.g., under compression level of 0.5 MPa, decreasing the element size by half only results in 2 % difference for the maximum shear resistance. Unless indicated otherwise, the mesh resolution for other simulations in this paper are as outlined at the beginning of this paragraph.
To further validate our 2D simulation, a full 3D analysis is carried out for the compression level of 0.5 MPa. In the 3D mesh, 70 0 344 10-node tetrahedron elements are used. The material parameters are the same as in Table 1. Figure 6 compares the forcedisplacement relation and the inserted cohesive elements at maximum shear resistance. The maximum shear resistance obtained by the 2D analysis is only 2.4% lower than the maximum shear resistance obtained with the 3D simulation. The crack patterns at maximum shear resistance are also similar. It can therefore be concluded that 3D effects are negligible and that the performance evaluation of the test setups can be based on the results of 2D simulations. In the previous calibration, the interface properties are assumed to be uniform along the interface. However, in reality, the material properties will vary along the interface due to the natural aleatoric variability of material properties. The variation of the interface properties will affect the initiation and propagation of the crack, and further more will influence the obtained maximum shear resistance. The objective of this section is to investigate the influence of random variation of material properties on the values estimated for cohesion and friction coefficient on the interface. The cohesion and the friction coefficient, are calculated by evaluating the maximum shear resistances at different normal stress levels [10]. Here five normal stress levels r N;i are selected and the maximum shear resistances s i is obtained from simulation. A linear regression is then conducted, with regard to the estimated cohesion and friction, c est and l est (Fig. 7).
in which s i is the maximum shear resistance corresponding to r N;i , calculated by the following equation where F i;max is the maximum shear load under normal stress r N;i , A is the nominal cross-sectional area of a specimen parallel to the bed joints, e.g., for the current specimen A ¼ 200 mm Â 100 mm. With uniform material properties on the interface, the maximum shear resistances under five normal stress levels are indicated by ''w/ mean properties'' in Fig. 7. The corresponding linear regression of estimation is noted as ''estimation/mean''. From the linear regression, we have c est ¼ 0:95 MPa and l est ¼ 0:83.
For considering randomness, the assumption is adopted that the variation of interface properties follows a Gaussian random distribution. Spectral representation is used for generating random samples [24]. An alternative approach is stochastic harmonic function representation method [25]. The critical stress, fracture energy, and friction coefficient of the cohesive elements are considered to be random variables and are assumed to be fully correlated. The standard deviation for critical stress and the fracture energy is assumed to be 0.3, while for the friction coefficient, the standard deviation is assumed to be 0.2. The correlation length is assumed to be 0.0125m, 1/16 brick length.
Under each normal stress level, eight samples are generated and the maximum shear resistances are indicated by ''w/random properties'' in Fig. 7. The corresponding linear regression of estimation is noted as ''estimation/mean'' from which it is obtained that c est ran ¼ 0:85 MPa and l est ran ¼ 0:73. Therefore, failing to consider randomness will lead to a 17% overestimation of cohesion.
We then compare the estimated cohesion with the local cohesion. Due to the set backs of the mortar joint, the effective interface area between mortar and brick is smaller than the nominal cross section of the brick. The mean value of the local cohesion is calculated by in which A is the nominal cross-sectional area defined in Eq. 15, A net is the net cross-sectional area of a specimen parallel to the bed joints, e.g., for the current case A net ¼ 185 mm Â 100 mm, b and r c are from Table 1. Therefore we have c ¼ 1.11MPa. Since variation of material properties is unavoidable in reality, a larger correction factor, compared to the assumption with uniform material properties, should be used when calculating local cohesion value from estimated cohesion. For example, for Van der Pluijm's test, which was analyzed here, it is obtained that c ¼ Max. shear resistance (MPa) w/ mean properties estimation / mean w/ random properties estimation / random 1:31c est with the specified random field and c ¼ 1:17c est with uniform material properties. When considering random properties along the interface, another typical failure mode, where the failure plane switches from one interface to the other, can also be retrieved (failure mode shown in the lower right corner of Fig. 7). This failure mode is also often observed in experiments [7].

Compare the ability of different test setups to estimate local material properties
In this subsection, the stress distributions, estimated cohesion and friction coefficient for the three test setups in Fig. 1, are compared. To simplify the discussion, the material properties is assumed to be the same as in Table 1 and random variation is neglected. The cohesion and the friction coefficient are estimated by Eq. 14 while the maximum shear resistance s i is calculated by the following equation where the parameters are already defined in Eq. 15. For comparison, the Mohr-Coulomb law in Fig. 8 is also plotted, which is obtained from the material properties specified in the analysis: in which r N is the average normal stress, l is the friction coefficient, and c is the local value of the cohesion, defined by Eq. 16. Figure 8 shows that: (1) for all three test setups, the local cohesion is underestimated. More specifically, the underestimation for the Van der Pluijm setup, the triplet test setup, and the test setup by Lourenço is 14, 13, 32%, respectively; (2) the error of the friction coefficient is 3, 4, 10% for the three test setups respectively; (3) Fig. 8 also shows that the error related to estimating the cohesion and friction coefficient is also influenced by the range of the normal stress considered when conducting a linear regression (Eq. 14), which represents the Mohr-Coulomb law.
The maximum shear resistances obtained from the Van der Pluijm test setup and the triplet test setup are similar. However, the stress distributions are in fact quite different. For a normal stress level of 0.5MPa, Fig. 9a shows the distance from the insertion criteria at damage initiation and also at maximum shear resistance (indicated by thick lines, for which ''From insertion criteria'' is only calculated at positions where cohesive elements are not yet inserted ). For the triplet test and for the test setup by Lourenço, the damage initiates from the edge of the interface at a very early stage, 35-40% of maximum shear resistance. All other regions are still far from insertion at this moment. For the Van der Pluijm test setup, damage initiates quite late, almost at 90% of maximum shear resistance.
However, at maximum shear resistance, a large part of the interfaces has been damaged in all three setups; the damaged parts are indicated by the thin line in Fig. 9 and by the hatched region in Fig. 9a. The normal and shear stresses along the interface corresponding to the maximum strength for different test setups are shown in Fig. 9b, in which thick lines indicate again the positions where cohesive elements have not yet been inserted. As shown in Fig. 9b, the criteria 3, which says that tensile stresses should be avoided, is violated even for the Van der Pluijm test setup, and the advantage of Van der Pluijm's test setup over the other test setups is less obvious at maximum shear resistance than at the point of damage initiation (Fig. 9). Thus, even if the stresses are evaluated at the interface level, elastic analysis does not yield stress distributions that are representative of the stress distribution at maximum shear resistance and non-linear analysis is required to fully evaluate the validity of a test setup. At maximum strength, for the test setup proposed by Lourenço, cohesive elements have not yet been inserted over a large part of the interface, which explains why the estimated cohesion is smaller compared to other two test setups.

Discussion on triplet test
Despite serving as a standard test method and been widely used [13,26,27], the interpretation of triplet tests is difficult because the two joints do not fail at the same time [28]. In reality, the mortar properties of the two joints are not the same. An absolutely symmetrical force boundary condition, as shown in Fig. 1c, is also often difficult to realize due to constraints in the test setup. Sometimes a displacement boundary condition is used on one side and a force boundary condition on the other [27]. The effects of the different material properties of the two joints and different boundary conditions with regard to the x-axis in Fig. 1c have not been well understood. To examine the validity of the triplet test and to deepen our understanding, the influence of three factors on the maximum shear resistance obtained are studied here, i.e., the influence of (1) the boundary conditions, (2) a difference in the elastic modulus of mortar of the two joints and (3) a difference in the interface properties of the two joints. To simplify the discussion, only the setup for the intermediate compression level with a normal stress of 0.5MPa is analyzed.
In Fig. 10, two different boundary conditions are considered: The term ''force boundary'' refers to two Neumann boundary conditions, fixed to 0.5MPa, applied on either side of the triplet (Fig. 1c). For the ''displacement boundary'' case, which was used in the previous sections, a force boundary is applied only on one side (right side) of the specimen while the y displacement is fixed on the left side. Hence the rotation is restraint. With constant and equal material properties along both joints (Table 1), the maximum strength obtained for the force boundary condition is 5% lower than the one obtained for the displacement boundary condition.
The influence of differences in the interface properties of the two joints is studied next. The strength, i.e., critical stress, fracture energy and friction coefficient, is decreased for the two interfaces of the right joint. The average of the two interface strengths, indicated by ''predicted strength'' in Fig. 10a, is expected to be in which sðaÞ is the mean maximum shear resistance, and f intf1 and f intf2 are the two weakest interface (a) (b) Fig. 9 Comparison of stress distributions for different test setups at damage initiation and at maximum shear resistance (thick lines indicate undamaged parts). a Distance from the insertion criteria. b Stresses at maximum shear resistance maximum resistance on each side. In Eq. 19, it is assumed that a indicates the reduction of the interface strength on the right side, and f intf is the maximum strength obtained with constant and equal material properties, e.g., under normal stress 0.5MPa, for ''displacement boundary'' f intf ¼ 1.41MPa (Fig. 8), while for ''force boundary'' f intf ¼ 1.34MPa. The maximum shear resistance is normalized by sðaÞ=f intf . The maximum resistance obtained from simulation is shown in Fig. 10a. While interface strength of the right joint is reduced to 60% of the initial value (a is equal to 0.6), the estimated strength reduces from f intf to sð0:6Þ ¼ 0:8f intf (Eq. 19), i.e., from 1.41 to 1.13 MPa for ''displacement boundary'' and from 1.34 to 1.07 MPa for ''force boundary''. For the displacement boundary condition, the strength is generally overestimated, with the maximum error 7% reached at a ¼ 0:1 (Fig. 10a). While for the force boundary condition, the strength is generally underestimated. The error increases with the decrease of a and reaches the maximum 18% for a ¼ 0:6 ( Fig. 10a). With small difference between the two mortar joints, i.e., 0:8\a\1:0, ''force boundary'' exhibits a higher accuracy. While with large difference between the two mortar joints, i.e., 0:6\a\0:7, ''displacement boundary'' is preferable (Fig. 10a). It is expected that in real setups the boundary conditions will often fall in between the displacement and force boundary condition; for this reason, the error of estimation is likely to be acceptable.
Another factor that may influence the result is a different elastic modulus of the mortar in the two joints. To investigate this factor, the mortar elastic modulus on the right side is reduced by up to 40% of the initial value. The maximum strengths obtained for the two boundary conditions are plotted in Fig. 10b. It can be seen that its influence on the maximum shear resistance is negligible (error within 2%) for both boundary conditions.
The typical failure modes for the two boundary conditions (with reduced mortar elastic modulus for the right joint) are shown in Fig. 10b. For force boundary conditions, the failure surface is along the left mortar layer where the mortar layer is stiffer (the failure mode on the right). For the displacement boundary condition (the failure mode on the left), the failure mode is different because the rotation is restrained. Despite different failure modes, it is interesting to notice that the two specimens have almost the same maximum shear resistance (Fig. 10b).

Conclusion
The purpose of the current study is to re-evaluate three common setups for shear tests. Based on our recent development on cohesive elements, for the first time, the force-deformation characteristic can be fully captured in 3D (up until the residual shear resistance) and typical failure modes are retrieved for shear tests on couplet or triplet samples. This study has shown that: (1) when assessing different shear test setups, due to the non-negligible bending moment caused by the mortar joint thickness, the stresses should be directly evaluated at the interfaces between mortar and unit instead of at the center line of the mortar joint; (2) elastic analysis, commonly conducted in previous studies, does not provide a realistic stress distribution at maximum resistance; it is only representative up to damage initiation, which starts at a force around 35-40% of the maximum resistance for the triplet and Lourenço's test setup; for Van der Pluijm's test, although the damage initiates at 90% of the maximum resistance, the stress distribution at damage initiation is still significantly different from the stress distribution at maximum resistance. (3) contrary to the common belief, initiation of damage at the extremity of the interface does not significantly influence the maximum shear resistance; (4) for the shear tests analyzed here, 3D effects are negligible and the performance evaluation of the test setups can be based on the results of 2D simulations; (5) the random variation of the material properties has a non-negligible effect on the estimation of the cohesion, e.g., for Van der Pluijm's test setup, it was found, for the material distribution assumed here, a 17% difference in the estimated mean values of the cohesion using constant material properties and using random properties; (6) while the accuracy of estimation for friction coefficient is rather good (error less than 10%), the local cohesion is underestimated by all three test setups. The error obtained for the case studies analyzed here were between 13 and 32% with the smallest error obtained for Van der Pluijm's test setup and a similar error for the triplet test; (7) the maximum shear resistance obtained from triplet tests is influenced by the boundary condition and difference of the interface properties in the two mortar joints; however, the error introduced by these effects is limited. This paper justifies the use of triplet tests for determining the cohesion and friction coefficient of mortar-unit interfaces. However, as with other test setups that have been proposed for this purpose, the estimated cohesion from the triplet test is lower compared to the local cohesion. Therefore, a correction needs to be introduced, if the estimated cohesion is to be used in a detailed micro-modeling approach. The correction factor of the estimated cohesion for the case studies analyzed here is around 1.15, with constant material properties. Note that the correction factor depends on parameters such as specimen dimensions, material parameters, and material randomness. Further studies are required to generalize such factors. In addition, the proposed traction-separation law (Sect. 2.2) is based on the classical one proposed by Camacho and Ortiz [20], which is only valid for monotonic loading conditions. Dilatancy effect is not considered by the current tractionseparation law, related information can be found in [29][30][31]. Another limitation of the traction-separation law is the fixed dependency on normal stress which can affect the simulation accuracy. Further work needs to include a dynamic dependence of the tractionseparation law on the normal stress.