An orthotropic damage model for masonry walls with consistent damage evolution laws

This research presents a constitutive model for the macro scale simulation of masonry structures. The model is containing an orthotropic plane stress assumption, which appears as an appropriate assumption for the in-plane analysis of masonry walls. The material model is based upon damage mechanics, split into tensile and compression parts. The novelty herewith is the consistent mapping of the damage evolution laws. Aim of this research is to develop a simple but accurate constitutive law, suitable to simulate large scale structures and buildings with reasonable computational costs. The developments are presented and compared with available literature examples from laboratory testings.


Introduction
Ever since the use of masonry in buildings and monuments, it has been a challenging task to accurately predict the behaviour and durability of the construction material.With the introduction of finite element methods, numerical methods have eased the assessment and simulation significantly [1].Simultaneously, different scales of modeling have been raised, aiming to analyze different details of the structure.Detailed micro models, which are containing all constituents separately have been developed to show the heterogeneous composite behaviour.Macro models on the other end shall help estimating the structural behaviour of the material.
Macro models overlook the combined kinematics of bricks, mortar and their interface to aim a computationally efficient modelling of the large scale structures.The constitutive behaviour might be obtained from reference volume elements (RVE), where the properties can be gathered experimentally e.g.[2,3], from national codes as the EC 6 [4], or on the local dimensional fidelity requirements.Even though solid and shell elements can better represent the system geometry of masonry structures, dimensionally reduced 1D and spring elements are extensively used for seismic and blast simulations because of their optimized computational costs ( [11][12][13][14][15] for some seismic masonry examples).Continuum-based masonry macro models are typically accurate in the structural simulation (examples for in-plane loading include references [8,16] and out-of-plane loading [17][18][19][20]), however, estimating non-homogeneous material properties as unstructured bonds and detailed crack predictions is, up-to-date, complex.Apart from analytical and continuum based approaches, mesh free or particle based methods find increasing application in masonry classification as well (e.g.discrete element modeling [21], or latice discrete particle model [22]), which facilitates the latter problem of difficulties in crack propagation simulation, however, the computational costs can be tremendous.
This research aims to facilitate simulations capable to express complex geometries and structures with justifiable computational costs for full scale buildings.Auxiliary elements and beams typically have increasingly challenges when modelling complex structures.Whereas, solid and particle based modeling allows almost infinite possibilities in the geometrical specification of constructions, typically the computational efforts prevent the application in structural analyses processes.Thusly, the larger scope shall be on the description of surface-based shell structures, which allow the expression of many geometrical shapes, including reasonable computation costs.
Continuum based damage models define a common way to describe the structural properties of masonry buildings.Initially, one damage parameter was applied to describe the softening of the construction material [23].Where later, those models have been split up into compression and tension damages to allow deeper insights in highly loading dependent materials [24][25][26].Those concepts have been enhanced further with damage parameters per direction [27].Other concepts do not make the damage split within the stress state but the direction [28] with interdependent damage parameters, which is a powerful approach but less appropriate for the apparent masonry properties.Following, within this publication the focus shall be kept on two parameter d þ /d Àdamage models for simplicity with less parameters and fewer internal dependencies.
The mentioned material models are adequate in the analysis of isotropic materials with sensible triggers for the softening zones.Thusly, in masonry, its applications are so-far mostly for the simulation within continuum-based micro models of bricks (without holes) and mortar.However, as the structural properties of masonry show significant correlating orthotropic behaviour it is insufficient to fulfill the apparent direction-dependent estates.Rather than updating the yield surface with anisotropic surfaces as the Hill-type criterion [29][30][31] or Tsai-Wu [32] it is preferred to base upon a mapping from the real anisotropic space to an ideal mapped isotropic stress state.This concept is generally simpler with the possibility to keep former developments and mostly leaner computational costs.The concept of the implicit orthotropic yield criterion was initiated by Betten [33] and later advanced by Oller et al. [34,35].Pela `et al. [16,36] promoted those models to the simulation of masonry structures.The most significant improvement from the macro model of Lourenc ¸o et al. [31] has been achieved by the split in compressive and tensile damages, which allowed to model an exceedingly vast amount of problems.
Within this research the focus is on damage-based models, while plastic zones would be expressed by the damage parameter.This is a fast and generally reliable approach once a one directional loading is applied.For seismic and changing loading cases, plasticity effects should be considered separately for a better accuracy (see e.g.[23,27,37]).
Subsequently, a model which merges the latest improvements of the damage evolution law from Petracca et al. [14,38] and Pela `et al. [16] into a novel material model is presented.The presented research proposes a consistent mapping of the corresponding damage evolution laws with more sensitive triggers in the softening within an appropriate orthotropic behaviour.
This publication is structured as following: • Section 2 presents the applied continuum based d þ /d À -damage material model.• Section 3 shows the damage model with the orthotropic to a mapped isotropic state.• Section 4 gives some explanatory examples with the proposed constitutive model.
• Section 5 summarizes the applied approaches and gives an outlook to future studies.

