An Orthotropic Elastic-Plastic Constitutive Model for Masonry Walls

Abstract

The use of a continuum structural model for the analysis of masonry structures in the plane stress state is discussed in this paper. Attention is paid to orthotropic masonry at the material level and validation of the model after its implementation in a proprietary finite element method (FEM) system via user-supplied subroutine. The constitutive relations are established in the framework of the mathematical elastoplasticity theory of small displacements and deformations. Based on the orthotropic failure criterion that was originally proposed by Hoffman in the spatial stress state, the model includes a generalization of the criterion in the plane stress. As it is the case for isotropic quasi-brittle materials, different yield surfaces are considered for tension and compression, which are both of Hoffman type.

1. Introduction

The mechanical behavior of masonry is subjected to various influencing factors, mostly resulting from its complicated mesoscopic and microscopic structure and two basic materials used. Therefore, different modeling approaches are available for the numerical simulation of the mechanical behavior of masonry structures. Extensive research has been conducted on the development of advanced numerical modeling and the analysis of historical masonry structures for several decades [1,2,3,4]. However, a macroscopic approach is required for the viable analysis and the prediction of global structural behavior. The macroscopic composite behavior of masonry can be described assuming a homogeneous continuum and an anisotropic material with directional properties. This alternative of the constitutive modeling proved to be promising in two-dimensional problems, especially for models with closed form of the failure, yield, or limit surface. Some constitutive macro-models that are relied on for the finite element method (FEM) are widely used in the last decades with a different degree of complexity and idealizations—from initial models with simply isotropic, linear elastic behavior to advanced models with non-linear orthotropic behavior that are recently developed in the framework of modern concepts of continuum damage mechanics.

More considerable interest in the biaxially loaded masonry began more than forty years ago, in the late seventies of the last century, both by experiments and theory. In general, experiments were mainly carried out at that time concerning the failure of shear walls. They were usually related to the proposed failure criterion for the masonry in the plane state of stresses in the representative form of the particular tests that had to be involved [5,6]. The representative biaxial test for walls built with any kind of masonry elements was conducted very seldom because of the technical difficulties. The biaxial test data of Page [7,8] were interpolated in [9] in the form of a failure criterion. The criterion was written using three failure surfaces in the form of elliptic cones expressed by a second-order polynomial of stress components in the reference axes coaxial with masonry layers and without any reference to the observed distinct failure modes of masonry panels. On the other hand, an important composite failure criterion which was derived based on the four different failure mechanisms was proposed in [10]. This criterion was next used in the FEM-based constitutive model in [11].

There have been few attempts to use a single failure surface in constitutive models of masonry panels because of the non-acceptable fit of experimental values. The general orthotropic failure criterion of Tsai-Wu [12] was already available for composite materials since 1971. The use of this criterion for masonry was attempted in [13] in the form of a polynomial function of the second degree in stresses. A criterion in the form of a double pyramid with a rectangular base and the slope angle equal to the internal friction angle of the material was assumed in [14]. The non-acceptable fit of Page’s experimental values resulting from the Hoffman [15] single surface criterion was discussed in [16], although the criterion itself as a single limit condition is quite flexible and attractive to use. The phenomenological single-surface constitutive model cannot also distinguish between different failure modes. In the framework of computational plasticity, the use of a composite yield criterion containing several failure criteria seems to overcome this drawback. The choice of two failure surfaces, one for the tension regime and the second one for the compression regime, provides better agreement with Page’s experimental results. To describe the orthotropic behavior of masonry, two orthotropic yield criteria were used in [17], where a hydrostatic pressure insensitive material of a Hill-type as the yield criterion for compression together with a Rankine-type yield criterion for tension were assumed. The composite criterion containing three failure surfaces, each of them being for tension, compression, and shear regimes, allows distinguishing between different failure modes [18]. Contours consisting of a few surfaces are characterized good compatibility with the experimental results—see, e.g., [19,20] and [18,21]. However, this approach seems to be too demanding a computational task of plasticity.

Concerning masonry inelastic behavior, the closed-form macro-models are more efficient and suitable for complex structural computations. Some constitutive models that are recently developed in the framework of modern concepts of continuum damage mechanics are based on the assumption that the masonry axes of the bed and head joints are also damage principal axes. Usually, scalar damage parameters are assumed in each direction of the fixed axis, see for instance [14], where two independent scalar parameters in each direction of the material axis were used and their evolution is described by the energy-based approach. A similar approach in the framework of computational plasticity was used in [17], where the principal directions of damage are fixed and aligned with the initial orthotropy axes and softening/hardening relationships were adopted for the stress–strain diagrams in tension and compression, with different fracture energies along the axes of each material. The energies were coupled by a single scalar internal parameter used in the plasticity algorithm to measure simultaneously the amount of softening/hardening in two material axes. The main drawback of the closed-form orthotropic macro-models is the identification of the material parameters. To estimate macro-scale properties from mortar and brick parameters and their bonding, homogenization techniques can be used, both in elastic as well as plastic behavior, e.g., [22,23,24]. An alternative option is to transfer the identification problem to the level of masonry constituents by using multi-scale methods, e.g., [18,20,21], or to use methods developed recently for modeling of anisotropic quasi-brittle fracture, e.g., the so-called phase-field fracture model in the diffusive damage mechanics [25,26,27].

A continuum damage model in which the orthotropic behavior is simulated using a mapping relationship between the orthotropic behavior and an auxiliary model has been recently published in [28]. By using the concept of the mapped stress tensor the problem can be more efficiently solved in the mapped space and the results can be transported to the real field. Two distinct isotropic failure criteria are assumed in the mapped space and two stress transformation tensors are adopted. In the paper, the computational representation of complex failure loci obtained by experiments on orthotropic masonry is also presented.

This paper discusses an extension of the approach presented in [17]. Although the extension can be done both in a plane stress state as well as in a spatial stress state, the constitutive model with the generalization of the Hoffman criterion in a plane stress state is discussed here. Two orthotropic failure criteria are used, which are formulated in the framework of the representation theory of orthotropic tensor functions based on the Hoffman criterion [15]. The use of two Hoffman-type failure criteria as the yield criteria in the plasticity model seems to be particularly attractive and may give the better fit of the experimental values [29,30]. The composite masonry is treated as a homogenized orthotropic continuum. Since the failure criteria are scalar-valued functions of the stress tensor, the invariant representations of these criteria are dependent only on orthotropic invariants of the stress tensors. It is also the purpose of the paper to show the possibility of formulating robust numerical algorithms of the model implementation into a commercial finite element code at the integration point level using user-defined subroutines. Some tests of the proposed numerical algorithm for an anisotropic continuum are presented in the paper, both at the single element level and at the structural level and in the plane stress state.

2. The Orthotropic Hoffman-Type Failure Criteria in an Invariant Form

To model such effects as a marked difference observed between strengths in tension and compression, Hoffman [15] proposed a fracture criterion for brittle orthotropic materials as an extension to the Hill yield criterion. The criterion was proposed in the spatial stress state and in the $\left\{{\mathbf{m}}_i\right\}$ frame of reference that coincides with the axes of orthotropy. The criterion was originally described by the function with nine material constants $C_i$ that are dependent on the three uniaxial tensile strengths $Y_{ti}$ and the three uniaxial compressive strengths $Y_{ci}$, along with the orthotropy directions i and also on the three shearing strengths $k_{ij}$ on the planes of material orthotropy. In the case of the plane stress state, when the normal vector to the stress plane is coincided with the axis of orthotropy $\mathbf{n}={\mathbf{m}}_3={\mathbf{b}}_3$ (Figure 1), stress components ${\sigma }_{33}={\sigma }_{31}={\sigma }_{23}=0$ and the Hoffman criterion:

$$\begin{array}{c}{C}_1{\left({\sigma }_{22}-{\sigma }_{33}\right)}^2+C_2{\left({\sigma }_{33}-{\sigma }_{11}\right)}^2+C_3{\left({\sigma }_{11}-{\sigma }_{22}\right)}^2+\\ C_4{\sigma }_{11}+C_5{\sigma }_{22}+C_6{\sigma }_{33}+C_7{\sigma }_{23}^2+C_8{\sigma }_{13}^2+C_9{\sigma }_{12}^2-1=0,\end{array}$$

(1)

where:

$$\begin{array}{c}{C}_1=\frac{1}{2}\left(\frac{1}{Y_{t2}{Y}_{c2}}+\frac{1}{Y_{t3}{Y}_{c3}}-\frac{1}{Y_{t1}{Y}_{c1}}\right),C_2=\frac{1}{2}\left(\frac{1}{Y_{t3}{Y}_{c3}}+\frac{1}{Y_{t1}{Y}_{c1}}-\frac{1}{Y_{t2}{Y}_{c2}}\right),\\ C_3=\frac{1}{2}\left(\frac{1}{Y_{t1}{Y}_{c1}}+\frac{1}{Y_{t2}{Y}_{c2}}-\frac{1}{Y_{t3}{Y}_{c3}}\right),C_4=\frac{1}{Y_{t1}}-\frac{1}{Y_{c1}},C_5=\frac{1}{Y_{t2}}-\frac{1}{Y_{c2}},\\ C_6=\frac{1}{Y_{t3}}-\frac{1}{Y_{c3}},C_7=\frac{1}{k_{23}^2},C_8=\frac{1}{k_{31}^2},C_9=\frac{1}{k_{12}^2},\end{array}$$

(2)

takes the following form with the six material constants:

$$C_1{\left({\sigma }_{22}\right)}^2+C_2{\left({\sigma }_{11}\right)}^2+C_3{\left({\sigma }_{11}-{\sigma }_{22}\right)}^2+C_4{\sigma }_{11}+C_5{\sigma }_{22}+C_9{\sigma }_{12}^2-1=0,$$

(3)

where the constants $C_1÷C_3$ are dependent on uniaxial strengths $Y_{t3}$ and $Y_{c3}$ in the direction perpendicular to the stress plane.

The theory of tensor functions together with the theorems on their representations has been recognized to be an efficient mathematical tool for the formulation of constitutive relationships with both the desirable analytical clarity and required generality.

For some other recent applications of tensor functions see, e.g., [31,32]. It also allows accounting straightforwardly the invariance requirements of the principle of the space isotropy and the material symmetries so that the orientation of the material in space does not affect on its constitutive relation. Using this theory with Boehler’s results [33], we can assume that the orthotropic criterion (3) is a particular case of the more general scalar-valued orthotropic function of three invariants $tr\left({\mathbf{M}}_1\mathsf{\sigma }\right)$, $tr\left({\mathbf{M}}_2\mathsf{\sigma }\right)$, $tr{\mathsf{\sigma }}^2$ of the following form:

$$f\left(tr\left({\mathbf{M}}_1\mathsf{\sigma }\right),tr\left({\mathbf{M}}_2\mathsf{\sigma }\right),tr{\mathsf{\sigma }}^2\right)-1=0,$$

(4)

where $\mathsf{\sigma }$ is the symmetric plane stress tensor ($\mathsf{\sigma }\in T_2^s,dim{T}_2^s=3$) and ${\mathbf{M}}_{\alpha }$ are the parametric (structural) tensors defined as:

$${\mathbf{M}}_1={\mathbf{m}}_1\otimes {\mathbf{m}}_1,{\mathbf{M}}_2={\mathbf{m}}_2\otimes {\mathbf{m}}_2.$$

(5)

The unit vectors ${\mathbf{m}}_{\alpha }$ are the privileged directions of the orthotropic material, so they have to be perpendicular to each other. The invariants are very useful for the interpretation of the failure surface in any coordinate systems of the plane stress tensor that are different from the principal axes of orthotropy. Following the paper [34] or [35], we can choose another set of invariants $K_p$ in the form:

$$K_1=tr\left({\mathbf{M}}_1\mathsf{\sigma }\right),K_2=tr\left({\mathbf{M}}_2\mathsf{\sigma }\right),K_3=tr{\mathsf{\sigma }}^2-{\left(tr\left({\mathbf{M}}_1\mathsf{\sigma }\right)\right)}^2-{\left(tr\left({\mathbf{M}}_2\mathsf{\sigma }\right)\right)}^2,$$

(6)

where the symbol “$tr$” denotes the trace of a second order tensor ($tr\left(\mathbf{AB}\right)=A_{ij}{B}_{ij}$). The form (6) is very convenient because in the $\left\{{\mathbf{m}}_{\alpha }\right\}$ axes the invariants are:

$$K_1={\sigma }_{11},K_2={\sigma }_{22},K_3={\left(\sqrt{2}{\sigma }_{12}\right)}^2.$$

(7)

Using the invariants (7), the criterion (3) can be treated as a particular case of the criterion proposed in [34] and may be written in the following invariant form:

$$f\left(K_i\right)-1=a_{\alpha }{K}_{\alpha }+b_{\alpha \beta }{K}_{\alpha }{K}_{\beta }+c{K}_3-1=0,$$

(8)

where $\alpha ,\beta =1,2$, $i=1,2,3$ and the material constants are defined as:

$$\begin{array}{c}{a}_{\alpha }=\frac{1}{Y_{t\alpha }}-\frac{1}{Y_{c\alpha }},b_{11}=\frac{1}{Y_{t1}{Y}_{c1}},b_{22}=\frac{1}{Y_{t2}{Y}_{c2}},\\ b_{12}=\frac{1}{2}\left(\frac{1}{Y_{t3}{Y}_{c3}}-\frac{1}{Y_{t1}{Y}_{c1}}-\frac{1}{Y_{t2}{Y}_{c2}}\right),c=\frac{1}{2k_{12}^2}.\end{array}$$

(9)

Note that in the constant $b_{12}$ there are again uniaxial strengths $Y_{t3}$ and $Y_{c3}$ in the direction perpendicular to the stress plane.

The criterion (8) can be also written in the following invariant form:

$$\frac{1}{2}\mathsf{\sigma }·\mathbf{P}·\mathsf{\sigma }+\mathbf{p}·\mathsf{\sigma }-1=0,$$

(10)

where a dot means a double contraction of two tensors, $\mathbf{p}$ is the symmetric tensor function of the second-order:

$$\mathbf{p}=a_1{\mathbf{M}}_1+a_2{\mathbf{M}}_2,$$

(11)

$\mathbf{P}$ is the double symmetric tensor function of the fourth-order:

$$\mathbf{P}=2b_{11}{\mathbf{M}}_1\otimes {\mathbf{M}}_1+2b_{22}{\mathbf{M}}_2\otimes {\mathbf{M}}_2+2b_{12}\left({\mathbf{m}}_1\otimes {\mathbf{m}}_2+{\mathbf{m}}_2\otimes {\mathbf{m}}_1\right)+c\mathbb{M}$$

(12)

with the fourth-order tensor $\mathbb{M}$:

$$\mathbb{M}=4\mathbf{N}\otimes \mathbf{N},\mathbf{N}=\frac{1}{2}\left({\mathbf{m}}_1\otimes {\mathbf{m}}_2+{\mathbf{m}}_2\otimes {\mathbf{m}}_1\right).$$

(13)

Several criteria proposed in the literature for orthotropic materials are special cases of the quadratic limit surface (8), including an elliptic failure surface according to Tsai and Wu [12] and criteria discussed recently in [30]. However, a phenomenological single-surface model may give an insufficient description of the mechanical behavior. It does not permit easy identification of failure modes and thus renders the description of different post-failure mechanisms very difficult. At least two failure criteria should be taken into consideration, the one for the compression regime and the second for tension regime. Each of them may be of the form (8) as proposed in [34] where the failure criterion for orthotropic materials in the spatial stress state is represented by two quadratic functions of the six invariants of the stress tensor and parametric tensors. It may more accurately describe the failure data distribution than classical limit surfaces, although it may include fifteen independent material parameters for the description of failure surfaces. On the other hand, the concept of a smooth single-surface description seems to be attractive from a numerical point of view and also allows for modeling of different inelastic behavior by changing size, shape, and location of a quadratic state function according to the form (8) in orthotropic stress space.

It is possible to propose a generalization of the Hoffman criterion for the orthotropic material in the plane state of stresses. The determination of the six material parameters in the criterion (3) requires six strength tests. The five standard tests are uniaxial loading along the axes of orthotropy (two tests for tension and two for compression) and the shearing test in the plane of stresses. If the test of the uniform biaxial compression is used for determining the strength $Y_{cc}$ as the sixth test in addition to conventional tests, we will get the criterion in the invariant form (8) or (10), in which the material parameter $b_{12}$ is changed to the form:

$$b_{12}^{\left(1\right)}=\frac{1}{2Y_{cc}^2}-\frac{1}{2Y_{t1}{Y}_{c1}}-\frac{1}{2Y_{t2}{Y}_{c2}}+\frac{1}{2Y_{cc}}\left(\frac{1}{Y_{t1}}+\frac{1}{Y_{t2}}-\frac{1}{Y_{c1}}-\frac{1}{Y_{c2}}\right).$$

(14)

If as the additional test we will use the test of the uniform biaxial tension for determining the strength $Y_{tt}$, we will get the criterion in the invariant form with the following material parameter $b_{12}$:

$$b_{12}^{\left(2\right)}=\frac{1}{2Y_{tt}^2}-\frac{1}{2Y_{t1}{Y}_{c1}}-\frac{1}{2Y_{t2}{Y}_{c2}}-\frac{1}{2Y_{tt}}\left(\frac{1}{Y_{t1}}+\frac{1}{Y_{t2}}-\frac{1}{Y_{c1}}-\frac{1}{Y_{c2}}\right).$$

(15)

Note that now there are not the uniaxial strengths $Y_{t3}$ and $Y_{c3}$ in the material constant $b_{12}$. Finally, using two functions of the form (8) or (10), one with the parameter (14) for compression regime $f_1\left(K_i\right)-1=0$ and the second with the parameter (15) for tension regime $f_2\left(K_i\right)-1=0$, we can construct the composite failure criterion in such a way that the following set:

$$B=B_1\cap B_2\wedge B_{\alpha }\equiv \left\{\mathsf{\sigma }\in T_2^s|f_{\alpha }\left(K_i\right)-1<0\right\},$$

(16)

is convex and

$$f_{\alpha }\left(K_i\right)-1=a_{\alpha }^{\left(\alpha \right)}{K}_{\alpha }+b_{\alpha \beta }^{\left(\alpha \right)}{K}_{\alpha }{K}_{\beta }+c^{\left(\alpha \right)}{K}_3-1=0,$$

(17)

We can also assume that $c^{\left(1\right)}=c^{\left(2\right)}\equiv c$ and $a_{\alpha }^{\left(\alpha \right)}=a_{\alpha }$, $b_{\alpha \alpha }^{\left(\alpha \right)}=b_{\alpha \alpha }$ for simplification and we then have the following seven material parameters: the two uniaxial tensile strengths $(Y_{t\alpha }>0)$, the two uniaxial compressive strengths $(Y_{c\alpha }>0)$, the pure shear strength $(k_{12}>0)$, the biaxial uniform compressive strength $Y_{cc}>0$ and the biaxial uniform tensile strength $Y_{tt}>0$. This procedure is proposed in the paper [34], although the number of material parameters may be increased to 12 even though some of them can be used just to adjust failure surfaces to experimental data.

Figure 2 shows the criterion in the orthotropy axes with the following parameters for material: simple shear strength $k=0.4$ [MPa], vertical compression strength $Y_{c2}=6.34$ [MPa] and vertical tensile strength $Y_{t1}=0.31$ [MPa] perpendicular to the horizontal joint, compression strength $Y_{c1}=4.95$ [MPa] and tensile strength $Y_{t2}=0.14$ [MPa] parallel to the horizontal joint, biaxial uniform compression strength $Y_{cc}=9.0$ [MPa] and biaxial uniform tension $Y_{tt}=0.14$ [MPa]. Figure 2a shows general view and Figure 2b cross section of ${\sigma }_{12}=0$ (solid line) with contours every $0.2$ [MPa] (dashed lines). Figure 2c shows cross section of ${\sigma }_{22}=-6.34$ [MPa] (solid line) with parallel contours every $0.81$ [MPa] (dashed lines). Figure 2d shows cross section of ${\sigma }_{11}=-4.95$ [MPa] (solid line) with parallel contours every $0.81$ [MPa] (dashed lines).

An alternative connection of two surfaces is proposed, different from that used in [34], which is shown schematically in Figure 3. It allows for shifting a common edge. Figure 3 shows the method for constructing a boundary surface of the criterion in the normal stress components, in the orthotropy axes (at zero shear stresses). The surface in the tension range conventionally adopts compressive strength values $Y_{c\alpha }^t$ greater than twice those for the area of the compression range $Y_{c\alpha }^c$, that is $Y_{c\alpha }^t=2Y_{c\alpha }^c$. On the other hand, the surface in the compression range conventionally adopts tensile strength $Y_{t\alpha }^c$, by making them the value of compressive strength $Y_{c\alpha }^c$. It is assumed that they will be two-fold smaller, that is $Y_{t\alpha }^c=0.5Y_{c\alpha }^c$.

3. Comparison with the Experimental Results

The composite failure criterion is shown in Figure 4 in the orthotropy axes and in comparison with one of the most complete sets of experimental data of biaxially loaded masonry that was given by Page [7,8], who tested 102 panels of half-scale solid clay brick masonry. Tested elements were made on a scale of 1:2 with dimensions 360 × 360 × 50 [mm]. Tests were differentiated due to the rotation of the principal stress terms of the horizontal joint (axis of the material). The ratio of the principal stresses was changed so that the wall was considered in any possible state of stress. The following values of the parameters are adopted: axial tensile strength along the horizontal joint $Y_{t1}=0.43$ [MPa] and perpendicular to it $Y_{t2}=0.32$ [MPa], the axial compression strength along the horizontal joint $Y_{c1}=8.74$ [MPa] and perpendicular to it $Y_{c2}=8.03$ [MPa], the pure shear strength $k=0.33$ [MPa], uniform biaxial tensile strength $Y_{tt}=0.32$ [MPa] and compression strength $Y_{cc}=8.38$ [MPa]. Good agreement was found in the shape of the failure surface for principal stresses, which may mean that the criterion would appear to be sufficiently well validated for further investigations. Finally, the Rankine–Hill failure criterion [16] is shown in Figure 4 by dashed lines as a comparison to the proposed criterion.

Another discussed example is the comparison of the proposed failure criterion with the results of the experimental research from work [36]. The research program was carried out at the ETH Polytechnic in Zurich. Ganz and Thürliman’s team investigated biaxially loaded wall panels (designated K1–K12) with dimensions of 1200 × 1200 × 150 [mm${}^3$] and with different orientation of material axes to the principal stress directions. The results of experimental tests [36] are presented in

Table 1. Comparison of the proposed criterion with the experimental results of [36].

			Experimental			Numerical
Panel	${\mathit{\sigma }}_1/{\mathit{\sigma }}_2$	$\mathit{\varphi }$	[MPa]			[MPa]			Exp/Num
		[°]	${\mathit{\sigma }}_{11}$	${\mathit{\sigma }}_{22}$	${\mathit{\tau }}_{12}$	${\mathit{\sigma }}_{11}$	${\mathit{\sigma }}_{22}$	${\mathit{\tau }}_{12}$	[-]
K1	−0.09	22.5	$-0.08$	$-0.92$	$0.42$	$-0.09$	$-1.04$	$0.47$	$0.88$
K2	−0.05	22.5	$-0.17$	$-1.42$	$0.62$	$-0.17$	$-1.35$	$0.59$	$1.05$
K3	0.0	0.0	$0.00$	$-7.63$	$0.00$	$0.00$	$-7.61$	$0.00$	$1.01$
K4	0.0	90.0	$-1.83$	$0.00$	$0.00$	$-1.87$	$0.00$	$0.00$	$0.98$
K6	0.0	45.0	$-0.32$	$-0.32$	$0.32$	$-0.43$	$-0.43$	$0.43$	$0.74$
K7	0.0	22.5	$-0.39$	$-2.25$	$0.93$	$-0.37$	$-2.15$	$0.89$	$1.05$
K8	0.0	67.5	$-0.22$	$-0.04$	$0.09$	$-0.22$	$-0.04$	$0.08$	$1.01$
K10	0.33	0.0	$-2.11$	$-6.44$	$0.00$	$-2.02$	$-6.15$	$0.00$	$1.05$
K11	0.31	22.5	$-2.04$	$-4.49$	$1.23$	$-2.02$	$-4.46$	$1.22$	$1.01$
K12	0.30	45.0	$-2.03$	$-2.03$	$1.08$	$-2.01$	$-2.01$	$1.07$	$1.01$

excluding two tests (K5 and K9), which concerned samples of reinforced walls. The second column of the Table 1 contains the proportion of stresses, while the third column gives the angle by which the system of material axes has been rotated to the load directions. The ratio of principal stresses allows the direction of the load path to be determined. The intersection point of the load path direction with the criterion surface determines the stress state, the components of which are placed in columns 7–9 of the Table 1. The criterion surfaces are determined by the following parameters [MPa]: for a tension regime—$Y_{t1}=0.28$, $Y_{t2}=0.01$, $Y_{c1}=3.74$, $Y_{c2}=15.72$, $k_t=0.048$, $Y_{tt}=0.01$ and for a compression regime—$Y_{t1}=0.94$, $Y_{t2}=3.81$, $Y_{c1}=1.87$, $Y_{c2}=7.61$, $k_c=2.868$, $Y_{cc}=2.06$. The anisotropy ratio of the compressive strength is $Y_{c2}/Y_{c1}=4.07$ and is related to the arrangement of cores in ceramic masonry elements. The length of the load vector should be determined as $\sqrt{{\sigma }_{11}^2+{\sigma }_{22}^2+2{\sigma }_{12}^2}$. The last column of the Table 1 presents the values of the ratios of the length of the vectors resulting from the experiment and from the initial surfaces of the proposed criterion. The results show good compliance of the criterion with the experimental results.

Figure 5a shows the initial surface of the criterion with points marked on the surface, indicating the stress limits of individual experimental tests conducted by Ganz and Thürliman. Analyzing the contour lines, it can be seen that most of the safe area is limited by the tensile regime. This is more clearly seen in Figure 5b,c, in which the plane cross-sections of the criterion are shown. One can see cross-sections with assigned to them, the values of the ultimate stress obtained from the individual experimental tests which are marked with diamonds.

At the ETH Polytechnic in Zurich, a test program for walls made of concrete, hollow masonry elements ZSW1-ZSW12 was also carried out [37]. The results of the experimental tests are given in

Table 2. Comparison of the proposed criterion with the experimental results of [37].

			Experimental			Calculated
Panel	${\mathit{\sigma }}_1/\phantom{{\mathit{\sigma }}_1{\mathit{\sigma }}_2}{\mathit{\sigma }}_2$	$\mathit{\varphi }$	[MPa]			[MPa]			Exp/Num
		[°]	${\mathit{\sigma }}_{11}$	${\mathit{\sigma }}_{22}$	${\mathit{\tau }}_{12}$	${\mathit{\sigma }}_{11}$	${\mathit{\sigma }}_{22}$	${\mathit{\tau }}_{12}$	[-]
ZSW1	0.0	0.0	$0.00$	$-9.12$	$0.00$	$0.00$	$-9.21$	$0.00$	$0.99$
ZSW2	0.14	0.0	$-6.12$	$-0.83$	$0.00$	$-6.05$	$-0.82$	$0.00$	$1.01$
ZSW4	1.53	0.0	$-5.98$	$-9.13$	$0.00$	$-5.58$	$-8.52$	$0.00$	$1.07$
ZSW5	0.0	45.0	$-3.06$	$-3.06$	$3.06$	$-3.06$	$-3.06$	$3.06$	$1.00$
ZSW6	0.22	45.0	$-4.60$	$-4.60$	$2.93$	$-4.69$	$-4.69$	$2.98$	$0.98$
ZSW7	1.0	45.0	$-6.12$	$-6.12$	$0.00$	$-6.36$	$-6.36$	$0.00$	$0.96$
ZSW8	0.0	67.5	$-2.34$	$-0.40$	$0.97$	$-2.34$	$-0.40$	$0.97$	$1.00$
ZSW9	0.0	22.5	$-0.97$	$-5.66$	$2.35$	$-0.97$	$-5.66$	$2.35$	$1.00$

, and in Figure 6. Figure 6a shows failure surface of the criterion with the experimental results of [37]. Figure 6b,c shows cross-sections with planes of the constant normal stress. A position of the yield stress from experimental tests are marked with diamonds.

The criterion surfaces are determined by the following parameters [MPa]: for a tension regime—$Y_{t1}=0.01$, $Y_{t2}=0.01$, $Y_{c1}=11.52$, $Y_{c2}=18.42$, $k_t=0.01$, $Y_{tt}=0.01$ and for a compression regime—$Y_{t1}=2.88$, $Y_{t2}=4.61$, $Y_{c1}=5.76$, $Y_{c2}=9.21$, $k_c=3.98$, $Y_{cc}=6.36$. The anisotropy ratio of the compressive strength is $Y_{c2}/Y_{c1}=1.60$. The material has almost zero tensile strength. As before, the load vector length was determined for each limit point as $\sqrt{{\sigma }_{11}^2+{\sigma }_{22}^2+2{\sigma }_{12}^2}$. The maximum difference of vectors resulting from the experiment data and those of the initial surface compared criterion does not exceed 7 percentage points, which shows a very good agreement of the proposed model with the experiment results.

4. The Constitutive Relations and Implementation of the Model

The elastic-plastic orthotropic material is considered assuming an additive decomposition of the strain tensor into the elastic part ${\mathsf{\epsilon }}_e$ and the plastic part ${\mathsf{\epsilon }}_p$. The elastic part is defined by orthotropic Hooke’s law. The plastic part of the strain tensor is defined by a flow rule associated with the yield function given by the plasticity (failure) criterion written based on the forms (10) and (17) as:

$$f_{\alpha }\left(\mathsf{\sigma },z_{\alpha }\right)=\frac{1}{2}\mathsf{\sigma }·{\mathbf{P}}_{\alpha }·\mathsf{\sigma }+{\mathbf{p}}_{\alpha }·\mathsf{\sigma }-K_{\alpha }\left(z_{\alpha }\right)=0,$$

(18)

where the $K_{\alpha }\left(z_{\alpha }\right)$ are given functions with the real functional value from a closed interval [0, 1] that describes the type of hardening/softening (Figure 7), $z_{\alpha }$ are internal scalar hardening variables. The softening behavior is modeled with a smeared approach, where the localized damage is represented by the scalar, which is related by an equivalent length h to the released energy per unit cracked area, $G_f$. The length h should correspond to a dimension of the finite element mesh. As one can see in Figure 7, different fracture energies are introduced in the model, as additional parameters—the tensile fracture energy $G_{ft}$ and the compressive energy $G_{fc}$.

In the frame of reference coinciding with the orthotropy axes, we have the following matrix representations of tensors ${\mathbf{P}}_{\alpha }$ and ${\mathbf{p}}_{\alpha }$ in the Voigt notations for the plane stress $\mathsf{\sigma }\Rightarrow {\left[{\sigma }_{11}{\sigma }_{22}\sqrt{2}{\sigma }_{12}\right]}^{\mathrm{T}},$ ($\alpha =t$ for tension surface, $\alpha =c$ for compression surface):

$${\mathbf{P}}_{\alpha }\Rightarrow {\mathbf{P}}_{\alpha }=\left[\begin{array}{ccc}2b_{11}^{\alpha }\hfill & 2b_{12}^{\alpha }\hfill & 0\hfill \\ 2b_{21}^{\alpha }\hfill & 2b_{22}^{\alpha }\hfill & 0\hfill \\ 0\hfill & 0\hfill & 4c^{\alpha }\hfill \end{array}\right],{\mathbf{p}}_{\alpha }\Rightarrow {\mathit{p}}_{\alpha }=\left[\begin{array}{ccc}{a}_1^{\alpha }\hfill & a_2^{\alpha }\hfill & 0\hfill \end{array}\right].$$

(19)

For composite criteria, the subscript $\alpha$ also denotes the number of the active surface. Let us first assume that only one surface is active, which will allow this marking to be omitted. In the framework of the mathematical theory of the elastic-plastic material a permissible stress state is any state of stresses for which $f\le 0$. A stress state is called the elastic stress state if $f<0$. A plastic state refers to a stress state at the boundaries of the current elastic region for which $f=0$. The plastic part of the strain tensor is defined by a flow rule associated with the yield function given by the Equation (18). The flow rule defines the sign (direction) of plastic-strain increment in the following form:

$${\dot{\mathsf{\epsilon }}}_p=\gamma {\left.\frac{\partial f\left(\mathsf{\sigma },z\right)}{\partial \mathsf{\sigma }}\right|}_{\mathsf{\sigma }={\mathsf{\sigma }}^T}=\gamma \left(\mathbf{P}·\mathsf{\sigma }+\mathbf{p}\right)\equiv \gamma \mathbf{r},$$

(20)

where $\gamma >0$ is a plastic multiplier. After applying differentiation to orthotropic Hooke’s elastic law with the respect to the time and after substituting Equation (20) we obtain:

$$\dot{\mathsf{\sigma }}=\mathsf{C}·(\dot{\mathsf{\epsilon }}-\gamma \mathbf{r})\equiv {\mathsf{C}}_{ep}·\dot{\mathsf{\epsilon }},$$

(21)

where the elasto-plastic tangent operator ${\mathsf{C}}_{ep}$ can be calculated after the parameter $\gamma$ is known. Assuming that

$$\dot{z}=\gamma (\mathbf{r}·\mathbf{r})\equiv \gamma {∥\mathbf{r}∥}^2.$$

(22)

one can compute from the consistency condition

$$\gamma \dot{f}\left(\mathsf{\sigma },z\right)=0,\gamma >0,$$

(23)

the plastic multiplier

$$\gamma =\frac{\langle \mathbf{r}·\mathsf{C}·\dot{\mathsf{\epsilon }}\rangle }{\mathbf{r}·\mathsf{C}·\mathbf{r}+\partial K/\partial z{∥\mathbf{r}∥}^2},$$

(24)

and the operator ${\mathsf{C}}_{ep}$ in the following form:

$${\mathsf{C}}_{ep}=\mathsf{C}-\frac{\left(\mathsf{C}·\mathbf{r}\right)\otimes \left(\mathbf{r}·\mathsf{C}\right)}{\mathbf{r}·\mathsf{C}·\mathbf{r}+\partial K/\partial z{∥\mathbf{r}∥}^2}.$$

(25)

The double symmetric fourth-order tensor of elastic material constants $\mathsf{C}$ in the orthotropic Hooke’s law can be conveniently defined by the compliance tensor $\mathsf{S}\equiv {\mathsf{C}}^{-1}$, which in the frame of reference aligned with the orthotropic axes can be written in Voigt notation for the plane stress as

$$\mathsf{S}\to \left[\begin{array}{ccc}\frac{1}{E_1}& -\frac{{\nu }_{21}}{E_2}& 0\\ -\frac{{\nu }_{12}}{E_1}& \frac{1}{E_2}& 0\\ 0& 0& \frac{1}{G_{12}}\end{array}\right].$$

(26)

where we have five technical in-plane moduli: $E_1,E_2$ are Young’s moduli, $G_{12}$ is the shear modulus and ${\nu }_{12},{\nu }_{21}$ are Poisson’s ratios (${\nu }_{12}/E_1={\nu }_{21}/E_2$). It should be noted that two yield criteria are combined in the model into a composite yield surface and the intersection of different yield surfaces defines corners that require special attention in a numerical algorithm according to Koiter’s generalization.

4.1. Implementation into Finite Elements under the Plane Stress Condition

In this subsection, according to the convention adopted in many finite element programs the components of a symmetric second-order tensor are presented as a single column array, whereas fourth-order tensors are presented as two-dimensional arrays. The matrix representations of the tensors are shown in terms of the Cartesian components in the frame coinciding with the materialaxes of orthotropy.

The constitutive relationship in (21) is in the form of the “highly non-linear” differential equation which can be solved by the modified Euler method (usually the implicit Euler backward algorithm). Therefore, it is replaced by the incremental equation of the form:

$$\mathsf{\Delta }\sigma =\mathbf{C}(\mathsf{\Delta }\mathit{\epsilon }-\gamma \tilde{\mathbf{r}})\equiv {\tilde{\mathbf{C}}}_{ep}\mathsf{\Delta }\mathit{\epsilon },$$

(27)

where ${\tilde{\mathbf{C}}}_{ep}$ is called the operator consistent with the integration algorithm of constitutive relations. We assume that for each $t_n\in [0,T]$ the strain increment $\mathsf{\Delta }\mathit{\epsilon }={\mathit{\epsilon }}_{n+1}-{\mathit{\epsilon }}_n$ is known, thus the problem is strain driven, and we want to compute the stress state ${\mathit{\sigma }}_{n+1}$ for $t_{n+1}$. We assume:

$${\mathit{\sigma }}_{n+1}={\mathit{\sigma }}_{n+1}^{trl}-\mathsf{\Delta }\gamma \mathit{C}{\mathit{r}}^{trl},$$

(28)

where ${\mathit{\sigma }}_{n+1}^{trl}=\mathit{C}{\mathit{\epsilon }}_{n+1}$ is called the trial stress state and ${\mathit{r}}^{trl}$ is the gradient of $f\left({\mathit{\sigma }}_{n+1}^{trl},{\mathit{z}}_n\right)$. The calculation of the multiplier $\mathsf{\Delta }\gamma >0$ and the tensor function $\tilde{\mathit{r}}$ (Equation (27)) is significantly dependent on the realization of ${\mathit{\sigma }}_{n+1}^{trl}$. A detailed description of the numerical implementation into the FEM is given in [38] for the elastic-plastic material with the yield criteria of Huber–Mises–Hencky (isotropy) and Hill (orthotropy). The most important step is to calculate the multiplier using the quadratic equation of the variable $\mathsf{\Delta }\gamma$. This step is significantly different from the case of the Hill yield criterion.

Based on the consistency condition (23) in the algorithmic form:

$$f\left({\mathit{\sigma }}_{n+1}^T,z_{n+1}\right)=\frac{1}{2}{\mathit{\sigma }}_{n+1}^T\mathit{P}{\mathit{\sigma }}_{n+1}+{\mathit{p}}^T{\mathit{\sigma }}_{n+1}-K\left(z_{n+1}\right)=0,$$

(29)

we obtain after the substitution of (28) into (29) a following quadratic equation (indices $n+1$ suppressed):

$$A\mathsf{\Delta }{\gamma }^2+B\mathsf{\Delta }\gamma +C=0,$$

(30)

where:

$$\begin{array}{c}A=\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{trl}\right)}^T\mathit{P}\mathit{C}{\mathit{r}}_{n+1}^{trl}-\frac{{\partial }^{2}K}{\partial z^2}{\left[{\left({\mathit{r}}_{n+1}^{trl}\right)}^T{\mathit{r}}_{n+1}^{trl}\right]}^2,\hfill \\ B=-\left[\frac{1}{2}{\mathit{\sigma }}_{n+1}^{trl}\mathit{P}\mathit{C}{\mathit{r}}_{n+1}^{trl}+\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{trl}\right)}^T\mathit{P}{\mathit{\sigma }}_{n+1}^{trl}+{\mathit{p}}^T\mathit{C}{\mathit{r}}_{n+1}^{trl}+\frac{\partial K}{\partial z}{\left({\mathit{r}}_{n+1}^{trl}\right)}^T{\mathit{r}}_{n+1}^{trl}\right],\hfill \\ C=f^{trl}\left({\sigma }_{n+1}^{trl},z_n\right).\hfill \end{array}$$

(31)

The solution of the Equation (30) is particularly simple if the hardening/softening function $K\left(z\right)$ is the second degree polynomial of the form as shown in the Figure 7. The last part of (30) may become then equal to the yield function of the trial state $f^{trl}$. It can be seen that the linearization of the Equation (30) does not lead to significant errors. Determining the plastic multiplier from the linear part of Equation (30) in the form $\mathsf{\Delta }\gamma =-B/C$ does not lead to large errors, although it depends on the assumed strain increment step length.

A model based on two criteria complicates the procedure algorithm. If only one of the criteria is exceeded, e.g., the stresses are reduced to the exceeded surface according to the algorithm (28)–(30), as for a model based on only one criterion. It remains, to establish which criterion is actually exceeded.

A separate case is exceeding both conditions at the same time:

$$f_1^{trl}>0\wedge f_2^{trl}>0.$$

This is the case where the point is outside the edge joining both surfaces and requires a different procedure. A linear combination of gradients has been used here (compare Koiter’s law). The component of the plastic multiplier is calculated separately for the tensile criterion and for the compression criterion. Plastic strains and stress state components at time $t_{n+1}$ can be determined on the basis of the formulas:

$$\begin{array}{c}{\mathit{\epsilon }}_{n+1}^p={\mathit{\epsilon }}_n^p+\mathsf{\Delta }{\gamma }_1{\mathit{r}}_{n+1}^{1,trl}+\mathsf{\Delta }{\gamma }_2{\mathit{r}}_{n+1}^{2,trl},\\ {\mathit{\sigma }}_{n+1}={\mathit{\sigma }}_{n+1}^{trl}-\mathsf{\Delta }{\gamma }_1\mathit{C}{\mathit{r}}_{n+1}^{1,trl}-\mathsf{\Delta }{\gamma }_2\mathit{C}{\mathit{r}}_{n+1}^{2,trl},\\ {\mathit{C}}_{ep}=\mathit{C}-\frac{\mathit{C}{\mathit{r}}_{n+1}^{1,trl}{\left({\mathit{r}}_{n+1}^{1,trl}\right)}^T{\mathit{C}}^T}{{\left({\mathit{r}}_{n+1}^{1,trl}\right)}^T\mathit{C}{\mathit{r}}_{n+1}^{1,trl}+\partial K_1/\partial z_1{∥{\mathit{r}}_{n+1}^{1,trl}∥}^2}-\frac{\mathit{C}{\mathit{r}}_{n+1}^{2,trl}{\left({\mathit{r}}_{n+1}^{2,trl}\right)}^T{\mathit{C}}^T}{{\left({\mathit{r}}_{n+1}^{2,trl}\right)}^T\mathit{C}{\mathit{r}}_{n+1}^{2,trl}+\partial K_2/\partial z_2{∥{\mathit{r}}_{n+1}^{2,trl}∥}^2}.\end{array}$$

(32)

After similar transformations as before, the following system of quadratic equations is obtained.

$$\left\{\begin{array}{c}{A}_1{\left(\mathsf{\Delta }{\gamma }_1\right)}^2+B_1\mathsf{\Delta }{\gamma }_1+C_1+D_1{\left(\mathsf{\Delta }{\gamma }_2\right)}^2+E_1\mathsf{\Delta }{\gamma }_2+F_1\mathsf{\Delta }{\gamma }_1\mathsf{\Delta }{\gamma }_2=0\\ A_2{\left(\mathsf{\Delta }{\gamma }_2\right)}^2+B_2\mathsf{\Delta }{\gamma }_2+C_2+D_2{\left(\mathsf{\Delta }{\gamma }_1\right)}^2+E_2\mathsf{\Delta }{\gamma }_1+F_2\mathsf{\Delta }{\gamma }_1\mathsf{\Delta }{\gamma }_2=0\end{array}\right.,$$

(33)

where some of the coefficients are analogous to the problem with one active surface and are marked with letters $A_{\alpha },B_{\alpha },C_{\alpha }$ ($\alpha =1,2$):

$$\begin{array}{cc}\hfill A_{\alpha }& =\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{\alpha ,trl}\right)}^T{\mathit{P}}_{\alpha }\mathit{C}{\mathit{r}}_{n+1}^{\alpha ,trl}-K_{\alpha ,n}^{″}{\left[{\left({\mathit{r}}_{n+1}^{\alpha ,trl}\right)}^T{\mathit{r}}_{n+1}^{\alpha ,trl}\right]}^2,\hfill \\ \hfill B_{\alpha }& =-\frac{1}{2}{\mathit{\sigma }}_{n+1}^{trl}{\mathit{P}}_{\alpha }\mathit{C}{\mathit{r}}_{n+1}^{\alpha ,trl}-\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{\alpha ,trl}\right)}^T{{\mathit{P}}_{\alpha }\mathit{\sigma }}_{n+1}^{trl}-{\left({\mathit{p}}_{\alpha }\right)}^T\mathit{C}{\mathit{r}}_{n+1}^{\alpha ,trl}-K_{\alpha ,n}^{\prime }{\left({\mathit{r}}_{n+1}^{\alpha ,trl}\right)}^T{\mathit{r}}_{n+1}^{\alpha ,trl},\hfill \\ \hfill C_{\alpha }& =f_{\alpha }^{trl}\left({\mathit{\sigma }}_{n+1}^{trl},z_n^{\alpha }\right).\hfill \end{array}$$

(34)

The remaining coefficients depend on the data related to both boundary surfaces and the formulas can be expressed as follows:

$$\begin{array}{cc}\hfill D_1& =\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{2,trl}\right)}^T{\mathit{P}}_1\mathit{C}{\mathit{r}}_{n+1}^{2,trl},D_2=\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{1,trl}\right)}^T{\mathit{P}}_2\mathit{C}{\mathit{r}}_{n+1}^{1,trl},\hfill \\ \hfill E_1& =-\frac{1}{2}{\mathit{\sigma }}_{n+1}^{trl}{\mathit{P}}_1\mathit{C}{\mathit{r}}_{n+1}^{2,trl}-\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{2,trl}\right)}^T{\mathit{P}}_1{\mathit{\sigma }}_{n+1}^{trl}-{\mathit{p}}_1^T\mathit{C}{\mathit{r}}_{n+1}^{2,trl},\hfill \\ \hfill E_2& =-\frac{1}{2}{\mathit{\sigma }}_{n+1}^{trl}{\mathit{P}}_2\mathit{C}{\mathit{r}}_{n+1}^{1,trl}-\frac{1}{2}{\left(\mathit{C}{\mathit{r}}_{n+1}^{1,trl}\right)}^T{\mathit{P}}_2{\mathit{\sigma }}_{n+1}^{trl}-{\mathit{p}}_2^T\mathit{C}{\mathit{r}}_{n+1}^{1,trl},\hfill \\ \hfill F_1& ={\left(\mathit{C}{\mathit{r}}_{n+1}^{1,trl}\right)}^T{\mathit{P}}_1\mathit{C}{\mathit{r}}_{n+1}^{2,trl},F_2={\left(\mathit{C}{\mathit{r}}_{n+1}^{2,trl}\right)}^T{\mathit{P}}_2\mathit{C}{\mathit{r}}_{n+1}^{1,trl}.\hfill \end{array}$$

(35)

According to the above numerical algorithm, the several models with the failure criterion of the type (18) has been coded in the programming language FORTRAN and next implemented into a commercial finite element code DIANA [39]. There are standard Newton–Raphson and Riks algorithms for the solving nonlinear equilibrium equations in the program. Because of this, implementation of the model into a finite element code is done at the integration point level by means of user-defined subroutine USRMAT. The subroutine lets the user specify a general nonlinear material behavior by updating the state variables over the equilibrium step $n\to n+1$ within a framework of an incremental–iterative algorithm of finite element method. Both the return-mapping algorithm allowing the stresses to be returned to the yield surface and a consistent tangent stiffness operator have been coded. The implementation of the model is presented both at the single element tests (see Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12) and at the structural level test (next section).

4.2. Results of the Single-Element Tests

Tests in the plane stress were performed in the homogeneous stress state and with one isoparametric continuum finite element (four-node, linear interpolation and Gaussian integration). The displacement-controlled load diagram is shown in Figure 8. The following material data are adopted: the moduli of elasticity in two directions of orthotropy: $E_1$ = $E_2$ = 8 [GPa], the shear modulus: $G_{12}=3.478$ [GPa], the Poissonś coefficient: ${\nu }_{12}=0.15$, that is, as for the isotropic material. The strength parameters associated with the initial tensile failure surface of the criterion are: $Y_{c1}=17$, $Y_{c}2=17$, $Y_{t1}=0.35$, $Y_{t2}=0.25$, $Y_{tt}=0.22$, $k_{12}=0.296$ [MPa]. The strength parameters associated with the compressive surface are: $Y_{c1}=8.5$, $Y_{c2}=8.5$, $Y_{t1}=8.5$, $Y_{t2}=8.5$, $Y_{tt}=8.5$, $k_{12}=4.9$ [MPa]. The material parameters are typical for unreinforced walls in terms of the values and the proportion between them. The strength degradation curves are adopted with two internal scalar parameters: $z_t$ for the tensile softening and $z_c$ for the softening during compression (Figure 7). The curves are matched to set the fracture energies along the first orthotropic material axis as $G_{f1}=54.0$ [J/m${}^2$] during tension and as $G_{fc1}=20.0$ [kJ/m${}^2$] during compression.

In the test with the homogeneous deformation field, all eight degrees of freedom in the four-node finite element were fixed to force the desired linear deformation. To compare the behavior of our model in the tests, the model of the Rankine–Hill criterion was used, which is standard in the DIANA system and dedicated to the analysis of masonry structures. The following parameters of the Rankine–Hill model are adopted: $Y_{c1}=8.5Y$, $Y_{c2}=8.5Y$, $\beta$ = $-1.0$, $\gamma =3.0$, $Y_{t1}=0.35Y$, $Y_{t2}=0.25Y$, $\alpha =1.0$, where $Y=1$ [MPa] and the tensile fracture energies $G_{fX}=54.0$ [J/m${}^2$], $G_{fY}=18.0$ [J/m${}^2$] and the compressive fracture energies $G_{fcX}=20.0$ [kJ/m${}^2$] and $G_{fcY}=15.0$ [kJ/m${}^2$]. The algorithm of the own model was also programmed in MATHEMATICA. As a result, it was possible to control the correctness of the algorithm and its implementation in both computing environments. Results of the tests are presented in Figure 9, Figure 10, Figure 11 and Figure 12.

5. Tests at the Structural Level

Two examples of the proposed constitutive model are briefly presented here as the tests at the structural level. More information on the tests can be found in the works [41,42].

The first example is restricted to the numerical simulation of the load capacity tests of masonry structures in the plane stress state that were conducted experimentally in [43]. In laboratory tests, two similar wall fragments marked as J2G and J3G were tested, Figure 13. The masonry shear walls with an opening and the thickness of 10 mm were built from 18 courses of masonry cement-sand units 210 × 52 × 100 [mm]. The top and bottom courses were fully clamped in steel beams. The wall is shearing with the horizontal force F under displacement control. The top edge can move with a horizontal displacement (Figure 13). The wall was first vertically loaded through a steel beam to the value p = 0.3 kN/m${}^2$ that remains constant through the subsequent loading steps of the horizontal force up to the maximum horizontal displacement of $\mathsf{\Delta }$ = 16 mm. The beam movement was limited to the horizontal direction, i.e., the lower and upper edges of the wall were kept parallel. Figure 14 shows the forms of cracks obtained in the tests.

Cracks appeared diagonally between the opening and the lower and upper corners of the wall. In addition, tensile cracks appeared on the outer vertical edges of the wall, on both sides of the opening at the top of the left pillar and at the bottom of the right pillar, and their development ran from the outside to the center of the wall. The failure mechanism resulting from the wall tests is shown schematically in the Figure 14. As can be seen, the kinematics of the system focuses on the movement of four blocks connected by hinges. Due to the development of the material crushing zone, also marked in the drawing, the mechanism will activate the compression criterion.

For the numerical simulation of the wall failure using the model based on the own criterion, a 1989 mesh of four-node flat finite elements with an average side length of 23 [mm] was built. The values of material parameters used in the numerical simulation are shown in the

Table 3. Masonry shear wall. Material properties of the model.

Elastic Moduli
	$E_{11}$	$E_{22}$	${\nu }_{12}$	$G_{12}$
	[MPa]	[MPa]	−	[MPa]
	7520	3960	0.09	1460
Uniaxial, Biaxial and Shear Strengths in the Orthotropic Axes [MPa]
	$Y_{c1}^{\alpha }$	$Y_{c2}^{\alpha }$	$Y_{t1}^{\alpha }$	$Y_{t2}^{\alpha }$	$Y_{\alpha \alpha }$	$k_{12}$
Compression $\alpha =c$	$5.25$	$3.75$	$2.625$	$1.825$	$3.0$	$0.45$
Tension $\alpha =t$	$11.0$	$7.5$	$0.35$	$0.25$	$0.30$	$0.30$
Fracture energies in J/m${}^2$		$G_{fc1}$ = 2350		$G_{ft1}$ = 43.3

. The values given in the table were adopted on the basis of work [44], where the data were adopted based on the results of experimental tests [43], and the supplementary parameters result from the homogenization procedures and were taken from work [16].

The comparison between numerical and experimental load–displacement diagrams, for wall J2G and J3G, is given in Figure 15. Apart from own calculation, results obtained by Pelá in [44] are also shown. Good agreement is found in the elastic range and satisfactory agreement in the inelastic range, although slightly worse in terms of the load capacity than the results obtained in [44]. The maximum divergence is around 17%. It is possible to better match the results of the experiment, however it requires additional calculations and time-consuming calibration of strength parameters and fracture energies. The behavior of the wall at the horizontal displacement $\mathsf{\Delta }=12$ mm is depicted in Figure 16 and Figure 17 in terms of the maximum and minimum principal plastic strains, respectively. Both tensile and compressive failure zones are captured by the model. This indicates that the wall deformability and the general mechanism of its destruction are correctly reflected.

Figure 16 shows the principal plastic strains ${\epsilon }_1^p$ obtained from own simulations. They can be interpreted as the distribution of material failure zones as a result of exceeding the criterion in tension. These are also zones of development of cracks in the structure. Figure 16a shows the strain results in the form of maps plotted on a deformed finite element mesh. For a better effect, the deformation has been scaled 10 times.

In Figure 16b, the results in the vector form were superimposed on the mesh, wich allows to compare the length of the strain and the direction. Figure 17a shows the principal plastic deformations ${\epsilon }_3^p$, which in Figure 17b are shown in vector form. Areas with large ${\epsilon }_3^p$ strains can be interpreted as zones of material failure due to possible material crushing.

In the second example, the behavior of the masonry infill wall that is built within a reinforced concrete frame is numerically simulated (Figure 18). The wall was experimentally tested in [45]. The frame and the wall were first subjected to the vertical loads $P_2=97.8$ kN and $P_3=48.9$ kN that remain constant through the subsequent loading steps of the horizontal force $P_1$ up to the failure.

For the numerical analyses 4-noded quadrilaterals, linear plane stress, and continuum elements are utilized. The subsequent loading steps of the horizontal force up to the maximum horizontal displacement $\mathsf{\Delta }$ = 20 mm were analyzed under a displacement control. The material properties of the model are given in

Table 4. Characteristics of materials for orthotropic model of masonry.

Elastic Moduli
	$E_{11}$, MPa	$E_{22}$, MPa	${\nu }_{12}={\nu }_{21}$	$G_{12}$
	[MPa]	[MPa]	[-]	[MPa]
	13,790	13,790	0.16	3480
Characteristics of the Strength Surface in MPa
	$Y_{c1}^{\alpha }$	$Y_{c2}^{\alpha }$	$Y_{t1}^{\alpha }$	$Y_{t2}^{\alpha }$	$Y_{\alpha \alpha }$	$k_{12}^{\alpha }$
Compression ($\alpha =c$)	20.7	20.7	20.7	20.7	23.00	5.0
Tension ($\alpha =c$)	40.0	40.0	1.38	1.38	1.2	0.5
Fracture energies in J/m${}^2$		$G_{fc1}$ = 350		$G_{ft1}$ = 16.7

. The comparison between numerical and experimental load–displacement diagrams is given in Figure 19. The low initial vertical load combined with the confinement provided by the reinforced concrete frame yields extremely ductile behavior.

In Figure 20, a crack zone is plotted as the map of the tensile principal plastic strain. A good agreement is found with respect to the calculated collapse load value. The model gives the response that is close to the experimental equilibrium path. At the ultimate stage, a well-defined failure mechanism is formed with a final diagonal shear band going from one corner of the specimen to the other.

The precise localization of cracks can also be achieved using models based on a smeared crack approach [46,47], fracture-based models with advanced crack tracking techniques [44] or by using micro-modeling [16]. To the evaluation of historic buildings for dynamic actions methods based on ultimate load-bearing capacity are used. Especially in combination with the monitoring based on operational modal analysis (OMA) techniques, it is becoming a popular practice [48]. An alternative may also be the relatively young topic of multi-scale modeling [49] or to describe diffuse cracks a gradient-enhanced damage model based on nonlocal displacements and the extended finite element method (XFEM) for sharp cracks [50,51]. However, if we compare the numerical result presented here with the failure mechanisms shown in Figure 16 and Figure 20, we find a good agreement and confirm the ability of our own model to correctly reflect the behavior of the structure.

6. Conclusions

The mechanical behavior of the continuum material model can be described using constitutive relations based on the theory of tensor functions. This theory, together with the theorems on the representation of tensor functions, constitute an effective mathematical tool for the formulation of constitutive relations of the orthotropic material. The invariance requirements of the isotropy of space and the orthotropic symmetry of materials are easily considered, so that the orientation of the material in space does not affect its constitutive relation. It is possible to obtain the analytical transparency as well as to maintain the required universalism of constitutive equations, see, e.g., [30]. The composite orthotropic failure criterion of the proposed model is constructed from two square scalar functions that are dependent on three orthotropic invariants of the plane stress tensor. In general, the criterion needs to have the twelve material parameters to define the failure surface. However, in practice, only seven of them are used that are obtained from the appropriate uniaxial, biaxial and shear strength tests. The criterion is an example of the orthotropic failure criteria, which can be treated as a generalization of the well-known Hoffman failure criterion that is often used for brittle orthotropic materials in which compressive strengths and tensile strengths are significantly different.

The numerical tests confirm correctness of the implementation and the ability of the models to reproduce failure modes in the structural tests in certain situations. They also show that it is possible to incorporate strain softening into the proposed class of models with a single yield surface. At present, the implementation of the other models within the framework of multisurface plasticity is being tested. Although for the multisurface plasticity the intersection of the different yield surfaces defines corners that require special attention in a numerical algorithm, it has the advantage of engaging different hardening laws for each surface, which might be more physically realistic due to the distinction of tension and compression regimes. The choice of a hardening parameter is not crucial because the available experimental data are scarce.

The commercial version of Diana with the so-called user-supplied subroutine mechanism can be used as the development environment for computational materials research. We have used the user-supplied subroutine USRMAT to implement new material models for continuum spatial and plane stress elements. The new models can be applied to a variety of structural problems, to single element tests but also to simulate physical experiments, using different elements types and using standard features of the program such as advanced solution procedures, for instance indirect displacement control with full Newton–Raphson. Also, the use of the post-processing capabilities of the program is an advantage, although the post-processing of the user-defined status variables might be more user-friendly.