Send Orders of Reprints at Reprints@benthamscience.org Homogenization Approach for the Evaluation of Crack Patterns Induced by Foundation Settlement on an Old Masonry Building

A numerical FE homogenization approach for the interpretation of existing crack patterns induced by foundation settlement on old masonry buildings is discussed. The approach is quite general and may be applied to any case study. It relies on a 3D FE discretization of the entire structure by means of rigid infinitely resistant six-noded wedge elements and non-linear interfaces, exhibiting deterioration of the mechanical properties. Soil is modeled by means of elastic translational springs, with values derived from at hand simplified approaches. The case study analyzed is the so called " Siloteca " [1] building in Milan, Italy, which belongs to a more complex built aggregate , originally conceived by the French Napoleonic army as riding stable during the Cisalpina Republic. At present, the building is utilized as an archive within the Science and Technology Museum. The aggregate may be regarded as being subdivided into two separate blocks, with each one further subdivided into eight isolated buildings. Nowadays, only six stables of one of the blocks are still present. Two of the six structures now serve as a Museum deposit and are the objects of the present study, whereas the other four are in worst condition, partially roofed and collapsed and in a general decayed state. Siloteca exhibits meaningful crack patterns and an active overturning mechanism of the main façade. The masonry face texture is relatively regular and well organized; its section is constituted of three leaves of header bricks, one leaf being alternatively constituted by one-half brick. A quite large sub-vertical crack is present in the central long wall, at a distance of about 10 meters from the main façade, which is progressively rocking. The reason of the façade movements at the base is probably due to differential settlements of the foundation, as a consequence of a large excavation realized some decades ago to install large gas oil tanks for the museum. In this paper, for a direct mechanical interpretation of the reasons at the base of the formation of the crack pattern, a simple but effective fully equilibrated model is discussed, facilitating in the accurate prediction of the position of the cracks. The model is also utilized to estimate soil elastic vertical stiffness –within a Winker approximation-to be used in a second phase with the fully non-linear FE code. Once the soil constants are at disposal from such a procedure, a homogenized non-linear FE code recently proposed by the …


INTRODUCTION
Foundation settlement is a major cause of damage for existing masonry structures.In the most general cases, vertical movements in foundations are caused by downward movement of a wall or wall footings.Vertical movements may be different for different walls, due to the presence of different soft soils, excavations, mining subsidence, application of overloads to the structure, etc.
As a consequence of the very poor resistance of masonry in tension, a foundation settlement may result into the formation of meaningful crack patterns, up to the activation of a failure mechanism, associated with the collapse of part of the structure.
Apart from monitoring the evolution of the existing crack patterns, the use of numerical procedures appears quite interesting aimed at the prediction of the crack wide opening evolution.Such approach may become crucial when rehabilitation interventions are planned in order to interrupt cracks spreading or rehabilitate completely the structure.
To be fully predictive, numerical approaches should take into account both the non-linear behavior of masonry and the effect induced by the soil on the structure.When dealing with this second aspect, a 3D FE model for the soil, assumed eventually as a non linear material, should be used, with a so-called "per phases analysis" in the presence of excavations or mining subsidence problems.However, it is immediately clear that the adoption of such approach would result into prohibitive computational costs to perform even a single non-linear analysis for small and medium scale buildings.
The first aspect, i.e. masonry non-linearity, adds further complexity to the model.As a matter of fact, masonry is a composite material constituted by units, e.g.clay bricks and large natural blocks, and mortar joints.The variable dimension of units, joint width, the material properties of blocks and mortar, the arrangement of bed and head joints and the quality of workmanship, make the simulation of masonry structures extremely difficult.Moreover, an accurate masonry description needs a complete set of experimental data.Two main approaches have been developed to describe masonry behavior over the elastic limit -which usually is exceeded at very low levels of external loads-, usually known in the technical literature as macro-modeling and micromodeling.
Macro-modeling does not make any distinction between masonry units and joints, averaging the effect of mortar through the formulation of a fictitious continuous material.It has been widely used in the past (e.g.no-tension material [4]), because it allows the rough discretizations necessary for actual large scale structures.Nevertheless, it appears really difficult to take into account some distinctive aspects of masonry, such as anisotropy in the inelastic range and the postpeak softening behavior.In order to consider such aspects, some equivalent macro-models have also been presented [5][6][7], featuring orthotropic elastic-plastic behavior with softening.However, the mechanical properties required are derived from experimental data fitting.
The alternative micro-modeling consists of representing separately mortar joints and units.Usually, joints are reduced to interfaces [8][9][10] in order to limit the computation effort.Nevertheless, the need of modeling separately bricks and mortar limits its applicability to small panels.
At present, the analysis of masonry walls in the inelastic range still requires macro-scale computational approaches conducted, e.g., through Fes [11,12].Recently, efficient models based on homogenization have been presented [13][14][15][16][17][18][19][20], which allow non-linear analyses of large scale structures, still considering both the real disposition of bricks and the actual mechanical properties of the constituent materials.Clearly, the numerical models applied at a structural level should be sufficiently simple, reliable and efficient to allow a quick evaluation of (a) collapse loads, (b) displacements near collapse and (c) post peak behavior of the structures.
Homogenization consists of extracting a representative element of volume (REV) which generates the whole structure by repetition, in solving a boundary value problem on the REV and in substituting the assemblage of bricks and mortar at a structural level with a fictitious orthotropic equivalent material.The most straightforward procedure is the utilization of FEs [15], assuming either elasto-plastic or damaging constitutive laws for units and mortar.Nevertheless, FEM requires a great computational effort, since the field problem has to be solved numerically for each loading step, in all Gauss points.
To contemporarily take into account the actual masonry non-linear behavior at a structural level and foundation settlement in an approximate but effective way, with the application of a very limited computational effort, a simplified homogenization method is used to substitute the heterogeneous assemblage of blocks and mortar with a macroscopic equivalent material through a simplified averaging procedure.Soil is modeled within a Winkler approach, where the translational stiffness of the springs is evaluated by means of a simplified equilibrated procedure, where masonry walls are regarded as rigid bodies translating and rotating when subjected to vertical loads.
The FE approach is tested on a real scale masonry structure exhibiting critical crack patterns, induced by a founda-

Olivetani monastery
Siloteca building tion settlement.The building under consideration is called "Siloteca" and is represented in Fig. (1) and Fig. (2).It is a historical masonry building with quite regular shape located in Milan, Italy.Siloteca is part of the so called "Cavallerizze" built aggregate, which owes its name to its original destination, being conceived by the French Napoleonic army as riding stable during the Cisalpina Republic (beginning of 1800).At present, the building belongs to the Science and Technology Museum and is utilized as an archive.The aggregate may be regarded as being subdivided into two separate blocks, each one further sub-divided into eight isolated buildings.Nowadays, only six stables of one of the blocks are still present.Two of the six structures still present now serve as Museum deposit (Siloteca), whereas the other four are in worst condition, partially roofed and collapsed and in a general decayed state.
Within this research activity, a detailed photographic and graphical survey of the existent crack pattern has been conducted.From geometric and photographic surveys, shared with visual inspection, a globally uniform character of the structural system of the Cavallerizze complex has been ascertained.The masonry face texture is relatively regular and well organized; its section is constituted by three leaves of header bricks (thickness 60 cm, bricks dimensions approximately equal to Italian standard bricks, i.e. 250×120×55 mm 3 ), one leaf being alternatively constituted by one-half brick.The structural and architectural unitary of the complex, can be extended also to the part now demolished, originally placed between the two bodies, as demonstrated by the residual parts of the walls interconnected with the edge walls of the Cavallerizze.
Siloteca exhibits an important crack pattern, with a clear partial detachment of the façade with part of the perpendicular walls and its progressive rocking, probably due to differential settlements near the façade foundation, resulting from a large excavation realized some decades ago to install large gas oil tanks for the museum.
In this paper, a simple but effective fully equilibrated model is discussed, proving useful to interpret the causes at the base of the formation of the crack pattern, and presenting the ability to predict quite accurately the position of the cracks.The model is also utilized to estimate soil elastic constants -within a Winker approximation-to be used in a second phase within a fully non-linear FE code.The at hand model is constituted by a rigid block subjected to self-weight and dead loads actually present in the building, reactions forces of the soil and internal actions on the section subjected to cracking.Crack position is sidelined as unknown problem and is determined for checking the maximum state of stress present within masonry.The knowledge of the actual position of the cracks allows, by means of direct displacement identification with experimental data, an analytical evaluation of foundation stiffness.
Once soil elastic constants are at disposal from such a procedure, the overall Siloteca stable is considered and a homogenized non-linear FE code [2,3] is utilized to have an insight into the state of mechanical degradation of the structure.A hypothesis on foundation settlement is provided to justify crack maps exhibited by the structure.

DESCRIPTION OF THE MATERIALS AND STRUCTURAL DEGRADATION OF SILOTECA AND CAVALLERIZZE BUILDINGS
The crack pattern present in the Siloteca building (Figs. 3  and 4) is characterized by a quite meaningful sub-vertical crack on the central wall (front B and B'), combined by a series of cracks on the walls belonging to fronts E, E', E'', E'''.After a careful historical analysis of the construction, it has been deduced that, most probably, such crack pattern has been determined by a settlement of the foundation soil intersecting the zones near fronts E and E'''.In addition, a localized settlement intersecting the corner between fronts A' and E is visible.Differential settlement is consequent to an excavation realized near fronts E and E''' to host two large Diesel The failure mechanism on the building, at present almost active, is an overturning of the façade, with partial detachment and in-plane rotation of the perpendicular walls A, B and C. A wide crack is indeed present on central and lateral perpendicular walls, see Fig. (5), still active and now monitored experimentally to evaluate the damage evolution.
A schematic representation of the dimension of the excavated zones and the tanks, deduced from the documentation available is sketched in Fig. (6).
To further support the conclusion that the differential soil settlement is the main cause of the existent crack pattern, it is worth noting that front D and walls A, B and C far from the excavation zone do not present meaningful cracks.

ANALYTICAL INTERPRETATION OF THE CRACK PATTERN, WALL B
A simple analytical interpretation of the causes at the base of the formation of the crack pattern now visible on wall B, may be obtained evaluating, by means of simple equilibrium relations, the variation of the state of stress along  the wall consequent to a decrease of the pressure transferred by the soil to the front wall E.
With good approximation, such decrease near the excavation area results in a global redistribution of the pressures transferred by the soil, as schematically depicted in Fig. (6).In particular, in Fig. (6), pressure distributions assumed in the analytical model proposed in the following section before and after the excavation are compared.While such distributions are merely hypothetical, according to authors' experience in this field they seem quite reasonable.
The area where the pressure transferred is assumed null that has been geometrically determined knowing position, area and depth of the excavation.Assuming for the soil a friction angle Φ equal to 26°, all the soil contained between the excavated area and a plane inclined to the horizontal direction with angle of inclination equal to 2 / 45 Φ + ° is assumed to be unable to sustain structure gravity loads, signifying that a pressure equal to zero has to be assumed there.
From Fig. (6), where the zones unable to sustain masonry gravity loads are graphically determined, it is pretty clear how the soil under the principal front wall is totally unable to withstand vertical loads in the central region, however gradually increasing load carrying capacity near the edges.The pressure variation from zero (central point) to the undisturbed value (lateral edges) is hereafter assumed linear for the sake of simplicity.Such a redistribution of soil pressures requires that the main front wall is in equilibrium, thanks to the interlocking with perpendicular walls, but developing on the same perpendicular walls, through the schematization by means of a cantilever beam behavior, non-negligible bending.The existent crack pattern on such walls suggests that the visible pseudovertical cracks under the roof are consequent to a bending behavior of wall B, thus a-posteriori confirming the simplified hypotheses assumed here.Also the visible crack on the lateral part of the front wall A, while less wide, shows analogous origin.
On the contrary, no meaningful cracks are visible on edge wall C, probably for an eccentricity in the excavation zone present in reality.The different damage level of wall B with respect to the other transversal walls may be easily explained considering that central wall B is both subjected to higher dead loads and very close to the excavated area, especially at the intersection with the main front E.

SIMPLIFIED EVALUATION OF INTERNAL AC-TIONS ON MASONRY WALLS
As a consequence of the modified distribution of soil pressures, also internal stresses on masonry walls modify.A simplified evaluation of the internal stress distribution consequent to soil pressures re-distribution may be obtained assuming that the following hypotheses hold, see also Fig. (7):   • In the phase preceding the excavation, the equilibrium condition under dead loads (masonry self-weight and roof loads) is uniform, having assumed shear and bending actions acting in the plane of the walls equal to zero.Whilst this is a quite debatable hypothesis, it seems reasonable for the problem at hand, due to the regularity of the building.Elastic analyses conducted through a commercial code essentially confirm the effectiveness of such hypothesis.
• As a consequence of the excavation, the soil pressure near such region decreases.From a kinematic point of view, such a static assumption could be a-posteriori verified reducing the stiffness of the springs representing soil within a Winkler approximation.
• As a consequence of the different pressures transferred by the soil to the structure, it is further assumed that on vertical sections belonging to walls A, B and C shear and bending actions act.
• Masonry self-weigh in correspondence of the zone where soil pressure is kept equal to zero, is transferred to the contiguous walls by means of compressed struts belonging to wall E middle plane (see Fig. 6).
The vertical weight of the roof considered in the calculation is equal to 1.4 kN /m 2 , whereas for masonry a specific weight equal to 1.8 kN/m 3 is assumed.Both are typical values for existing masonry structures with light roofs.Within these hypotheses, a reduction of the soil pressure in correspondence to wall E and B intersection is assumed.To maintain walls in equilibrium, for any vertical section drawn on wall B, non-null bending moment and vertical shear are present.With simple equilibrium equations (vertical translation and rotation) it is possible to evaluate both bending moment M and shear V values at different positions of the vertical section with respect to wall E. Results of the analytical calculations, performed easily with an Excel sheet, are reported in Fig. (8).From an analysis of the figure, it is interesting to notice that M reaches a maximum at around 12 meters.Big bending moments are associated to big values of tensile stresses, meaning that M is a useful diagram to check the approximate position of the vertical cracks.Finally it is interesting to notice that V exhibits a maximum at 2.5 meters and is zero at 12 meters, in agreement with the bending moment diagram (M convex for x<2.5 and maximum value at 12 meters).
In Fig. (9), maximum tensile stress on vertical sections belonging to wall B at a distance x from the façade is depicted.For the sake of simplicity, the diagram is plotted assuming an elastic behavior for the sections, meaning that the maximum stress may be simply evaluated as M/J*H/2, where M is the bending moment, J the inertia and H the height of the section respectively.
The diagram of the maximum tensile stress in Fig. (9) shows visible peaks where actually the sub vertical crack activated.Local peaks correspond, obviously, to sections sssswwith reduced area, i.e.where the large stones used to The maximum traction value is largely higher than the masonry peak resistance along the horizontal direction, even considering the orthotropy of the masonry material [17], having for this kind of masonry typical values ranging between 0.4-0.9MPa, see following Section for homogenized numerical calculations.
After the activation of the crack at the extrados of the wall, the linear elastic hypothesis for the wall is not valid anymore and a more complex material model should be used.In any case, the simple equilibrium model proposed is able to predict with good approximation the position of the cracks.

Evaluation of the Soil Stiffness for Wall B
When the hypothesis to model the soil by means of equivalent elastic Winkler vertical springs is made, an analytical estimation of soil springs stiffness is possible by means of direct identification with experimental data on measured cracks wide.While this may be rather debatable, it is without any doubt the most straightforward and less expensive approach that may be proposed.
To estimate the soil stiffness for wall B to use within the non-linear numerical FE code, the following hypotheses have been assumed: • As a consequence of the excavation, the pressure change is assumed due to a decrease of the Winkler stiffness of the linear elastic springs representing the soil.The masonry building is considered much more stiff than the soil and, at a first glance, its movement on the foundation is characterized by only three degrees of freedom per wall, namely two in-plane translations and an inplane rotation.
• The stiffness of the soil springs under walls A and C is assumed always equal to the undisturbed state stiffness k 1 , with a progressive linear decrease up to zero under the façade subjected to overturning.
• The maximum measured opening of the sub-vertical crack on wall B (at the top under the roof) is assumed equal to w max =10 mm, in agreement with measured data.
Under such hypotheses, the rotation θ of the portion of wall B between the crack and wall E is equal to θ = w max L .For Hence, the result- ant bending moment with respect to axes origin, see Fig. (7) is constituted by two contributions, namely M 1 and M 2 , respectively the bending moment of the springs having stiffness equal to k 1 and M 2 , with linear stiffness variation from 0 to k 1 : Striking an equilibrium with dead and live loads M 1 +M 2 +M M equate gravity loads of the structure, being indicated by M M masonry resistant bending moment on the cracked section positioned at 11 m from the façade.From such equation, it has been calculated that the value of k 1 for the particular case under investigation is equal to 0.8 N/mm.

NUMERICAL HOMOGENIZED NON-LINEAR MODEL
The analytical solutions obtained with the equilibrated closed-form model previously described, especially when dealing with the values to assume for the stiffness of the equivalent Winkler springs have been utilized within existing FE non-commercial software, suitable for the analysis of complex masonry structures in the non-linear range.
The homogenization proposed, e.g.[16], pertains to running bond non-strengthened masonry, regarded as an assemblage of bricks interacting through interfaces (mortar joints).Bricks are assumed to be infinitely resistant, whereas for joints a Mohr Coulomb failure criterion with tension cut-off and compressive limited strength is adopted.In this way, a full description of the model can be given at a micro-scale considering a representative volume REV constituted by a generic brick interacting with its six neighbors, see Fig. (10).A sub-class of possible elementary deformation modes acting in the unit cell is a priori chosen in order to describe joints cracking under normal, tangential actions and bending.Then, a numerical procedure of identification between the 3D discrete system and a continuum 2D equivalent model is proposed, equating internal work expended by the two models.

Heterogeneous and Macroscopic Homogeneous Model
In the heterogeneous model, the whole REV is meshed through six-noded wedge elements interconnected by interfaces (internal brick-brick interfaces and mortar joints, see Fig. 10-a).The motion of a generic element E , see Fig.  ( )( ) When dealing with the continuous model, a standard Cauchy bi-dimensional continuum, Fig. x , 2 x and 3 x .
The displacement field of a point P (coordinates [ ] ) belonging to the equivalent continuum plate is given by fields ( ) T and 23 T ) contribute to the internal work.In particu- lar, the work dissipated by an equivalent plate model is simply: where E is the in-plane strain vector, χ the out-of-plane strain vector and γ the out-of-plane shear strain.

Simplified Homogenization
To substitute the heterogeneous material with the homogeneous equivalent 2D model, a simple, compatible identification model is proposed [21], where the work expended by the blocks model, is equated to the work (5) by the equivalent model.
Here, fields The application of equation ( 6) to the heterogeneous model permits to directly determine displacements to apply to the boundary surfaces of the REV.
Brickwork behavior in flexion is obtained by integration of in-plane actions at a structural level, therefore, at the micro-scale it is possible to limit the study to in-plane and outof-plane shear actions ( E , γ respectively).

Non-linear Behavior of the Interfaces at a Meso-Level
At the meso-scale level, joints behave basically as Lourenço & Rots [22] interfaces.Their behavior is essentially elasto-plastic and bi-dimensional, exhibiting softening.
The macro-scale behavior is obviously different because of the utilization of homogenization concepts to derive the macroscopic stress-strain relationship of the interfaces.
We denote with E b and E m the elastic moduli of bricks and mortar respectively and we consider a masonry pillar constituted by two half bricks (height: H/2) and a central mortar joint (thickness: e h ).It has been shown by Lourenço & Rots [22] that, by making the deformation of the actual pillar equal to a simplified one constituted by elastic bricks (height: H/2+e h /2) and a joint reduced to interface, the stiff- . The same consideration may be repeated for the shear stiffness t k , i.e.

( )
having defined b G and m G as brick and mortar shear moduli.
Brick-brick interfaces connect rigid elements having the same material properties.An elastic-perfectly plastic material is assumed for bricks (Mohr-Coulomb failure criterion).It is also assumed that the common edge 12 Γ is subdivided into rectangular small elements of area e A 12 Γ .Following Kawai [23], at each area For a mortar interface, the elastic domain is, in the most general case, bound by a composite yield surface that includes tension, shear and compression failure with softening (see Fig. (10-d).A multi-surface plasticity model is adopted, with softening in both tension and compression.
For multi-surface plasticity the form of the elastic domain is defined by each i-th yield function 0 . Let us consider a planar mortar interface and a point on the interface where a normal stress σ and two mutually perpendicu- lar shear stresses 1 τ and 2 τ act.
The i-th yield functions are of the form , where scalar i κ rules the amount of softening of the i-th yield surface and i Φ and i Ψ are generic functions representing respectively the initial i-th yield surface and the correction which accounts for the evolution of the strength during the inelastic deformation process.
The usual elasto-plastic equations for single surface plasticity hold; assuming the hypothesis of small deformations, the total strain rate ε  is decomposed into an elastic compo- For any corner of the proposed model two yield surfaces are active and the previous equations must be suitably stated for multi-surface plasticity.Details can be found in, e.g.[22].
In order to model the failure of the joint, a classical Mohr-Coulomb type strength criterion is used with a tension cut-off and a linear compression cap, Fig. (10-d).The parameters f t and f c are, respectively, the tensile and compressive Mode-I strength of the mortar or mortar-brick interfaces, c is the cohesion, Φ is the friction angle, and Ψ is the angle which defines the linear compression cap.
For the tension mode, exponential softening on the tensile strength is assumed according to the mode I experiments conducted by many authors.The yield function reads: where the yield value ( ) deteriorates in agreement with the following formula: ( ) where 0 t f is the initial joint tensile strength and I f G is the mode I fracture energy.An associated flow rule is assumed here.
When dealing with the shear mode, a Mohr-Coulomb yield function is adopted:  When dealing with the linearized compressive cap, a three function model as the one proposed in [22] is utilized.The typical hardening/softening behavior of the law adopted is shown in Fig. (10-d), where the subscripts e, m, p and r of the yield value c f denote respectively, the elastic limit, medium, peak and residual values.The peak value cp f equals the masonry compressive strength c f of the inter- face.

Numerical Simulations AT A Cell Level
This section provides an insight into the inelastic behavior of a masonry REV in running bond texture and is constituted by common solid clay Italian bricks (dimensions 250 × 125 × 55 mm).While the actual Siloteca texture is a triple head patter, here for the sake of simplicity we assume that the walls are constituted by three non-interconnected leafs.Obviously, this assumption is quite simplistic and a detailed survey of the actual interconnections along the thickness should be required.However, the assumption made, i.e. no interconnection between walls, allows to conclude that the simplified approach proposed provides, at least for out-ofplane actions, lower bound estimations of the actual strength of the structure.
Elastic and inelastic material properties are summarized in Table I.Two different values of fracture energy G I are assumed.The first corresponding realistically to existing masonry (Case A), the second assuming an almost perfect plastic behavior in tension (Case B).The behavior in uniaxial tension is depicted in Fig. (11-a) for horizontal and verti-cal tension.The anisotropy of the homogenized model is particularly evident and is mainly due to the contribution in horizontal tension of the bed joint, which fails in shear.In order to validate the results, the curves obtained using classic FE simulation [15] performed on a mesh with 384 elastic plane stress quadrilateral elements and mortar elasto-plastic interfaces are also represented, indicated as "FEM model".As it is possible to notice, the agreement is almost perfect, even in the softening range.This is not surprising because fracture lines concentrate on joints reduced to interfaces, as demonstrated by the REV deformed shape depicted in Fig. (11-b), where normal stress-shear masonry interfaces damage maps are also reported for the sake of completeness.A very similar behavior is experienced in horizontal bending, as can be noted by deformed shape and interfaces damage patch reported in Fig. (11-c).For compression loads, the anisotropy is less evident, due to the low shear strength of the joint when compared to the compressive strength.Hence, minor differences are expected when comparing the horizontal and vertical compression.Compressive behavior is not reported for the sake of conciseness and because it is much less important in foundation settlement problems.The reader is referred to [3] for a discussion of this issue.

Macro-Scale (Structural Level): A Simple Sequential Quadratic Programming -SQP-Approach
Under some general hypotheses holding for materials exhibiting an elasto-plastic behavior, as for instance that the plasticity condition is piecewise-linearized with r linearly elastic-plastic interacting planes in the space of superimposed stress and strain components, that the unloading of yielded stress-points does not occur and the continuum is discretized into finite elements, it has been shown in classic papers, e.g.[24], that the solution of an elasto-plastic problem can be achieved solving the following equivalent quadratic programming: where E D is the assembled elastic stiffness matrix, E ε ( E pl ε ) is the assembled elastic (plastic) part of the total strain vector E t ε , E N is the shape functions matrix of the used finite element, E λ is the plastic multiplier vector, E H is the harden- ing matrix and E σ the assembled stress vector.
The FE model utilized next to analyze in the non-linear range masonry structures relies on a discretization through six noded wedge elements, assumed rigid and infinitely resistant, and quadrilateral interfaces where all deformation occurs (linear and non-linear).
Each interface connects rigid elements representing the same homogenized materials.For each interface, three displacement and two rotational non-linear springs are utilized, as schematically shown in Fig. (12).The third rotational spring, acting along an axis normal to the surface, is assumed rigid and infinitely resistant.To properly take into account some distinctive aspect of masonry behavior in flexion (dependence of the flexural behavior by in-plane compression), but limiting to a great extent the number of optimization variables involved in the QP scheme, the procedure envisaged in Fig. ( 12) is adopted for each interface.
For the sake of brevity, the focus is laid exclusively on bending moment acting on an interface k (a similar procedure is adopted to handle torsion), at an iteration (i) of the loading process, bending rotation Φ n (i−1) and normal displacement of the interface centroid δ n (i−1) of the previous iteration (i-1).the displacement field δ n y t 2 ( ) is therefore immediately calculated along the interface thickness (abscissa y t 2 ).For each interface, depending on its orientation with respect to the blocks disposition, the homogenized stressstrain behavior is investigated from the meso-scale.At each assumed strain ε n , an interface displacement at the macro- scale is univocally associated simply by applying what is stated in [23].In particular, given ε n , the corresponding displacement in the discrete model on the interface k between elements M and N is δ n = 1 / 2 H M + H N ( )ε n , where all symbols have been already introduced.
For the interface k the homogenized stress-displacement relationship is therefore known for each point of the interface.By integration with a reasonable subdivision along the thickness into layers (authors experienced that the utilization of 10 layers represents a good compromise between numerical efficiency and accuracy) the compression load N (i−1) on the interface at the (i-1)-th iteration is calculated.At a fixed value of membrane normal force, the non-linear relationship moment-curvature is investigated again from the meso-scale, along with its linear stepwise constant approximation (necessary to use the sequential quadratic programming scheme discussed in the sequel).Again it has been noticed that the passage between curvatures and rotations, necessary when a discrete representation at a structural level is adopted, is trivial and is due to Kawai [23].
In this way, bending moment and torsion may be evaluated step by step during the deformation process simply by integration.
A database of moment-curvature diagrams at different levels of normal stresses is always at disposal from mesoscale computations before any structural non-linear simulation.When normal membrane force is within the range inspected but does not exactly match with the values investigated, a linear interpolation law for the diagrams is used.In order to utilize sequentially the QP approach [11] an approximation of the non-linear behavior through a linear piecewise constant function is used.
Following this procedure, the resultant mechanical model is thus composed of 5 elasto-plastic springs, Fig. (12).Within each iteration, an elastic-perfectly plastic approximation for each spring is utilized, signifying that 10 plastic multipliers (two for each spring, λ + and λ − , corresponding to posi- tive or negative kinematic variables) for each interface are needed.In this way, optimization variables entering into the QP problem are relatively small (10 plastic multipliers for each interface, 3 displacements and 3 rotations for each element).
Since each interface is modeled through elasto-plastic uni-dimensional springs independent from each other, it can be stated that H E in equation ( 11) turns out to be diagonal.
Within the FE model adopted, problem (11) may be rewritten as: Assuming that the structural model has n in interfaces and n el elements, symbols in equation ( 12) have the following meaning: 1. K el is a 6n el × 6n el assembled matrix, collecting elastic stiffness of each interface.Local elastic stiffness matrix of each interface is obviously diagonal, whereas the global stiffness matrix K el is generally not diagonal.It is worth remembering that elastic stiffness values are evaluated at the meso-scale, as discussed in the previous section.
3. K ep is a 10n in × 10n in assembled matrix built from diag- onal matrices of hardening moduli of the interfaces.
4. U el is a 6n el vector collecting the displacements and rotations of the elements.
5. F is a 6n el vector of external loads (forces and moments) applied on element centroids.
Typically, the independent variable vector is represented by element displacements U el and plastic multiplier vectors + λ and − λ .
As usually done in a non-linear structural analysis, QP problem (12) is solved in terms of displacement and plastic multipliers step increments.The initial trust independent variable vector is always represented by the solution at the previous step.
In the framework of the two-step approach proposed, the non-linear behavior of the springs is approximated using a linear-discontinuous piecewise constant function, as recommended in [3].The softening behavior is handled within a Sequential Quadratic Programming SQP scheme, already discussed in [2,3], where the reader is referred for further details.

NUMERICAL RESULTS AT A STRUCTURAL LEVEL
The numerical model previously described and already presented in [2,3] without the presence of translational elastic springs representing soil stiffness, has been here generalized in the presence of a Winkler model.
Before presenting results obtained with the FE approach proposed, here it is worth noting only that, in the model, the following hypotheses/simplifications have been made: • a homogenization approach is assumed for masonry, highlighting that a non-linear material exhibiting differ- (y ) known at iteration (i)

Bending and torsion
Normal axion and shear ent strength in tension and compression and along vertical and horizontal direction is adopted.Mechanical properties of such material are derived fitting as close as possible experimental data available for the masonry under consideration, through homogenization theory.The inelastic behavior of masonry, see Fig. (11), realistically reproduces the actual properties of the masonry material under consideration [1], also reflecting the expected orthotropy ratio for the texture considered and the actual disposition and geometry of the bricks (Italian bricks of dimension 250×120×55 mm).
• To realistically reproduce the experimental crack pattern, it is assumed that Winkler stiffness for the springs placed at the foundation level remains unchanged and that dead loads and self-weight are incremented step by step from a null initial value to their actual value at the last step.The non-linearity of the masonry material results into non-linearity of the global behavior of the structure, with consequent cracks formation on the mesh used.The resultant crack pattern found numerically should be compared with the existing degradation pattern, in order to establish if the approximate procedure proposed in the previous Section adequately estimates the overturning of the façade.
At this stage, a FE discretization with 5056 wedge elements is used, as schematically shown in Fig. (13) on the deformed shape of the structure obtained at the end of the non-linear static analysis.A detail of the crack pattern is sketched in Fig. (14).Finally in Fig. (15) point A crack horizontal opening vs vertical loads multiplier is represented.
Point A is positioned under the roof at a distance equal to around 10.50 m from the façade subjected to rocking.From numerical results obtained, the position of the vertical cracks in the central wall, the façade crack pattern and the opening of the crack under the roof seem in very good agreement with experimental evidences, meaning that the combined analytical numerical approach proposed may provide quite accurate predictions of the state of degradation of the building.In order to further assess numerical results obtained with the non-standard procedure proposed, in Fig. (15) an additional crack opening vs vertical load multiplier curve is ob-tained with the commercial code Strand 7, assuming for masonry an isotropic elasto-plastic material obeying a Drucker-Prager failure criterion (with tensile strength kept as the average of the vertical and horizontal tensile strength).Whilst this latter approach is only able to provide crude approximations of the actual masonry behavior in the inelastic range, it seems however the only alternative procedure available in common design to assess present numerical results.

CONCLUSIONS
A combined analytical numerical model regarding the evaluation of the causes at the base of the state of degradation of an ancient masonry built aggregate has been presented.The approach proposed is a two-step process.In the first phase, an equilibrated analytical model constituted by rigid elements subjected to overturning has been proposed.The model allows a reliable evaluation of the position of the cracks determining overturning of the façade and part of the perpendicular walls.By simple identification with data regarding the actual cracks opening, elastic stiffness of the soil is evaluated, to be used within a sophisticated FE approach accounting for damage in the masonry elements.The FE model so tuned provides quite accurate results in terms of both crack pattern and determination of crack opening.

Fig. ( 1 )
Fig. (1).Milan, part of the museum with the riding-stable "Cavallerizze" on the right and the "Siloteca stables" on the left.

Fig. ( 5 ).
Fig. (5).Detail of the sub-vertical crack at the top of lateral wall B, fronts B and B'.

Fig. ( 9 )
Fig. (9).Wall B, maximum tensile stress on vertical sections at a distance x from the perpendicular wall E.

( 10 )
, is described as a function of its centroid ( C E ) dis- placement vector u E (components u xx E , u yy E and u zz E ) and of its rotation vectorE Φ (components Φ xx E , Φ yy E and Φ zz E )around centroid.When two contiguous bricks M and N are considered, the displacement of a generic point P in a position 12 the common interface between the two elements) is:( ) Fig. (10).-a: FE discretization of the REV.-b: Rigid infinitely resistant six-noded wedge element used for the REV discretization.-c:Γ 12 interface between contiguous elements.-d: Modified Mohr-Coulomb criterion for the mortar joint reduced to interface (left) and hardening/softening law in compression (right) as a function of the inelastic parameter κ 3 .
as a combination of elementary deformations in the unit cell, corresponding to actual failure mechanisms occurring, according to experimental evidences, in the presence of running bond brickwork with weak joints reduced to interfaces.From a practical point of view, fields sub-class of regular motions are obtained assuming alternatively one component of vector E , χ or γ unitary and setting all the other components equal to zero, subsequently choosing the most simple polynomial expressions for procedure described, rotations and displacements of each element belonging to the REV in the heterogeneous model are determined solving a boundary value problem on the REV where displacements (or displacement increments) on the boundary are imposed.For instance, when only 11 χ ≠0 is applied on the REV, a choice for ( )

νL
Kawai [23], their stiffness may be evaluated as is the Poisson's ratio, M H and N H are heights of triangles[3] and all the other symbols have been already introduced.
nent el ε  and a plastic component pl ε  .The elastic strain rate is related to the stress rate by the elastic constitutive matrix D= plastic potential corresponding to the i-th yield surface (which rules the direction of pl ε  in the stress space) and non-associated plasticity i g may not coincide with i f . where

G
is the mode II fracture energy and r φ tan is the residual friction angle, hereafter kept always equal to 75% of the initial one.A non-associated flow rule is assumed here,

Fig. ( 11 ).
Fig. (11).Masonry under consideration.-a: Uniaxial response of the homogenization model along horizontal and vertical tension for two values of fracture energy.-b: REV deformed shape at collapse for horizontal tension (mesh used and magnified view) with indication of interface damage in horizontal tension (center) and vertical tension (right).-c: same as previous, but for horizontal bending.

Fig. ( 12 ).
Fig.(12).Evaluation of the non linear load-displacement (or moment-rotation) behavior of the interfaces at each load step.