Continuum tension/compression damage model
Within this section the continuum based d þ /d Àdamage model from [10,17,38] shall be discussed.The aim of this research is the modelling of thin walled masonry.Therefore, the following formulations are expressed in a plane-stress state.

Constitutive relations
The effective stresses r (with the notation Á used to mark the effective state) are obtained from the elastic constitutive tensor C elastic and the actual strains : r ¼ C elastic : ð2:1Þ Similarly to [24][25][26], the effective stresses are split into the positive r þ -tensile and negative r À -compressive constituents, while implying: The two parts can be computed as following: where r i is the i-th principal stress, with the corresponding unitary eigenvector p i .The Macaulay brackets hÁi selects the positive values, while being zero for negative r i .

Equivalent stress threshold
The stress thresholds r þ and r À are the maximum uniaxial equivalent stresses which can be carried within a system.The applied threshold yield surface shall be based upon a Drucker-Prager criterion, enhanced by Lubliner et al. [23] and thoroughly defined by Petracca et al. [38].The two equivalent stress surfaces in tension s À and compression s þ , are graphically displayed in Fig. 1 and discussed in the following: with the scalars a and b being used as subsequently: Furthermore, I 1 is the first invariant of the effective stress tensor and J 2 the second invariant of the deviatoric stress tensor.
r max is used as the maximal principal effective stress, respectively r min being the smallest principal effective stress.The shear reduction parameter k 1 is introduced to control the compressive strength in the shear state, which properties are shown Fig. 1.Petracca [14] and Petracca et al. [17] provide comprehensive studies on the structural effect of this parameter.Finally, the factor k b is defined by the ratio from the biaxial strength of the applied material (see Fig. 1).As the thresholds are defined for any stress state, in some cases it may happen that the damage surface advances for tension in uniaxial compression or in reverse.To avoid this, the Heavyside function HðxÞ ¼ is used to ensure following criteria: • tensile damage only evolves if at least one principal stress is positive, • compressive damage only evolves if at least one principal stress is negative.
The effect of the applied Heavyside is deciphered in Fig. 1 by the active and inactive threshold surfaces.

