CH of masonry materials via meshless meso-modeling

In the present study a multi-scale computational strategy for the analysis of masonry structures is presented. The structural macroscopic behaviour is obtained making use of the Computational Homogenization (CH) technique based on the solution of the boundary value problem (BVP) of a detailed Unit Cell (UC) chosen at the meso-scale and representative of the heterogeneous material. The smallest UC is composed by a brick and half of its surrounding joints, the former assumed to behave elastically while the latter considered with an elastoplastic softening response. The governing equations at the macroscopic level are formulated in the framework of finite element method while the Meshless Method (MM) is adopted to solve the BVP at the mesoscopic level. The work focuses on the BVP solution. The consistent tangent stiffness matrix at a macroscopic quadrature point is evaluated on the base of BVP results for the UC together with a localisation procedure. Validation of the MM procedure at the meso-scale level is demonstrated by numerical examples that show the results of the BVP for the simple cases of normal and shear loading of the UC.


INTRODUCTION
odeling of structures constituted by heterogeneous materials, as masonry and composite laminates, is still one of the most complex and challenging matters of the modern structural engineering.Anisotropy, asymmetry, geometrical and material non-linearities are the main features of the mechanical behavior to be considered in the formulation of the structural model.In masonry structures the most relevant kinematical and mechanical phenomena take place at a scale which is small if compared to the dimensions of the structure.On the other side, the structure is governed, in its peculiar overall response, by its global geometrical and morphological configuration.In literature, two different scales of interest are distinguished, directly linked to as many theoretical approaches: the mesoscopic approach and the macroscopic approach.The mesoscopic approach considers units and their interfaces individually.Most of the efforts in this research area are concentrated on the formulation of advanced interface constitutive laws capable to describe the damage evolution and the onset of irreversible strains and related coupled effects for different mechanical problems.In Giambanco et al. [1] the masonry material is considered as an assembly of units in contact by means of elastoplastic zero-thickness interfaces, focusing on the cohesive-frictional joint transition.Lourenço and Rotz [2] proposed a joint failure criterion which adopts a cap model to take into account the mortar compaction due to high compressive stresses.Alfano and Sacco [3] introduced a damage-plasticity interface model that discriminated a linear elastic undamaged zone from a damaged zone with an unilateral Coulomb friction law.More recently, Spada et al. [4] and Sacco and Lebon [5] concentrated in the formulation of interface laws which include the non-linear effects of stiffness M degradation, loss of adhesion and unilateral contact with related frictional phenomena.The contribution of internal stresses and strains at the interface level has been finally inserted by Giambanco et al. [6] to enrich the interface model.The mesoscopic approach is a rigorous method of analysis but many difficulties arise in the mesh creation and a fine discretization of the structure has to be used, which lead to prohibitive computational costs for large structures.The macroscopic approach considers the structure as constituted by a fictitious homogeneous and continuous material.In this approach simplified constitutive models are formulated in a phenomenological manner, intentionally neglecting the actual material composition but smearing the effects of the heterogeneities presence by appropriate mechanical assumptions [7].The multiscale techniques belong to the second approach and couple different scales of interest by means of specific transition laws capable to exchange information between different consecutive scales [8][9][10].According to [11] three homogenization techniques exist: the analytical homogenization, the numerical homogenization, and the computational homogenization (CH).In particular the CH consists in the following steps.A macroscopic strain or stress is used to apply kinematic or static boundary conditions to the UC (macro-meso scale transition).The equilibrium of the UC is obtained by solving a BVP at the mesoscale.Lastly, the macro-strains or macro-stresses are assumed to be the average on the UC of the corresponding meso-strains or meso-stresses (meso-macro scale transition).In a strain driven process, the aforementioned boundary conditions are of Dirichlet type and are chosen linear or periodic.The latter imply the repetition of inelastic phenomena for adjacent UCs.Despite it is known in literature that periodic conditions offer a better estimate of stiffness matrix with respect to linear conditions, in this paper linear conditions are used, with the aim to verify the difference in future works.Two different CH methods are distinguished.In the first order CH method Cauchy models are used at all scales.This method is based on the satisfaction of the Principle of Separation of Scales, which asserts that the characteristic size of the UC is much smaller then the size of the structure [12].Massart [13] proposed an enhanced multiscale model for masonry structures using nonlocal implicit gradient isotropic damage models for both constituents of the UC.A second CH method employs Cosserat or higher order continua to consider those cases in which strong strain or stress gradients arise at the macroscopic level or the mesoscopic dimensions of the UC constituents are comparable with the ones at the macro-level.Addessi and Sacco [14] analyzed masonry panels in the framework of transformation field analysis by using a two-dimensional Cosserat continuum for macro-scale and a non-linear damage contact-friction model for the mortar joints at the meso-scale.Many authors have also faced the problem of localization of deformations.For ductile materials with an elastic-plastic behavior sufficiently deformed into the plastic range, the abrupt changes in the deformation field occur across narrow zones (shear bands).This phenomenon is called localization of plastic deformation.Rice [17] and Ottosen and Runesson [18] studied the problem and found the conditions under which the localization occurs.Massart et al [16] inserted the localization procedure into a multi-scale problem, considering a gradient damage model at the mesoscale.In this paper a multiscale computational strategy for the analysis of masonry structures is presented.The applied technique belongs to the first order CH methods.The UC chosen at the meso-scale is composed by a brick and half of its surrounding joints.Brick is assumed elastic while mortar joints are modeled by means of zero-thickness elasto-plastic interface elements.The BVP is solved applying the Meshless Method (MM) and the UC is subjected to linear boundary displacements.Unlike FEM, in MM the shape functions are only approximants and not interpolants.These are obtained by defining the influence of nodes on a fixed point by means of weight functions, which have a compact support.The support size is dependent by the so-called dilatation parameter or smoothing length.The smoothing length is critical for the solution accuracy and stability, it plays the role of the element size in the finite element method.The MM is employed to overcome the problem of computational cost associated with the use of a classical FEM analysis at the meso-scale.This method is also able to easily treat the problem of crack formation and propagation.It will therefore be applied to describe crack formation and propagation in the brick, that is one of the future developments of the present work.Another advantage is related to the possibility to use higher-order continuous weight functions.The macro-scale consistent tangent stiffness matrix and localization of deformation are derived.Numerical examples show the behavior of the UC and results of localization for two different cases, i.e. when boundary displacements enslave the UC to normal and shear strains.

