Fast brick-based homogenized limit analysis for in- and out-of-plane loaded periodic masonry panels

Computers and Structures 231 (2020) 106206

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Fast brick-based homogenized limit analysis for in- and out-of-plane loaded periodic masonry panels Simone Tiberti, Gabriele Milani ⇑ Department of Architecture, Built Environment and Construction Engineering, Technical University of Milan, Piazza Leonardo Da Vinci, 20133 Milan, Italy

a r t i c l e

i n f o

Article history: Received 19 July 2019 Accepted 13 January 2020

Keywords: Masonry Homogenization Limit analysis Kirchhoff-Love plates In-plane behavior Out-of-plane behavior

a b s t r a c t This paper presents a 3D brick-based model that aims at the derivation of in- and out-of-plane homogenized failure surfaces for masonry. The considered masonry panel is discretized into regular parallelepiped 3D finite elements that are supposed to be rigid; the model here introduced uses a Kirchhoff-Love plate kinematics for the out-of-plane description of the displacement rate field of the elements. A linear programming problem formulated in standard form is scripted into Matlab to derive the in- and out-of-plane homogenized failure surfaces, also enabling the extraction of failure modes for the considered masonry element. The constraints of the linear programming problem come from the combination of an upper bound limit analysis problem and a homogenization-based approach. The proposed model is validated for two separate case studies: a running bond masonry test-window and an English bond masonry test-window. The homogenized failure surfaces resulting from the current model show good correspondence to those presented in three distinct works available in literature. Also, a few relevant failure modes are derived for the two case studies, and they are consistent with the expected deformed shapes at collapse for their related load conditions. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction Over the last quarter of century, homogenization has established itself as a strong, reliable tool to model the behavior of masonry. Based on the concept of Representative Element of Volume (in short, REV), which is the smallest cell of material containing all the characteristics needed for its full representation, homogenization represents a good in-between with respect to the other two modelling strategies usually employed for masonry: macro-modelling and micro-modelling. The former substitutes the composite material consisting of unit and mortar with an equivalent homogeneous material, which can incorporate orthotropy [1] and damage [2]; the latter uses separate models for units and mortar [3–6], and sometimes even for the physical interfaces between them [7]. In fact, homogenization counters the critical issues of both strategies: it does not need extensive experimental tests to identify the material properties, yet it still enables to derive the characteristics of masonry with little computational effort starting from the accurate micromechanical model of the REV, which takes

⇑ Corresponding author. E-mail address: [email protected] (G. Milani). https://doi.org/10.1016/j.compstruc.2020.106206 0045-7949/Ó 2020 Elsevier Ltd. All rights reserved.

into account the units’ arrangement and the mechanical properties of the constituents. Homogenization has been widely used for deriving the elastic properties of masonry [8,9], as well as for modeling the postpeak behavior of masonry aiming at performing non-linear analyses [10–12] or at including the presence of damage and degradation [13–15]. A specific field of application for homogenization consists of its combination with the two theorems of limit analysis. This represents an interesting development because it allows the derivation of the masonry behavior at collapse, pairing the swiftness of homogenization to the results that can be achieved by limit analysis in the determination of collapse loads/stress distributions and/or mechanisms of failure. This path has been extensively explored in several works, mainly devoted to the assessment of the collapse behavior of masonry under in-plane load conditions. A series of two papers by Milani and co-workers [16,17] presents a micro-mechanical approach for the homogenized limit analysis of masonry walls subjected to in-plane loads. The macroscopic resistance to in-plane load conditions is there assessed by extracting homogenized inplane failure surfaces for the considered masonry elements, in the trail of both the upper and lower bound theorems of limit analysis. Further practical applications are presented in a solo work by Milani [18]. A different approach is used in a work by Milani and

2

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

Taliercio [19], which is based on the Method of Cells that is usually employed for fiber-reinforced composites. Other approaches are presented by Stefanou and co-workers [20] and Godio and coworkers [21]; the latter employs the Cosserat continuum theory for modeling masonry. Finally, some works deal with the extraction of in-plane homogenized failure surfaces for non-periodic masonry, where the units are assembled in the masonry element with different rates of randomness. Homogenization can still be exploited in these cases, which apparently undermine the very application of this technique, provided that statistical tools are used for selecting a suitable REV. A work presented by Milani and Lourenço [22] employs a method based on the equivalence between a discrete 3D block model and a 2D continuum model, and makes use of Monte Carlo simulations to derive the homogenized failure surfaces for a huge number of random assemblages. Conversely, Tiberti and Milani [23] introduce a full 2D model that is based on the so-called ‘‘pixel strategy”: an accurate finite element mesh is directly generated from the sketch of the considered non-periodic masonry panel, and is obtained by converting each pixel into an element. A simple procedure is implemented for a quick identification of a statistical REV, and the in-plane homogenized failure surfaces are eventually derived by solving a linear programming problem that corresponds to an upper bound limit analysis problem coupled to homogenization. When it comes to the out-of-plane homogenized failure surfaces for masonry walls, the available literature is not as abundant as for the in-plane case. Very few approaches have been presented over the years, most of them addressing the out-ofplane behavior with ad-hoc simplifications. Milani and coworkers [24] use an approach that considers the masonry cell under investigation as consisting of several layers over the thickness. Sab and co-workers [25] employ a procedure that is based on a Reissner-Mindlin model for determining the yield strength domain of periodic brick masonry, which expands upon a LoveKirchhoff plate model. A series of papers by Cecchi and Milani employs a procedure named ‘‘kinematic identification”: a 3D system of blocks connected by interfaces is kinematically identified with a 2D Reissner-Mindlin plate, for which the out-of-plane homogenized failure surfaces are extracted despite not directly using a rigorous homogenization approach. Applications are oriented on a running bond masonry cell [26] and an English bond masonry cell [27]. Milani and Taliercio [28] use a simplified kinematics and again ground their work on the aforementioned Method of Cells. This brief literature review highlights the absence of a full 3D model combining limit analysis and homogenization. This paper presents a homogenized limit analysis model that is capable of tackling both the in- and the out-of-plane collapse behavior of masonry, which represents its main novelty. A Kirchhoff-Love plate model is used for representing the kinematic field of the considered REV, and a Matlab script is employed for the derivation of both in- and out-of-plane homogenized failure surfaces, as well as the extraction of the deformed configurations at collapse (‘‘failure modes”). Regarding the structure of this paper, Section 2 offers a detailed presentation of the upper bound limit analysis problem as formulated for this application, with a complete insight on the kinematic model here employed and the constraints of the linear programming problem. The procedure for deriving the homogenized failure surfaces and the failure modes is also described in this section. Section 3 presents the results obtained for a running bond masonry REV and an English bond masonry REV, which are validated against data available in literature. Specifically, the results are compared in terms of homogenized failure surfaces, also commenting relevant failure modes obtained for the chosen REVs.