Damage evolution
The previously described equivalent stress thresholds shall be used as parameters for the irreversible damage evolution of the material.It is expected that the tensile softening is defined by the strength and the fracture energy with a non-linear softening, whereas the compression shows a more complex behaviour with multiple inputs for the analysis (see Fig. 2 for possible damaging zones).The tensile damage growth is obtained by the subsequent equation: where the thresholds are applied as It shall be noted that always the maximal s þ should be considered, to avoid damage relaxation.The discrete softening parameter H dis is defined as following: The characteristic length l ch is used to regularize the fracture zone and is dependent on the geometry of the finite element.Multiple approaches exist to approximate the characteristic length [10,[39][40][41], where it is possible to accurately calibrate the fracture dissipation length within structured quadrilaterals, more complex non-aligned geometries and high order form functions will typically produce discrepancies.
As proposed by Petracca et al. [17,38], the compression damage progression is analyzed on an ad hoc basis under a segmented Be ´zier curve, which is controlled by the elastic limit strength f c0 , the yield strength f cp , its corresponding yield strain cp and the residual strength f cr .The damage evolution is defined by the hardening variable RðnÞ with its strain-like counterpart n ¼ r À E (follow ''Appendix 1'' for a detailed extraction of the Be ´zier curve):

Total stresses
Including the convoluted damage parameters d þ /d À the total stresses are defined by:

Tangent constitutive tensor
With the projection operator P: the fourth order tangent constitutive tensor C tangent can be computed by: This tensor can directly be used within most finite element implementations.Checkout [42] for alternative approaches for a more consistent implementation of the constitutive tensor.

Fig. 2 Compression and tension fracture zones
In the following the recently described isotropic damage model shall be enhanced to enable orthotropic behaviours.Rather than updating the yield surface with anistropic surfaces as the Hill type criterion [29][30][31] or Tsai-Wu [32] it is preferred to base upon a mapping from the real anisotropic space to an ideal mapped isotropic stress state [33][34][35], as this is the simpler and thusly computationally more efficient approach.Furthermore, it allows the application of advanced threshold surfaces for isotropic materials, as e.g. the one discussed in Sect. 2.

Transformation between actual orthotropic and mapped isotropic space
The concept of the implicit orthotropic yield criterion was initiated by Betten [33] and lateron advanced by Oller et al. [34,35].Pela `et al. [16,36] brought those models to the application in simulation of masonry structures, with enhancing the previous models with the split in compressive and tensile parts.The models are based upon a transformation between the actual anisotropic, or in this case orthotropic stress state to an arbitrary mapped isotropic space (notation ðÁÞ Ã ) and its corresponding strain-based reverse mapping, as shown in Fig. 3.
The mapping from the actual stress state to the isotropic space is defined by the transformation tensor A r is defined as following: r Ã ¼ A r : r; ð3:1Þ considering the split into tensile r þ and compressive r À the mapping relations read: ð3:2Þ The transformation tensors A rþ and A rÀ map from the actual orthotropic space and its actual limit strength to the mapped isotropic space.As the mapped isotropic space is an independently chosen space, f þÃ and f ÀÃ can be selected arbitrary.However, Pela `L et al. [36] suggests to use f AEÃ ¼ f AE 11 1 to simplify the transformation, which are defined by following matrices according to Voigt notation (see [34] for a detailed derivation of the constituents): Additionally to f AEÃ , this model also requires f AEÃ 12 , which shall impress the mapped space shear force.Accordingly, the last entries of the mapping matrices shall devote the relation between the two shear strength. 2Considering Eq. 2.2, the stress split is still valid in the mapped isotropic space: Proceeding, the uniaxial stresses can be computed with the mapped isotropic effective stresses by Fig. 3 Mapping between actual anisotropic and isotropic space (adapted from [34]) 1 Within this section f AE 11 is the generic form for f t0;1 or f cp;1 . 2 Pela `et al. [36] considers a Drucker-Prager damage surface, which itself modifies the threshold r À 0 with the mapped strength: Þf ÀÃ .The Lubliner damage surface is reducing the uniaxial stress with the contribution of k b (see Sect. 2).This makes the mapping simpler and more consistent.
Eqs. 2.4 and 2.5 and with those the damage parameters by Eqs.2.7 and 2.9.With the updated damage parameters the total stresses per loading are defined as following: Those updated stresses are transformed back into the anisotropic space, where they yield with the appearing stress strain relation. ð3:7Þ In the actual anistropic stress space, the stress split in total stresses is equivalent to the stress split with effective stresses, as per Eq.2.2.That means that the scalars d þ and d À do not need to be mapped back, but are consistently valid in the actual space.Accordingly, the computation of the constitutive tensor as per Eq.2.12 yields similarly.

Consistent mapping of evolution law damage parameters
The major challenge of the proposed constitutive material model is the consistent mapping of all damage evolution law parameters, which shall be presented in the following.The transformation towards a mapped isotropic stress state dependent on the chosen mapping limit strength f AEÃ 1.If the material length per direction and in the mapped space are similar, (see Eq. 2.8 for the computation of l mat in the tensile regime and Eq.6.14 for the material length in the compressive regime) as: then, if no further mapping is introduced, this condition results that the fracture energies 3 in direction 1 and direction 2 would be dependent upon following restriction: To avoid this limitation, additional parameters will need a direction dependent mapping, which is dependent on the angle to the main stress h. 4 Within following association the material lengthens can therefore be interpolated consistently: 5l AEÃ mat ðhÞ ¼ l AE mat;1 cos 2 ðhÞ þ l AE mat;2 sin 2 ðhÞ: ð3:12Þ Similarly to f AEÃ , E Ã can be chosen arbitrarily.However, it again simplifies the transformation operations significantly, if it yields E Ã ¼ E 1 , as otherwise this would imply a more involved update of G Ã f , according to Eq. 3.8.Therefore, the mapped G Ã f can be computed by Pela `et al. [36]: ð3:13Þ The tensile post crack behaviour is dependent upon E, G t and f t0 , which are consistently mapped by Eq. 3.13.However, the applied compressive softening is requiring additional parameters as f c0 , f cr and cp .The stresses are given by a direct interpolation of the mapped stresses: The strain cp (see Sect. 2 and ''Appendix 1'') is mapped according to its corresponding undamaged strain i (see Fig. 13) at the limit strength, as per: following to: 3 Within this scope G AE f represents the generic fracture energy, which is either G c or G t . 4The angle to the main stress is being calculated as: cp;2 sin 2 ðhÞ: ð3:18Þ In the scenario that f ÀÃ ¼ f À 11 and E Ã ¼ E 1 , follows that Eq. 3.18 can be simplified to: The Be ´zier controllers are mapped subsequently 9 : ð3:20Þ equally as the shear reduction parameter k 1 and the factor k b defined by the ratio of biaxial strength 9 :

Numerical example
This section demonstrates the applicability of the proposed model using two numerical examples.First, this model is used to represent the directional masonry material features from [3].Next, the shear wall tests from [43] are examined with varying pre-compression and different geometries (without and with hole).Within this publication the parameter convention in the directions parallel and normal to bed joints is given in Fig. 4. E 1 corresponds to the direction e 1 and E 2 to the direction e 2 .Same accounts for the remaining parameters.e 3 refers to the out of plane direction, which is not discussed within this scope.See also [3], which refers to the same parameter definition.

Uniaxial numerical tests of [3]-masonry
In this study, the uniaxial behaviours within both directions shall be examined with the apparent material model, including a perspective to the damage growth.This is conducted to prove that the orthotropic mapping is working properly, with a focus on the variated strengths within compression and tension.
Within Table 1 are summarized the material properties including the additional Be ´zier curve parameters (see ''Appendix 1'').
The obtained results for compression and tension are presented within Fig. 5.It demonstrates that the damage evolution in both directions can be expressed appropriately.Therefore, all control parameters are mapped consistently between the apparent spaces.It shall be noted that the graphs act independently as the additional mapping is employed.If this were not provided, the relaxation would always conform dependently in the two directions.The additionally  presented damage growth shall communicate the reader an impression of how the material model is operating.Figure 5b is providing the respective selection of the tensile part.

Shear wall
This study considers a masonry wall which is first subjected to a vertical load (I) at different magnitudes.In a second step it is fixed vertically at this deformation (IIa) and sheared horizontally (IIb).Geometrically this problem appears within two different variations: once as a whole wall (see Fig. 6a) and once with a window opening (see Fi. 6b).Both walls have a dimension of 1 Â 0:99 m.
This example [43,44] is well known in literature, among others shall be mentioned: • Lourenc ¸o [45]: micro and macro models • Lourenc ¸o et al. [5]: meta models • G. Milani [46]: homogenization macro model • Pela `et al. [16,36]: orthotropic mapping damage model • Petracca et al. [10]: FE 2 model • Fu et al. [47]: orthotropic plastic damage model • Abdulla et al. [19]: simplified micro model • Bilko et al. [48]: orthotropic elastic plastic model • Pulatsu et al. [49]: discrete element model • D'Altri et al. [50]: two step adaptive limit and pushover analysis Due to the high amount of references, the example delivers great possibilities for comparisons and thus, it appears as a good benchmark.Most of the in-plane masonry phenomena are apparent: significant orthotropic material, shear failure and biaxial stiffness effects.This helps to justify if the model is suitable to express the major in-plane masonry characteristics.
It shall be noted that the apparent literature, typically focuses on either of the geometrical options: full wall or wall with window.This shows that covering this example entirely remains challenging.
The material properties of the tests are summarized in Table 2. Some of the macro scale parameters are Fig. 5 Stress-strain relation of compression and tension tests in two directions and according damage parameter growth for the apparent properties of [3].The subscript 1/2 is thereby used to denote the according direction, which is either 0 or 90 to the first material direction Fig. 6 Vermeltfoort and Raijmakers [43] shear walls.The experimental program is split into 2 phases: the pre-compression (I) then holding this deflection (IIa) and shearing (IIb) known directly from test data t ðÁÞ by Vermeltfoort and Raijmakers [43] or analyzed by micro scale analyses [45].While others are acquired through speculation s ðÁÞ and reverse engineering.These parameters are either obtained from literary works or extrapolated with the apparent demands.It should be noted that certain speculated values change across several publications since the employed models make modest adjustments and apply parameters in various manners.Unfortunately, not sufficient experimental data is available to mechanically justify either of the selected parameters.Additionally, various parameter sets could be providing the correct outcome for the apparent model.As a further challenge shall be noted that the material parameters seem to vary significantly between the experimental tests.However, within this numerical study only one parameter set shall be considered.Although there may be larger disparities compared to some of the results, this should be a step in the direction of homogenization, allowing predictions rather than just result fitting by calibration.

Wall with varying pre-pressure
The initial example is considering the wall without window, as per Fig. 6a.The experiments are discussed with 2 different lateral pre-loadings: 0.3 MPa and 1.21 MPa.As per Vermeltfoort et al. [43] the pre-loading with 0.3 MPa appears in 2 consecutive experimental trials.
The stress-strain graphs of this example are presented within Fig. 7.The strengthens of the 0.3 MPa pre-loading examples are matching well.The experimental test at 1.21 MPa indicates material characteristics that are noticeably different since the first elastic zone seems stiffer than in the previous test.That is quite usual with this heterogeneous material.In both test setups, the softening of the specimen seems Table 2 Masonry material properties from experimental tests [43] or evaluated and updated within [16,31,47]  to be faster and more significant in the numerical tests than in the experimental tests.The results of the experimental study show that the overall material behaviour can be estimated correctly.With an increased pre-loading the maximum resistance increases.While overestimating the ultimate strength for the bigger pre-compression, the overall reaction of the numerical model appears to be accurate.

Wall with hole
The second example uses the identical wall sample as the first, with the exception of a window opening.It will examine if this model can handle diverse forms.
The problem is evaluated with the same material parameters.

Quantitative analysis of stress-strain behaviour
Figure 8 displays the findings and compares them to several existing references.The graph shows good agreement between both test results.The first softening is adequately addressed by the suggested model, however the first major softening is exaggerated.
The comparison to some literature in Fig. 8 show that most results contain an almost continuous softening after reaching the peak.The presented solution is showing a jump in the softening.It is believed that this is in relation to the shear cracking which is not covered in such detail in the other references.This damaging would allow a movement, which otherwise would not be possible.However, this is based on speculation and it is suggested to conduct further research on this topic.Furthermore, it shall be mentioned that with certain calibrations this jump can be smeared.
It should be emphasised that while each model from the literature that has been cited has some merit.However, typically distinct zones are highlighted.Some results match mostly the initial limit state.Others focus more on the accurate representation of the softening.

Qualitative analysis of results
Within Fig. 9 are shown the displacements in reference to the aimed current top displacements.Furthermore are provided the tensile damage and compressive damage contour plots and the resulting maximum principal strains of the digital model.Within Fig. 11 are presented the corresponding Cauchy stresses, the principal stresses in their direction and the max and min principal stresses in the system.The colours of the principal stresses are mostly kept similar, however, peaks are highlighted slightly, which can be compared to the neighboring figures.
The evolving damage results from Fig. 9 denote crack initiations at the corners of the window, that are opening towards the supports of the walls amid increasing displacements.At the top left and basal right part of the specimen is shown a damage, which is not displayed within the referenced crack pattern.This is showing a separation of the specimen from the concrete top and bottom bars, which can also be observed within the experimental results.Furthermore, in the damage patterns can be seen that in both, the experimental tests and the numerical model, a crack happens in the left upper third and right bottom third.This proves quite effectively that a certain trust can be given to the macroscopic models and that not only micro models would display this local phenomena.Finally, the cracks which initiate at the corners of the windows do not go towards the corners of the specimen, but towards a close point at top or bottom.The reason of this is that most of the forces, which are supplied from top to bottom go through the parts left and right of the window.Accordingly those areas are largely under compression and thusly do not fail by tensile forces.Loure ´nc ¸o et al. [31] and Pela `et al. [16] are reporting a compressive damage in the upper left and lower right corner of the specimen.This is not observed in the current model.The apparent compressive stresses in those corners are higher than in the .This change is associated to different assumptions of the compressive damage criteria in previous works compared to the current research.It shall be noted that the original reports [43,44] do not specifically report this damage.
The compression damage zones are situated in regions where the tension damages are likewise apparent (see Fig. 9).This is a result of the shear reduction controlled by k b (see also [38]) and is therefore a damage which occurs due to shear in those areas.It demonstrates that the shear damage generated by this model functions as intended.The additional reduction of strength and therefore the allowance of a deflection due to shear may be one explanation to the lower compressive stresses in the corners on the opposite side of the tensile damages.
Within the Cauchy stresses and the principal stress plots can be observed that the numerical model develops a stress distribution, where most of the energy distributes around the window.The areas of the cracks do not carry any load and split the masonry wall into some structural pieces.This is a significant result of such a continuum damage analysis, as the force rearrangement might give suggestions about remaining resistances of a damaged structure.Additionally, it may indicate that other parts of the undamaged zones or supports might get higher stresses, which could lead towards a failure of connected parts or supports.
The measured displacements show that the anticipated separation into two main pieces, with some additional local failures is apparent.After a final crack opening, those pieces would mostly have a rigid body movement with a remaining sliding resistance.

Refinement and robustness of the simulation
The mesh element size and the number of elements denote important parameters in the accurateness and sometimes robustness of the simulations.Accordingly, a qualitative mesh refinement study is presented within Fig. 12.With increased refinement the appearing crack zones shrink, which is towards the expectations as the cracks are actually not zones but thin lines.Additionally, certain cracks cannot be displayed in the coarse mesh as e.g. in the 10 Â 10 mesh the middle cracks on both edges are not apparent.
It shows that the damage plots vary by an increased refinement.However, a reasonable fine mesh appears to be on the safe side.It shall be noted that refinement and regularization plays a significant role towards a robust and reliable solution.As this is a separate topic, which is valid for all types of constitutive associations, it is not be discussed further within the scope of this research.It is suggested for future research to quantitatively discuss the mesh refinement and different mesh types.
Within the discussion of this example a mesh of 50 Â 50 elements has been chosen.This is selected as some other literature appear to contain a similar refinement.Furthermore, the result plots are visually appealing, while showing all required details.However, it shall be noted that other refinement levels are equally valid, while containing small quantitative changes in the force-displacement graphs.

Conclusion and outlook
The proposed constitutive material law is being advanced from the constituents of the orthotropic mapping model introduced in [16,36] and the damage model from [14,38], which itself is based on the Lubliner threshold surface [23]  processes compared to micro models.The computational costs seem to be fairly low, whereby a comprehensive quantitative comparison would exceed the scope of this publication, as it would require a comparison to micro models with comparative implementations.
One of the improvements of the proposed model is the capability to cover the biaxial stress states and the resulting internal forces in an orthotropic manner, which allows the analysis of masonry in different preloading scenarios.Furthermore, the model has a large bandwidth of possibilities to trigger the stress strain behaviour, which brings stability within the simulation if the material response is known.As contrast, comparatively many parameters need to be obtained, which is up-to date infrequent, as typically only few parameters are studied within a series of experimental tests.It has to be noted that a large amount of material properties needed to be speculated as not sufficient test data was provided.Seeing the performance of the law within a fully comprehensive test setup [3] would be required to specify the entire generality of the model.It shall be noted that the model contains 17 parameters per direction, making it 34 in total.Unfortunately, mostly barely that much information is known about the material.Once the uniaxial evolution law would be simplified, this number could be reduced.However, this would limit the possibilities for a correct adjustment of the model.
The model was studied within examples which cover uniaxial compression and tension, shear, biaxial stress states with a simple and a more complex geometry.The test cases have been in-plane 2D problems in a plane-stress state, while the expected behaviour could be simulated properly.However, more complex shapes would need to be considered to allow a more generic evaluation of the proposed model with the behaviour of the sophisticated construction material masonry.Accordingly, an avenue for future research shall be the application of the constitutive law within 3D solid formulations for complex shapes with largely varying thickness.Another possibility would be the imposition within thin-walled or thick/solidbased shell formulations to investigate the bending and the out-of plane performance.Thereby, shells with a thickness integration would consider multiple entities of the material per in-plane integration point (see [51,52] for eventual applications).

Declarations
Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix 1: Uniaxial compression Be ´zier-curve based relaxation
Within this appendix, the uniaxial compression Be ´zier-curve based relaxation introduced by Petracca et al. [17], which is shown within Fig. 13 shall be deciphered.The model itself requires as input E, f c0 , f cp , cp and f cr .Furthermore, the locations of the remaining required stress/strains can be controlled by a set of scalars, namely the Be ´zier controllers c 1 , c 2 and c 3 .However, if preferred, the input can also contain the remaining stresses and strains, rather than controllers.
With the input of the Be ´zier controller inputs, the complete set of initial control parameters shall be defined as following: 13 Uniaxial compression segmented Be ´zier-curve based relaxation (adapted from [17]) Within, the damage model the strain-like counterpart shall be defined as following: which is the basis for the evaluation of the hardening variable RðnÞ: RðnÞ ¼ 0 for n 0 Bðn; 0 ; i ; p ; r 0 ; r i ; r p Þ for 0 \n p Bðn; p ; j ; k ; r p ; r j ; r k Þ for p \n k Bðn; k ; u ; r ; r k ; r u ; r r Þ for k \n r r u for r \n 8 > > > > > > < > > > > > > : ð6:7Þ whereby, each Be ´zier curve segment is defined as following: BðX; x 1 ; x 2 ; x 3 ; including the the parameters: Furthermore, the fracture energy G c , which is the shaded area under the Be ´zier curve (see Fig. 13) can be computed as following: r p p 2 G c;2 ¼ Að p ; j ; k ; r p ; r j ; r k Þ G c;3 ¼ Að k ; u ; r ; r k ; r u ; r r Þ ð6:10Þ whereby A is the respective area under the Be ´zier curve: Aðx 1 ; x 2 ; x 3 ;y 1 ; y 2 ; ð6:11Þ Considering that the explained model is an ad hoc formulation, G c is not the total area as it will need to be regularized to comply with the grey shaded area, being G c l dis (see Fig. 13).The regularization is applied with a stretching factor S, defined subsequently: The stretching factor S is applied on the constituents j ; k ; u and r : n ¼ n þ Sð n À p Þ for n ¼ j; k; uandr ð6:13Þ It shall be noted that the stretching factor should not be lower than -1, which would otherwise lead into snap back and accordingly a discontinuity in the softening.
To avoid S\ À 1, the dissipated length should be smaller than the respective plastic material length:

Fig. 4
Fig. 4 Parameter definition within this scope

Fig. 9
Fig. 9 Observed displacements, damage, and strains at [43]-wall with window tests.Computed with 50 Â 50 elements . The major contribution of this work lays in the development of a consistent mapping of the parameters of the respective damage evolution laws.It shows great potential with the examples from Sect. 4 as it provides correct responses on the load capacities of the studied masonry structures.Furthermore, as a continuum based macro model, it is easier within the modelling

Fig. 10
Fig. 10 Observed damage at experimental tests of [43]-wall with window in comparison to the numerical solutions.Red is the damage from the first experimental test, blue from the second.(Colour figure online)

Acknowledgements
The authors thank the reviewers whose constructive comments improved the quality of the article.The authors gratefully acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) as part of the project ''Eine Methode zur effizienten Simulation großer Mauerwerksscheiben unter exzentrischer und/oder zyklischer biaxialer Beanspruchung auf Grundlagewirklichkeitsnaher Kleinko ¨rperversuche'' (BL306/34-1) and the Australian Government Research Training Program Scholarship and The University of Queensland as part of the ''Global Strategy and Partnerships Seed Funding Scheme, 2019''.Data availability All simulations have been processed with Kratos Multiphysics [53, 54] and the implementation of the corresponding constitutive law is available under a open source license.The pre-and post-processing has been done with the Rhino plugin Cocodrilo [55, 56] by employing the CADintegrated analysis workflow from [57].