THE COMPUTATIONAL HOMOGENIZATION (CH) SCHEME
et us consider, in the Euclidean space  3 referred to the orthonormal frame 1 2 3 ( , , , ) O i i i , a structure  constituted by an heterogeneous material, at the mesoscale (Fig. 1).L Figure 1: Mechanical scheme of a structure constituted by an heterogeneous material.
Our aim is to derive the structure-properties, or the properties at the macroscopic level (M), based on the kinematical and mechanical phenomena occurring at the mesoscale level (m) of the heterogeneous material.The boundary of the structure is divided in two parts  u and  t where kinematic and loading conditions are specified, respectively.The external actions are the body force f for unit volume and the surface tractions t .The structural response at the macroscopic level is given by the displacement M u , strains M ε and stress M σ fields.The equilibrium of the structure can be assessed in a weak form making use of the principle of virtual displacement (PVD), thus where M  ε are the virtual strains satisfying the following kinematical compatibility conditions expressed under the small displacement hypothesis: with C the kinematical compatibility matrix.The above reported variational formulation of the structure equilibrium can be solved in an approximated way making use of the finite element method.The continuous structure is discretized in a finite number of subdomains or elements e  (Fig. 2).The equilibrium of the discretized structure requires the equilibrium of all the elements.It can be expressed in a weak form as follows: The displacement at any point within the element is approximated making use of appropriate interpolation functions: being N the matrix of interpolation functions and M  U the virtual displacements vector of the element nodes.
In the CH methods, the macroscopic response is evaluated as the average of the response of the UC.The scale transition is based on the Hill-Mandel principle [20,21] which establishes that the virtual work density at the macro-scale must be equal the volume average of the virtual work at the meso-scale: From Eq. ( 3) and ( 5) it comes out that the equilibrium of the discretized structure constituted by the heterogeneous material can be solved if the equilibrium of the UC is assessed under specific boundary conditions.The downscaling operation, also known as macro-meso transition, consists in the evaluation of the boundary conditions to be applied to the UC.In the first-order computational homogenization the macro-meso transition is usually "deformation driven" because this procedure can be directly fit into a displacement-based finite element framework.Therefore, the macroscopic deformation tensor is calculated at the material point of the finite element and it is next used to formulate the kinematical boundary conditions to be imposed on the UC associated to this point.In literature the periodic boundary conditions are suggested for periodic materials.In the present work, considering that the UC mechanical response in the non linear stage is strongly affected by the decohesion process occurring at the material interfaces, firstly the linear boundary conditions are imposed to the UC, thus on where, for 2D applications:

D
is a matrix that depends on the x-y coordinates of the nodal point once a reference system is fixed.m u are the prescribed values of the displacements for the point of position x located on the boundary UC  of the UC.This choice could be easily modified on the basis of the obtained numerical results.Since the body forces can be neglected and the whole boundary is constrained, the variational form of the UC equilibrium can be expressed as follows: and it has to hold in presence of the kinematical constraint on The kinematical constraint (8) can be incorporated in the variational equality (7) making use of the Lagrange multiplier method, thus where r is the vector of Lagrangian multipliers and m  ε the virtual strain tensor of the UC related to the virtual displacements through the following compatibility condition: Imposing the Hill-Mandel principle and substituting ( 6) in ( 9), after some algebra the following equation is obtained: Since the previous equation has to hold for any value of  r and M  ε , it provides the following conditions: Eq. ( 12) is the macro-meso transition condition here adopted for the multiscale problem.Eq. ( 13) represents the equilibrium condition of the UC but also relates the macroscopic stress to the mesoscopic mechanical response in terms of reaction forces arising along the boundary of the cell and, in this sense, it can be considered as the meso-macro transition or upscaling equation.

UNIT CELL BOUNDARY VALUE PROBLEM
or the case of running bond masonry, according to [22] and [13], the smallest UC is constituted by a single block surrounded by half of the head and bed mortar joints and the periodicity is defined by the two directions 1 i and 2 i since the heterogeneous material can be constructed by repeating the UC along these directions (Fig. 3).In the present work the material of the blocks is considered indefinitely elastic and the interfaces constitutive laws, expressed in terms of contact tractions i σ and displacement discontinuities at the interface   i u , are developed in the framework of elastoplasticity for non standard materials.The following general assumptions are adopted to describe the mechanical behavior of the interfaces:  the traction components are continuous at the interface  the strain components along the thickness h of the interface are uniform and evaluated on the basis of the values assumed by the displacement discontinuity components  the total strain and the total displacement discontinuities are decomposed in the elastic (e) and inelastic (p) parts F With reference to the assumed hypotheses, the material of the block is considered elastic and inelastic displacement discontinuities arise at the zero-thickness interfaces, thus where b E and i E are the elastic matrices of the block material and of the interfaces, respectively.Irreversible discontinuous displacements occur when the interface stress state reaches a limit condition.The elastic domain is defined by two convex limit surfaces intersecting in a non-smooth fashion: the Coulomb bilinear limit surface and a tension cut-off (Fig. 4).The limit functions, reported in the stress space take the following form: where mi τ and mi  are the tangential stress vector and the normal stress component of the contact stresses,  is the friction angle, 0 c and 0  the cohesion and tensile strength of the virgin interfaces.
p  is a static variable which is associated to the internal variable p  which regulates the isotropic hardening-softening interface behavior: with p h the hardening-softening parameter.The inelastic displacement discontinuities develop according to a non-associative flow rule: where 1 p   and 2 p   are the plastic multipliers which satisfy the complementarity conditions 0, 0, 0, 0, 0 The plastic potential related to the limit condition (19) is expressed by the following function: with   0,    dilatancy angle and r an arbitrary material constant selected in such a way to satisfy Eq. ( 24) for any stress state.
The weak form of UC equilibrium (9) leads to the following conditions:

Meshless solution of the BVP
Generally the solution of the UC boundary value problem is approached in an approximated way making use of the finite element method and, since two finite element meshes are used at the macroscopic and mesoscopic levels, the method is known as FE 2 .The present study approaches the numerical solution of the mesoscopic model by means of a meshless strategy which is described in the following for plane stress conditions.In two-dimensional models of masonry structures the plane stress or plain strain assumption is adopted depending on the geometry of the structural element.A specific study on the range of validity of the two different assumptions has been developed by Anthoine [23].The principal conclusion is that in the linear range both provide satisfactory results while in the non linear range the plane stress assumption may lead to inaccurate results.In this case the plane stress assumption is adopted because the implementation is finalized to the study of masonry walls where the thickness is quite small if compared to the other dimensions.The plane stress condition is also used to develop numerical simulations on masonry specimens under the assumption that the non linearities are concentrated at the joint level and the block material is indefinitely elastic.In view of this last assumption inaccurate results related to the plane stress condition can be considered of the same order of the numerical approximations.
with   p x a monomial basis function and   a x a vector of coefficients which are functions of the spatial coordinates x .In the present case the basis function adopted is linear, , and the coefficients are obtained by performing the minimization of the following functional: where N is the number of nodes and x x is the weight function depending on the distance between the sampling point and the node j.Minimizing functional (30) the result is a system of linear equations where , , , Substituting the result (31) into (29) the MLS approximation function is finally obtained as where the shape function The weight function plays an important role in the performance of the MLS approximation.The quartic spline function is adopted in this work, where / j j r d   x x , being d the radius of the domain of definition of the point x (smoothing length).Making use of the MLS approximations for the two components of the displacement field, the kinematical compatibility conditions for the block and for the interfaces are where b Φ and i Φ are the approximation functions, b S and i S are the selectivity matrices, evaluated for the block and for the interface domains, respectively.C is the classical compatibility matrix for the plane stress case.In addition, the reaction forces on the boundaries k  are approximated in terms of nodal values R as where   Ψ x are the approximation functions for the reaction forces.Considering the kinematical Eq. ( 38), the approximation (39) and the elastic constitutive laws for the block (17) and for the mortar joints (18), the weak form of the UC equilibrium can be particularized in the following expression where i K is the stiffness matrix for interfaces, defined as Since Eq. ( 45) has to hold for any arbitrary m  U and  R , the governing equations of the UC are obtained, in the following matrix form: In order to evaluate the integrals (41)-( 44), the Gauss quadrature rule is adopted with nine sample points in the block domain and three sample points for each interface.Furthermore, it is assumed i i  Ψ Φ .In order to be introduced in a step-by-step algorithm, the governing Eq. ( 46) need to be rewritten in a discrete formulation.A Newton-Raphson procedure is employed at the mesoscale to find the solution of the BVP for a single load step.Each iteration of the load step is divided in an elastic trial predictor stage and a plastic corrector stage.

Consistent tangent stiffness matrix evaluation
In the time step It is clear how nodal displacements and reactions are dependent on two contributes that can be assumed equal to an elastic part (function of m U ) and a plastic part (function of p F ).Recalling ( 13), (39), (44) and (48) it is simply obtained that 22 21 1 1 The macroscopic tangent stiffness matrix is defined as: In the previous equation the second term has to be explicitly evaluated.After mathematical manipulation it is possible to write that where From Eq. (47) which, substituted in (51), furnishes: The consistent macroscopic tangent stiffness matrix is therefore written as the sum of an elastic and a plastic part:

Localization of deformation at a single macroscopic quadrature point
Let us consider, at a macroscopic quadrature point, a discontinuity band characterized by an incipient loss of strain continuity (Fig. 6).
Figure 6: Body split by a potential discontinuity band.
The band splits the volume in the two subdomains V  and V  and is identified by the unit vector n oriented towards V  .The traction continuity condition imposes that where  σ  and  σ  denote the rates of stress vectors on the negative and positive sides of the discontinuity band, respectively.n C is the kinematical compatibility matrix particularized for the surface having unit normal vector n .A kinematical condition has also to be imposed, concerning displacement continuity.This condition leads to the equation In (59)  ε  and  ε  are the rates of strain vectors on the negative and positive sides.The following position has also been used: with c representing the jump magnitude, while m is a unit vector called polarization vector.
The angle between m and n unit vectors characterizes the failure mode: when m is aligned with n a splitting mode I process is present; when m is perpendicular to n a shear mode II is present.In general, V  and V  can be constituted by different materials and the following constitutive relations can be written: with ct  E and ct  E the tangent stiffness matrices in V  and V  .Substituting the expressions ( 61) and ( 62) on (58), the following equation is obtained: If eq. ( 59) is also introduced we finally derive: In the case of continuous bifurcation (when the tangent stiffness matrices in V  and V  are the same: ), the condition of incipient macroscopic localization occurs if it happens that In the last condition is called localization matrix (or acoustic matrix).Since m is a unit vector, the condition of incipient localization can be written as Besides, localization occurs when the stiffness matrix loses its positive definiteness.This means that the set of orientations n for which localization may appear is determined by the inequality [16]: In a localization procedure n is not known in advance.Therefore, scanning for all possible unit vectors n , the discontinuous band is chosen looking for all n directions that comply with (68) and whose minimum eigenvalue of L matrix is the absolute minimum one.The corresponding eigenvector finally identifies the polarization vector.