2. Homogenized upper bound limit analysis problem formulation This section presents the mathematical formulation of the problem that aims at deriving homogenized failure surfaces for masonry panels. In particular, the out-of-plane behavior of the masonry panel is assessed introducing a Kirchhoff-Love plate model into the displacement rate field. An upper bound limit analysis problem is combined with a homogenized approach and is formulated as a standard form linear programming problem in Matlab [29]. For the purposes of this problem, any investigated masonry panel must be discretized into a mesh consisting of rigid, regular parallelepiped 3D elements bereft of rotation rate (/_ ¼ 0). This enables to handle a higher number of elements and entails an overall reduction of the problem unknowns. 2.1. Kinematics of Kirchhoff-Love plate model Considering the elements as rigid and devoid of rotation rate, the kinematics of each element is completely determined by the
 displacement rate field of its centroid u_ x ; u_ y ; u_ z , where axis Z represents the transversal direction. As anticipated in the previous section, a Kirchhoff-Love plate model is here used for a full representation of each element’s kinematics, so that the three components can be expressed as:

u_ x ¼ u_ x;per þ E_ xx xG þ E_ xy yG þ v_ xx zG xG þ 0:5v_ xy zG yG

ð1Þ

u_ y ¼ u_ y;per þ E_ xy xG þ E_ yy yG þ v_ yy zG yG þ 0:5v_ xy zG xG

ð2Þ

u_ z ¼ u_ z;per  0:5v_ xx xG 2  0:5v_ yy yG 2  0:5v_ xy xG yG

ð3Þ

where fxG ; yG ; zG g are the coordinates of the element’s centroid with respect to a reference system located at the center of the masonry
 panel, u_ x;per ; u_ y;per ; u_ z;per are the periodic velocities of the element, n o E_ xx ; E_ xy ; E_ yy the components of the average strain rate tensor (with n o E_ yx equal to E_ xy for symmetry), and v_ xx ; v_ xy ; v_ yy the components of the average curvature rate tensor (again, with v_ yx equal to v_ xy for symmetry). The choice of such kinematics is compatible with the hypothesis of a Kirchhoff-Love plate model; the requirements of 2 2 @ 2 u_ z v_ xx ¼  @@xu_2z , v_ yy ¼  @@yu_2z , and v_ xy ¼ 2 @x@y are all satisfied. 2.2. Velocity jumps and plastic flow constraints The hypothesis of rigid elements means that plastic dissipation is possible only across their mutual interfaces, which also implies a reduction in terms of unknowns. Moreover, the use of regular parallelepiped elements entails that their six sides are already oriented according to the reference system chosen for the masonry panel. The problem is further simplified by introducing three accessory hypotheses:  no dissipation occurs across interfaces that are orthogonal to the Z axis;  no shear dissipation occurs along the Z axis of all the other interfaces;  no condition is enforced with respect to the elements’ velocity along the Z direction. Therefore, the only active interfaces are those orthogonal to the X or Y axis, and the shear component along the Z axis is neglected. Each active interface can then only experience tangential and normal

3

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

Fig. 1. Velocity jumps for a generic interface normal to axis X.

Fig. 3. Mohr-Coulomb failure criterion.

Since the expressions of Eq. (9) are already oriented along the reference system chosen for the masonry panel, they are not subject to any modification. Therefore, the velocity jumps are simply expressed as: 2 X

q @F q k_ I ¼ @ rn q¼1

Du_ n ¼

Du_ t ¼

Fig. 2. Velocity jumps for a generic interface normal to axis Y.

velocity jumps, which are evaluated according to Eqs. (4), (5) or (6), (7) depending on the considered interface (see Figs. 1 and 2):

    j i Du_ n  Du_ x ¼ u_ x;per  u_ x;per þ E_ xx xjG  xiG þ v_ xx zjG xjG  ziG xiG

ð4Þ

    j i Du_ t  Du_ y ¼ u_ y;per  u_ y;per þ E_ xy xjG  xiG þ 0:5v_ xy zjG xjG  ziG xiG ð5Þ     l k Du_ n  Du_ y ¼ u_ y;per  u_ y;per þ E_ yy ylG  ykG þ v_ yy zlG ylG  zkG ykG

Du_ t  Du_ x ¼

l u_ x;per



k u_ x;per

2 X

q @F q k_ I ¼ @s q¼1

q 1 2 k_ I Aqn ¼ k_ I tan / þ k_ I tan /

ð10Þ

q¼1 2 X

q 1 2 k_ I Aqt ¼ k_ I  k_ I

ð11Þ

q¼1

where the velocity jumps Du_ n and Du_ t coincide with either Du_ x or Du_ y , depending on the considered active interface. A custom-built algorithm is implemented in the Matlab script aiming at identifying which type of interface is considered, and then assigns the correct expression for both Du_ n and Du_ t . Afterwards, Eqs. (10) and (11) are set equal to Eqs. (4) and (5) or (6) and (7). For instance, considering a single interface I between elements i and j and normal to axis X ,their combinations become: j i u_ x;per  u_ x;per 

2 X

    q k_ I Aqn þ E_ xx xjG  xiG þ v_ xx zjG xjG  ziG xiG ¼ 0

q¼1

ð12Þ

ð6Þ j i u_ y;per  u_ y;per 

    þ E_ xy ylG  ykG þ 0:5v_ xy zlG ylG  zkG ykG

2 X

2 X

    q k_ I Aqt þ E_ xy xjG  xiG þ 0:5v_ xy zjG xjG  ziG xiG ¼ 0

q¼1

ð7Þ

ð13Þ

As shown in [30], a kinematically admissible velocity field must satisfy constraints given by an associated flow rule. For instance, if a Mohr-Coulomb failure criterion is employed, its related associated flow rule is expressed in terms of normal stress rn and tangential stress s as:

Using a matrix formulation to compact Eqs. (12) and (13), these become:

jsj  c  rn tan /

ð8Þ

Two straight lines are needed to describe the bounding yield surface for the Mohr-Coulomb criterion (Fig. 3), and their expressions are linear both in s and rn :





s þ rn tan /  c ¼0 F ðs; rn Þ ¼ s þ rn tan /  c

ð9Þ

in which / is the friction angle and c the cohesion of the active interfaces. It must be noted that the bounding yield surface must always be expressed with the general form Aqn rn þ Aqt s  C qI ¼ 0 regardless of the chosen failure criterion.

2



1 1 0 0

2

j u_ x;per

6 6 _ i 6 ux;per 6 j 1 1 6 6 u_ y;per 4 i u_ y;per 0 0

 xjG  xiG 0 0 6   þ4 0 0 xjG  xiG

3

7 2 3 7

 1 7  tan /  tan / k_ I 7þ 4 5 7 2 1 1 7 k_ I 5

3 E_ xx 7 6 36 E_ yy 7   7 6 6 E_ 7 0  zjG xjG  ziG xiG 0 0 76 xy 7   56 7¼ 6 v_ xx 7 0 0 0 0:5 zjG xjG  ziG xiG 7 6 6 v_ 7 4 yy 5 v_ xy 2

ð14Þ

4

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

This matrix formulation can be written even more compactly:

_ ij Aeq;ij 11 uper

þ

_ Aeq;ij 13 kI

þ

_ Aeq;ij 14 D

¼0

ð15Þ

From Eq. (15) it is clear that the periodic velocity field of the ij

two adjoining elements i and j (collected in the vector u_ per ), the plastic multiplier rates of interface I (k_ I ), and the components of the average strain and curvature rate tensors (both collected in _ are among the unknown variables of the linear prothe vector D) gramming problem. It must be noted that these are not the only unknown variables of the overall problem, as some more are introduced in the following section. After some simple assemblage operations, the constraint in its global form becomes: eq _ eq _ _ Aeq 11 uper þ A13 kI;ass þ A14 D ¼ 0

ð16Þ

2.3. Master-slave relations for unit elements The kinematic field of the finite elements that pertain to a masonry unit is governed by master–slave relations that link the kinematic field of a single finite element and that of the masonry unit to which it belongs (see Fig. 4). This is enforced in order to reduce the computational effort needed for solving the linear programming problem; the kinematics of unit finite elements is expressed in such a way that their interfaces are not considered in the overall solution, mirroring what happens in real case studies where failure within masonry units is unlikely to occur [23]. The master-slave relations are expressed according to Eqs. (17)–(19):

the X and Y axes of the master macroelement, respectively, so that its kinematics is enriched enabling macroscopic rotations that are representative of the deformed shape associated to an out-ofplane load condition. It must be remarked that the finite elements do not actually rotate, as previously stated in Section 2.1, but the general deformed shape of a masonry unit is indeed able to simulate rotations in a ‘‘pixeled” way. The expansions of Eqs. (17)–(19) lead to the following expressions:

  S   S  M S M S M M _ _ u_ x;per  u_ x;per þ h_ yy zM G  zG þ Exx xG  xG þ Exy yG  yG  S S    M S S M M _ þ v_ xx zG xG  zM G xG þ 0:5vxy zG yG  zG yG ¼ 0

ð20Þ

  S   S  M S M M M _ _ u_ y;per  u_ y;per þ h_ xx zSG  zM G þ Exy xG  xG þ Eyy yG  yG  S S    M S S M M _ þ v_ yy zG yG  zM G yG þ 0:5vxy zG xG  zG xG ¼ 0

ð21Þ

   M S M S M _M S u_ z;per  u_ z;per þ h_ xx yM G  yG þ hyy xG  xG h  2   2 i h  2   2 i þ 0:5v_ xx xM  xSG  ySG þ 0:5v_ yy yM G G   M S S þ 0:5v_ xy xM G y G  xG y G ¼ 0

ð22Þ

These can also be written in a matrix formulation, which is here omitted for sake of brevity. The global compact formulation of the constraints coming from the master–slave relations is: eq _ eq _ _ Aeq 21 uper þ A22 R þ A24 D ¼ 0

ð23Þ

 M S M u_ x ¼ u_ x þ h_ yy zSG  zM G

ð17Þ

Eq. (23) gives the remaining set of unknown variables of the linear programming problem, represented by vector R_ that includes the periodic velocity fields and macroscopic rotations of all the masonry units of the considered masonry panel.

 M S M u_ y ¼ u_ y  h_ xx zSG  zM G

ð18Þ

2.4. Periodicity boundary conditions

   M S M M _M S u_ z ¼ u_ z þ h_ xx ySG  yM G  hyy xG  xG

ð19Þ

M M M u_ x , u_ y , and u_ z are the components of the displacement rate field of a single masonry unit, which acts as the master macroelement, and M M they are expressed according to Eqs. (1)–(3). xM G , yG , and zG are the coordinates of the centroid of the considered masonry unit, which S S S are automatically calculated by the Matlab script. u_ x , u_ y , and u_ z

are the components of the displacement rate field of a generic finite element - in this case, the slave element – that belongs to the considered masonry unit; they are also expressed according to M M Eqs. (1)–(3). The quantities h_ and h_ represent the rotations about xx

The homogenization approach requires the inclusion of periodicity constraints applied to the velocity field at the boundaries of the investigated panel. These must be enforced on elements lying at the external boundaries of the panel that are orthogonal to axes X and Y; the elements here involved are located at the opposite sides of each of those boundaries (Fig. 5). a b a b u_ x;per ¼ u_ x;per ) u_ x;per  u_ x;per ¼ 0

ð24Þ

a b a b u_ y;per ¼ u_ y;per ) u_ y;per  u_ y;per ¼ 0

ð25Þ

yy

Fig. 4. Graphical representation of a masonry unit (M) and a unit finite element (S).

Fig. 5. Periodicity boundary conditions.

5

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206 a b a b u_ z;per ¼ u_ z;per ) u_ z;per  u_ z;per ¼ 0

ð26Þ

Z Pint ¼ A

The same applies for the generic couple of elements c and d. The global compact formulation of this constrain is:

_ Aeq 31 uper

¼0

ð27Þ

Two sets of out-of-plane load conditions are defined through two separate ‘‘loading angles”, as shown in Fig. 6. Angle wflex is named ‘‘flexural loading angle”: it represents the arctangent of the ratio between the macroscopic horizontal and vertical bending moments M yy and M xx , and varies between 0 and 2p. Similarly, angle wtors is named ‘‘torsional loading angle”: it represents the arctangent of the ratio between the macroscopic torsional moment and the vertical bending moments Mxy and Mxx , again ranging between 0 and 2p. Also, an in-plane load condition is included and defined in the same way as shown in [23], where the reader is referred to for a complete description. Overall, the dissipated external power is expressed as the summation of the products between the macroscopic stresses and bending moments and their associated average strain and curvature rate tensors components. Eventually, the dissipated external power is normalized and set equal to 1 to enforce a restriction in terms of all the possible collapse mechanisms that are associated to the collapse load, which is unique:

Pext ¼ Rxx E_ xx þ Ryy E_ yy þ Rxy E_ xy þ Mxx v_ xx þ M yy v_ yy þ M xy v_ xy ð28Þ

In the end, Eq. (28) represents another constraint for the linear programming problem, concisely expressed in the form:

_ Aeq 44 D

¼1

ð30Þ

Substituting Eqs. (10) and (11), this becomes:

Z A

rn

2 X

q k_ I Aqn þ s

q¼1

2 X

! Z X 2  q q k_ I Aqt dA ¼ k_ I rn Aqn þ sAqt dA

ð31Þ

A q¼1

q¼1

The dissipated internal power can then be written as:

2.5. Normalization of dissipated external power

¼1

ðrn Du_ n þ sDu_ t ÞdA

ð29Þ

2.6. Power dissipation in velocity discontinuities The expression for power dissipated across a velocity discontinuity whose area is A can be written as:

Pint ¼ A

2 X

q k_ I C qI

ð32Þ

q¼1

For the global problem the assembled final expression in matrix form is:

Pint ¼ C TI;ass k_ I;ass

ð33Þ

2.7. Assembly and solution of the linear programming problem The homogenized limit analysis problem is eventually formulated as a linear programming problem. The objective function that ought to be minimized is the dissipated internal power as expressed in Eq. (33), and the equality constraints are given by Eqs. (16), (23), (27), and (29). For a smoother solution, this computational problem is formulated in a standard form, which requires that each unknown variable of the problem is greater or equal to zero. In particular, for this specific problem the following unknowns must be modified: the elements’ periodic velocity field, the masonry units’ periodic velocity field and macroscopic rotations, and the components of average strain and curvature rate tensors. All of them are expressed as the difference of two nonnegative quantities to satisfy the requirements of the standard form: S S;þ S; u_ i;per ¼ u_ i;per  u_ i;per

i ¼ x; y

ð34Þ

M M;þ M; u_ i;per ¼ u_ i;per  u_ i;per

i ¼ x; y

ð35aÞ

M M;þ M; h_ ii ¼ h_ ii  h_ ii

i ¼ x; y

ð35bÞ

þ  E_ ij ¼ E_ ij  E_ ij

i; j ¼ x; y

ð36aÞ

Fig. 6. Angles defining the load condition.

6

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

v_ ij ¼ v_ þij  v_ ij

i; j ¼ x; y

ð36bÞ

The standard form for this linear programming problem is:

Minimize C T X

ð37Þ

Subject to AX ¼ B

ð38Þ

X0

ð39Þ

where

2 6 6 A¼6 4

Aeq 11

Aeq 21 Aeq 31

Aeq 11

6 6 6 6 6 6 X¼6 6 6 6 6 6 4

Aeq 14

0

Aeq 22

Aeq 22

0

0

0 0

0

0

0

0 Aeq 44

0 2

Aeq 13

0

Aeq 21 Aeq 31

Aeq 24 0

Aeq 14 Aeq 24 0

3 7 7 7 5

ð40Þ

Aeq 44

3 þ u_ per 7  u_ per 7 7 7 R_ þ 7 7 7 R_  7 _kI;ass 7 7 7 7 D_ þ 5

ð41Þ

3. Validation of the proposed approach

D_ 

3 0 607 6 7 B¼6 7 405 2

ð42Þ

1 2

0

6 0 6 6 6 0 6 C¼6 6 0 6 C I;ass 6 6 4 0

3 7 7 7 7 7 7 7 7 7 7 5

ð43Þ

0 2.8. Construction of the homogenized failure surfaces The Matlab script includes a final part where the construction of homogenized failure surfaces for the chosen panel is actually performed. In the framework of the upper bound theorem of limit analysis, the kinematic limit multiplier v is the minimum among those computed for each kinematically admissible collapse mechanism, and may be expressed as the ratio between the dissipated internal and external powers:

lk ¼

The flexural and torsional out-of-plane homogenized failure surfaces are constructed as follows: after selecting a value for the loading angles wflex (or wtors ), the linear programming problem is solved and the kinematic limit multiplier v for that specific outof-plane load condition is obtained. The collapse bending (or torsional) moments are then calculated by multiplying v to the initial macroscopic moments ½ Mxx M yy  (or ½ M xx M xy  for the torsional failure surface). Globally, 41 different values of wflex (or wtors ) are investigated, ranging from 0 to 2p with a sampling step of p=20. Each pair of collapse moments vM xx and vM yy (or vM xy ) represents a pair of coordinates in the resulting out-of-plane homogenized failure surface, which is then piecewise linear and consists of the segments linking each adjacent pair of collapse moments. An analogous procedure is employed for extracting the in-plane homogenized failure surfaces. The post-processing phase also enables to plot the failure mode for a specific in- and out-of-plane load condition with the Matlab command patch: the position of each finite element is reconstructed using the solution values of the remaining variables of the linear programming problem (namely, those _ and D) _ and using Eqs. (1)–(3). included in vectors u_ per , R,

In this section, the proposed approach is validated by drawing comparisons with three distinct models used in previous works to describe the in- and out-of-plane behavior of masonry at collapse. In particular, the results concerning homogenized failure surfaces are critically discussed to assess the reliability of the present approach in correctly representing the in- and out-of-plane collapse behavior of common masonry bonds. In the first example, discussed in Section 3.1, a running bond masonry test-window is investigated in terms of homogenized in- and out-of-plane failure surfaces. The results are then compared to those obtained by Milani and Taliercio in two distinct works, one devoted to the in-plane behavior [19], the other to the out-of-plane case [28]. As previously mentioned in Section 1, they use the so-called ‘‘Method of Cells” as a homogenizing approach, which is typical of problem involving fiber-reinforced composites. Relevant failure modes in bending and torsion are also presented and discussed. In the second example, studied in Section 3.2, an English bond masonry test-window is investigated only in terms of homogenized out-of-plane failure surfaces. The results are then compared to those obtained by Cecchi and Milani [27]; as already discussed in Section 1, the authors do not use a rigorous homogenization approach, instead they enforce the equivalence between a 3D system of blocks connected by interfaces and a 2D Reissner-Mindlin plate through a kinematic identification. Also in this case, relevant failure modes in bending and torsion are presented and discussed. 3.1. Running bond masonry REV

Pint Pext

ð44Þ

Introducing the normalization of the dissipated external power of Eq. (29) implies that the denominator is equal to 1, which means that:

lk ¼ Pint ¼ C TI;ass k_ I;ass 



ð45Þ 

v ¼ min lk ¼ min C TI;ass k_ I;ass



ð46Þ

Hence, the kinematic load multiplier is directly determined by the minimization of the dissipated internal power, which is correctly the objective of the linear programming problem here developed.

The first case study deals with a running bond masonry REV that consists of standard Italian bricks (25  12  5.5 cm3) and mortar joints whose thickness is 1 cm; the in-plane layout is shown in Fig. 7a. Two meshes are used for validating the proposed model: one is extremely refined, consisting of about 260,000 elements and displaying 48 elements over the thickness (Mesh 1, Fig. 7b), the other is extremely coarse, consisting of about 17,000 elements and displaying 12 elements over the thickness (Mesh 2, Fig. 7c). 3.1.1. In-plane homogenized failure surfaces For the validation in terms of in-plane homogenized failure surfaces, a Mohr-Coulomb failure criterion with cut-offs in tension

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

7

Fig. 7. (a) In-plane layout of the running bond masonry REV; (b) Mesh 1; (c) Mesh 2.

and compression is used. The mechanical properties of the interfaces are the same as those used in [19] and are listed in Table 1, with the tensile and compressive strength f t and f c evaluated according to the formulas there reported. Since Milani and Taliercio consider mortar as being governed by a Mohr-Coulomb failure criterion with no further comments on possible differences regarding the unit-mortar interfaces, these values are employed for both unit-mortar and mortar-mortar interfaces. It must be remarked that the unit-unit interfaces are not considered since the unit elements obey the previously introduced master-slave relations. Fig. 9 shows the comparison between the in-plane homogenized failure surfaces obtained by Milani and Taliercio [19] and those resulting from the present model for Mesh 1; in particular, Fig. 9a–c shows the failure surfaces in the tension-tension range for three different values of angle h (which represents the inclination of the principal directions with respect to the reference system of the REV, see Fig. 8), whereas Fig. 9d–f shows the failure surfaces in the compression-compression range.

Table 1 Mechanical properties for the material employed in comparison with in-plane MilaniTaliercio. Cohesion c [MPa]

Friction angle / [°]

Tensile strength f t [MPa]

Compressive strength f c [MPa]

0.1

36

2c cos / 1þsin /

2c cos / 1sin /

Fig. 8. Graphical description of angle h.

As far as the tension-tension range is concerned, the results obtained from the present model are thoroughly consistent with those obtained by Milani and Taliercio, despite some minor differences for the case of h equal to 45°. Conversely, the results in the compression-compression range show some slightly more marked discrepancies, which are rather evident for the case with h equal to 22.5°: these may be due to the different homogenization approach employed in the present model. It must be noted that, regardless of the value assumed by h, the shape of the failure surfaces in the compression-compression range is more similar to those obtained by Milani and co-workers [16], albeit for different material properties. In fact, that work considers mortar joints reduced to interfaces, which more closely resembles the assumption of the present model (i.e. dissipation only occurs at the interfaces between adjacent elements). Nonetheless, the homogenized collapse stresses for uniaxial macroscopic load conditions (namely, uniaxial tension and compression, vertical and horizontal) are almost always coincident for any value of h. Four failure modes are also shown in Fig. 10 for relevant load conditions applied to Mesh 1. As it can be seen, splitting occurs in both bed and head joints for uniaxial horizontal tension, whereas for uniaxial vertical tension it only takes place in bed joints; crushing occurs for both uniaxial horizontal and vertical compression, with the former also displaying some splitting in bed joints. Therefore, these results are consistent with the expected deformed configurations at collapse for the four considered cases. 3.1.2. Out-of-plane homogenized failure surfaces For the validation in terms of out-of-plane homogenized failure surfaces, a simple Mohr-Coulomb failure criterion is used; the values of cohesion c and friction angle / are 0.132 MPa and 27°, respectively, to be consistent with the values used in [28]. Fig. 11 shows the comparison between the out-of-plane homogenized failure surfaces obtained by Milani and Taliercio and those resulting from the present model for both Mesh 1 and Mesh 2; the latter is included to investigate the influence of the mesh size on the results. The results coming from the present model display a satisfying correspondence to those obtained by Milani and Taliercio; only the shape of the out-of-plane homogenized failure surface in the Mxx  M xy plane is slightly different, but this may depend on the different kinematics of the two models. The use of coarser Mesh 2 does not affect the results in a significant way, producing only a small reduction of the out-of-plane failure surfaces. Three failure modes are also shown in Fig. 12 for relevant out-of-plane load conditions applied to Mesh 1. As expected, the failure mode coming from the application of vertical overturning bending moment M xx presents vertical cracks in the head joints, whereas the failure mode coming from the

8

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

Fig. 9. Comparison of in-plane homogenized failure surfaces: (a–c) tension-tension range; (d–f) compression-compression range.

application of horizontal overturning bending moment M yy presents horizontal cracks across the bed joints. The failure mode coming from the application of the twisting moment Mxy is also consistent with the expectations, displaying cracks corresponding to a clear torsional deformed shape.

3.2. English bond masonry REV The second case study deals with an English bond masonry REV that also consists of standard Italian bricks (25  12  5.5 cm3) and mortar joints whose thickness is 1 cm. Its finite element mesh

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

9

Fig. 10. Failure modes for in-plane load conditions on Mesh 1: (a) uniaxial horizontal tension; (b) uniaxial vertical tension; (c) uniaxial horizontal compression; (d) uniaxial vertical compression.

Fig. 11. Comparison of out-of-plane homogenized failure surfaces with Milani-Taliercio: (a) Mxx  Myy plane; (b) Mxx  M xy plane.

counts about 35,000 elements (25 over the thickness), and is shown in Fig. 13; it must be noted that the chosen REV is not exactly the same as the cell used by Cecchi and Milani in [27], but it has been selected in order to ensure a rigorous application of homogenization. Also for this case study, a Mohr-Coulomb failure criterion is used; the mechanical properties for the interfaces are listed in Table 2 and are the same as those used in the work by Cecchi and Milani. Fig. 14 shows the comparison between the out-of-plane homogenized failure surfaces obtained by Cecchi and Milani and those resulting from the present model. For the M xx  M yy plane, a satisfying overall correspondence is noted, despite a difference in terms of inclination of the lateral boundaries of the failure surface; it is

highly possible that this discrepancy comes from the difference in terms of approach used for the evaluation of the out-of-plane failure surfaces between the two models. For the M xx  M xy plane, the failure surface coming from this model is smaller than the one obtained by Cecchi and Milani; this can also be due to the different approaches used in the two applications. In fact, the work by the aforementioned authors does not consider the actual thickness of mortar joints; instead, since their cell just consists of masonry units, dissipation can only occur across the surfaces shared by adjoining units. Furthermore, their approach also enables dissipation to occur along the out-of-plane direction: this has an impact on the results obtained when the macroscopic torsional moment Mxy is significant, since in this case the dissipation between the

10

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

Fig. 12. Failure modes for out-of-plane load conditions on Mesh 1: (a) M xx ; (b) Myy ; (c) Mxy .

Fig. 13. English bond masonry REV used for validation against data from CecchiMilani [27].

Table 2 Mechanical properties for the material employed in comparison with Cecchi-Milani. Cohesion c [MPa]

Friction angle / [°]

Tensile strength f t [MPa]

0.132

27

c tan /

two central blocks in correspondence of the mortar joint at the cell’s midplane becomes relevant. In order to accurately assess how the results in the M xx  M xy plane are affected by the differences between the two approaches, the upper bound limit analysis problem presented in this work is suitably modified to include dissipation along the out-of-plane direction for those interfaces that are located in correspondence of the midplane mortar joint. Moreover, a second finite element mesh is created in which the mortar joints are extremely thin, aiming at mirroring more closely the cell’s configuration of Cecchi and Milani (Fig. 15a): this consists of about 240,000 elements, 51 over the thickness. To account for the dissipation also in the out-ofplane direction, it is necessary to express the Mohr-Coulomb failure criterion with a 3D surface in the rn  s1  s2 plane, which in principle is a cone. However, since linear programming requires the bounding yield surfaces to have linear expressions in their variables, a piecewise linear approximation of the failure criterion is needed in the s1  s2 plane. Fig. 15b shows the very refined one used by Cecchi and Milani; to keep the swiftness of the present model, two rougher approximations are instead used in this

modified upper bound limit analysis problem, namely a boxshaped approximation of the surface (Model A, Fig. 15c) and one that circumscribes the actual surface (Model B, Fig. 15d). Fig. 16a and Fig. 16b show the out-of-plane homogenized failure surfaces in the M xx  M xy plane for Model A and Model B, respectively, compared against the results obtained by Cecchi and Milani and against those coming from the approach originally used in this work. It can be noted that the finer mesh offers more accurate results for both models, especially for Model B – even though it still slightly underestimates the boundaries of the homogenized failure surfaces in correspondence of M xy . To this regard, it must be recalled that the REV used in the present work is different than the cell used by Cecchi and Milani, which displays three pairs of central blocks instead of just one as in this case. Finally, it must be remarked that the out-of-plane homogenized failure surfaces in the Mxx  M yy plane do not show meaningful differences among the various cases and finite element meshes, hence they are here omitted for sake of conciseness. Eventually, three failure modes are extracted for the coarser mesh in relation to the initial model of Fig. 14, and they are shown in Fig. 17 for relevant out-of-plane load conditions. Similarly to the running bond masonry case, the three resulting deformed configurations for the selected out-of-plane load conditions are consistent with the expectations, thus confirming the reliability of the initial model. The failure modes for M xy related to Model A and B and for the coarser mesh are here omitted for sake of brevity, since they display no relevant differences to those shown in Fig. 17. 4. Conclusions and future developments In this paper, a full 3D model devoted to the derivation of inand out-of-plane homogenized failure surfaces for masonry has been presented, based on an upper bound limit analysis problem coupled with homogenization. This represents a major breakthrough in the field of masonry modeling, since at present no full 3D model based on limit analysis and homogenization is available in literature. A Kirchhoff-Love plate model has been employed for a complete description of the kinematic field of the considered masonry element, which must be discretized into rigid, regular brick finite elements. The model is then implemented into a Matlab script that consists of a linear programming problem in standard form, which enables the extraction of the sought in- and out-ofplane homogenized failure surfaces as well as of failure modes for the considered masonry element. Two different case studies are considered for validating the proposed model, one dealing with a running bond masonry REV, the other with an English bond masonry cell; the validation consists in comparing the results in terms of homogenized failure surfaces obtained from this model

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

11

Fig. 14. Comparison of out-of-plane homogenized failure surfaces with Cecchi-Milani: (a) Mxx  M yy plane; (b) M xx  M xy plane.

Fig. 15. (a) Finer mesh for the English bond masonry REV; (b) piecewise linear approximation for the 3D Mohr-Coulomb failure criterion used by Cecchi-Milani; (c) Model A, box-shaped approximation in the s1  s2 plane; (d) Model B, circumscribing approximation in the s1  s2 plane.

to those coming from previous investigations on the same case studies, which employ different methodologies. It is shown that there is a good match between the results of this model and those presented in the aforementioned previous works, expect for some discrepancies displayed by the results coming from the English

bond masonry cell under macroscopic torsional moment. However, these have been satisfactorily justified in light of more in-depth investigations concerning the differences between the present model and the one against which it is validated - namely, the absence of mortar joints and the possibility of dissipation along

12

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206

Fig. 16. Comparison of different out-of-plane homogenized failure surfaces in the M xx  M xy plane: (a) Model A; (b) Model B.

Fig. 17. Failure modes for out-of-plane load conditions on English bond masonry REV: (a) M xx ; (b) M yy ; (c) Mxy .

the out-of-plane direction that are included in the latter model. Also, some relevant failure modes are extracted for selected load conditions, and they are shown to be consistent with the expected outcomes. Future developments of this work will focus on the extraction of out-of-plane homogenized failure surfaces for non-periodic masonry panels and for multi-leaf masonry walls. These have been scarcely investigated in literature so far and represent applications of extreme interest since both are commonly featured in several historical buildings in Italy. In particular, multi-leaf walls usually display an insufficient resistance against out-of-plane actions such as those generated by a seismic shock. One upcoming work will also present a procedure to directly generate a suitable 3D FE mesh from the sketch of the considered non-periodic or multi-leaf masonry wall, transforming each pixel into a voxel (that is a 3D pixel) and then into a single finite element. Declaration of Competing Interest The authors declared that there is no conflict of interest. References [1] Pelà L, Cervera M, Oller S, Chiumenti M. A localized mapped damage model for orthotropic materials. Eng Fract Mech 2014;124–125:196–216. https://doi. org/10.1016/j.engfracmech.2014.04.027.

[2] Gatta C, Addessi D, Vestroni F. Static and dynamic nonlinear response of masonry walls. Int J Solids Struct 2018;155:291–303. https://doi.org/10.1016/ j.ijsolstr.2018.07.028. [3] Baraldi D, Cecchi A, Tralli A. Continuous and discrete models for masonry like material: a critical comparative study. Eur J Mech A-Solid 2015;50:39–58. https://doi.org/10.1016/j.euromechsol.2014.10.007. [4] Sarhosis V, Lemos JV. A detailed micro-modelling approach for the structural analysis of masonry assemblages. Comput Struct 2018;206:66–81. https://doi. org/10.1016/j.compstruc.2018.06.003. [5] Baraldi D, Reccia E, Cecchi A. In plane loaded masonry walls: DEM and FEM/ DEM models. A critical review. Meccanica 2018;53(7):1613–28. https://doi. org/10.1007/s11012-017-0704-3. [6] Addessi D, Sacco E. A 2D finite element based on an enriched kinematics for nonlinear analysis of masonry walls. Int J Mason Res Innov 2019;4(1– 2):97–112. https://doi.org/10.1504/IJMRI.2019.096819. [7] Macorini L, Izzuddin BA. A non-linear interface element for 3D mesoscale analysis of brick-masonry structures. Int J Numer Meth Eng 2011;85 (12):1584–608. https://doi.org/10.1002/nme.3046. [8] Zeman J, Novák J, Šejnoha M, Šejnoha J. Pragmatic multi-scale and multiphysics analysis of Charles Bridge in Prague. Eng Struct 2008;30(11):3365–76. https://doi.org/10.1016/j.engstruct.2008.05.012. [9] Di Nino S, Luongo A. A simple homogenized orthotropic model for in-plane analysis of regular masonry walls. Int J Solids Struct 2019;167:156–69. https:// doi.org/10.1016/j.ijsolstr.2019.03.013. [10] Milani G, Bertolesi E. Quasi-analytical homogenization approach for the nonlinear analysis of in-plane loaded masonry panels. Constr Build Mater 2017;146:723–43. https://doi.org/10.1016/j.conbuildmat.2017.04.008. [11] Milani G, Bruggi M. Simple homogenization-topology optimization approach for the pushover analysis of masonry walls. Int J Archit Herit 2018;12 (3):395–408. https://doi.org/10.1080/15583058.2017.1323248. [12] Reccia E, Leonetti L, Trovalusci P, Cecchi A. A multiscale/multidomain model for the failure analysis of masonry walls: a validation with a combined FEM/ DEM approach. Int J Multiscale Com 2018;16(4):325–43. https://doi.org/ 10.1615/IntJMultCompEng. 2018026988.

S. Tiberti, G. Milani / Computers and Structures 231 (2020) 106206 [13] Zucchini A, Lourenço PB. A coupled homogenisation-damage model for masonry cracking. Comput Struct 2004;82(11–12):917–29. https://doi.org/ 10.1016/j.compstruc.2004.02.020. [14] Rekik A, Lebon F. Homogenization methods for interface modelling in damaged masonry. Adv Eng Softw 2012;46(1):35–42. https://doi.org/ 10.1016/j.advengsoft.2010.09.009. [15] Berke PZ, Peerlings RHJ, Massart TJ, Geers MGD. A homogenization-based quasi-discrete method for the fracture of heterogeneous materials. Comput Mech 2014;53(5):909–23. https://doi.org/10.1007/s00466-0130939-3. [16] Milani G, Lourenço PB, Tralli A. Homogenised limit analysis of masonry walls. Part I: Failure surfaces. Comput Struct 2006;84(3–4):166–80. https://doi.org/ 10.1016/j.compstruc.2005.09.005. [17] Milani G, Lourenço PB, Tralli A. Homogenised limit analysis of masonry walls. Part II: Structural examples. Comput Struct 2006;84(3–4):181–95. https://doi. org/10.1016/j.compstruc.2005.09.004. [18] Milani G. Simple homogenization model for the non-linear analysis of in-plane loaded masonry walls. Comput Struct 2011;89(17–18):1586–601. https://doi. org/10.1016/j.compstruc.2011.05.004. [19] Milani G, Taliercio A. In-plane failure surfaces for masonry with joints of finite thickness estimated by a Method of Cells-type approach. Comput Struct 2015;150:34–51. https://doi.org/10.1016/j.compstruc.2014.12.007. [20] Stefanou I, Sab K, Heck J-V. Three dimensional homogenization of masonry structures with building blocks of finite strength: a closed form strength domain. Int J Solids Struct 2015;54:258–70. https://doi.org/10.1016/j. ijsolstr.2014.10.007. [21] Godio M, Stefanou I, Sab K, Sulem J, Sakji S. A limit analysis approach based on Cosserat continuum for the evaluation of the in-plane strength of discrete

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29] [30]

13

media: Application to masonry. Eur J Mech A-Solid 2017;66:168–92. https:// doi.org/10.1016/j.euromechsol.2017.06.011. Milani G, Lourenço PB. Monte Carlo homogenized limit analysis model for randomly assembled blocks in-plane loaded. Comput Mech 2010;46 (6):827–49. https://doi.org/10.1007/s00466-010-0514-0. Tiberti S, Milani G. 2D pixel homogenized limit analysis of non-periodic masonry walls. Comput Struct 2019;219:16–57. https://doi.org/10.1016/ j.compstruc.2019.04.002. Milani G, Lourenço PB, Tralli A. Homogenization approach for the limit analysis of out-of-plane loaded masonry walls. J Struct Eng 2006;132(10):1650–63. https://doi.org/10.1061/(ASCE)0733-9445(2006)132:10(1650). Sab K, Dallot J, Cecchi A. Determination of the overall yield strength domain of out-of-plane loaded brick masonry. Int J Multiscale Com 2007;5(2):83–92. https://doi.org/10.1615/IntJMultCompEng.v5.i2.20. Cecchi A, Milani G, Tralli A. A Reissner-Mindlin limit analysis model for out-ofplane loaded running bond masonry walls. Int J Solids Struct 2007;44 (5):1438–60. https://doi.org/10.1016/j.ijsolstr.2006.06.033. Cecchi A, Milani G. A kinematic FE limit analysis model for thick English bond masonry walls. Int J Solids Struct 2008;45(5):1302–31. https://doi.org/ 10.1016/j.ijsolstr.2007.09.019. Milani G, Taliercio A. Limit analysis of transversally loaded masonry walls using an innovative macroscopic strength criterion. Int J Solids Struct 2016;81:274–93. https://doi.org/10.1016/j.ijsolstr.2015.12.004. MATLAB Release 2018b. Natick (Massachusetts, United States): The MathWorks, Inc.. Sloan SW, Kleeman PW. Upper bound limit analysis using discontinuous velocity fields. Comput Methods Appl Mech Eng 1995;127(1–4):293–314. https://doi.org/10.1016/0045-7825(95)00868-1.