NUMERICAL EXAMPLES
o show the effectiveness of the Meshless procedure for BVP solution and of the localization technique, the response of a UC under the two displacement conditions showed in Fig. 7 is analyzed.The UC is composed by a brick having length equal to 204 mm and height equal to 50 mm, surrounded by a mortar layer of thickness equal to 5 mm.The mechanical parameters adopted for the two materials are reported in Tab.

Example 1. Horizontal strain applied
In this case the external boundary nodes on 2  and 4  interfaces are forced to move along the 1 x direction.As a consequence, the inner block is stretched and results in a biaxial tensile stress condition.The UC average stress versus average strain diagram is reported in Fig. 8. Different stages can be distinguished.At the beginning the response is elastic.The slope of the curve changes when plasticity starts to appear on vertical interfaces, at around 0.8 MPa of stress level.At the end of this first nonlinear branch the plasticity level at each Gauss point of the interfaces is showed by histograms of Fig. 9a.An histogram ranges from 0 (absence of plasticity, white colour) to 1 (complete loss of cohesion-adhesion, dark colour).Increasing the strain, plasticity develops on horizontal interfaces also until the maximum stress is reached.A softening branch follows this stage until all Gauss points, unless those on the symmetry axis, have their maximum plasticity level (Fig. 9c,d).
T A residual stress is finally observed due to the UC confinement.The results of localization procedure is also showed at four stages of Fig. 9.The localization matrix has a positive determinant until the maximum stress is obtained and no localization can be identified before.Once the maximum load is reached det 0  L for 0 60 65        but the minimum eigenvalue is obtained for 0    (Fig. 10), that is the direction of n vector.Eigenvector analysis of L matrix furnishes an m vector aligned with n , that confirms the mode I failure of the UC.

Example 2. Shear strain applied
In this second example the UC is subjected to a monotonically increasing shear strain.In this case plasticity appears at two opposite corners of the UC due to shear stresses at each interface Gauss point.The average stress-strain curve is nonlinear and evolves from an elastic stage towards a residual line governed by Mohr-Coulomb limit condition and dilatancy angle (Fig. 11).It is important to highlight that after the localization procedure the n direction is close to the diagonal of the UC, while m evolves during the test and, if at the beginning it seems to be closer to n , at the end the two vectors are almost perpendicular.In this case the failure mode is of mixed type with a tendency to became of mode II (Fig. 12a-d).
The angles of n and m vectors with respect to the 1 x axis at the final configuration are of 15.9° and 86.4° respectively.
In Fig. 13 the determinant of localization matrix and minimum eigenvalue curves versus orientation  of n vector are reported.

CONCLUSIONS
his paper focuses on a first order CH method based on a FE-Meshless multiscale procedure to be employed to simulate running bond masonry materials in plane stress regime.The UC is chosen equal to a brick surrounded by half mortar joint layer and the boundary value problem is solved by means of a Meshless approach, which is the principal novelty of the paper.Linear displacements are applied on the UC boundary.The results will be compared with those obtained using periodic boundary conditions in a future work.A simple system represents the solving equations of T the BVP.Also, the consistent macroscopic tangent stiffness matrix is evaluated and a localization procedure to find discontinuity bands at the macroscale is taken into account.Two numerical examples focusing on the UC response only are reported.The results of the BVP are analyzed to confirm the goodness of the Meshless Method to be employed to solve the problem at the mesoscale.The localization procedure for the two examples correctly identifies the failure mode type.Future developments regards the extension of this strategy to other heterogeneous periodic or non-periodic materials also for 3D analyses and the introduction of a non linear behavior for bricks.Particular attention will be paid on meshdependency for the response at the macroscale level.

Figure 2 :
Figure 2: Mechanical scheme of the CH strategy.

Figure 3 :
Figure 3: Unit cell extracted from running bond masonry.

Figure 4 :
Figure 4: Bilinear plastic limit condition represented in the plane stress space.

Figure 5 :
Figure 5: Unit cell meshless model.With reference to Fig. 5 the UC is divided in five integration domains: the first ( b  ) corresponds to the volume occupied by the block.The others integration domains are the interfaces ( 1, ,4) k k    .Totally seventeen nodes are used to discretize the UC but different choices are possible in order to refine the numerical solution.The displacement field inside each sub-domain is obtained from the nodal displacement values m U by the Moving Least Square approximation (MLS).The approximated value of a displacement component is expressed as a polynomial function:

Figure 7 :
Figure 7: Boundary conditions imposed to the UC in the three examples analyzed.
1, while interfaces' parameters are reported in Tab. 2. The UC is constrained for its entire boundary surface.The tests are run under displacement control.The elastic tangent stiffness matrix is: