Determination of the Numerical Parameters of a Continuous Damage Model for the Structural Analysis of Clay Brick Masonry

Models based on the continuous damage theory present good responses in representing the nonlinear behavior of reinforced concrete structures with loss of strength and stiffness of the material. However, damage theory is rarely employed in the analysis of masonry structures and numerical simulations are currently performed mostly by Finite Element Method formulations. A computational program was designed to determine the numerical parameters of a damage model of the physical properties of masonry components, solid clay brick and mortar. The model was formulated based on the composition of tensile and compressive surface strengths in the plane stress state. The numerical parameters, the corresponding curves of the activation surfaces and the evolution of the surfaces are presented. The results were fed into the computational program based on the Boundary Element Method (BEM) for the simulation of masonry walls, and two types of masonry were simulated. The results confirm the good performance of the model and the program based on the BEM.


Introduction
Numerical modeling of the behavior of structural masonry, allied to laboratory results, can be extremely useful for understanding and predicting the performance of masonry structures.Numerical modeling requires the characterization of the component materials of masonry and mathematical modeling of the behavior (linear or nonlinear analysis, aspects of homogenization and anisotropy).The structural behavior of masonry considering geometrical aspects (load versus bed joint directions) has been studied since the pioneering works by Page 1 and Dhanasekar et al. 2 .
Masonry walls can be built in several ways (Figure 1a-e), resulting in different dimensions and properties of masonry.
However, a standard geometrical arrangement (Figure 1) can be employed to estimate the elastic behavior of the brick-mortar unit on each axis, since this arrangement of units is repeated, simplifying the mathematical treatment for structural analysis.Thus, the homogenization technique was used here to translate the behavior of brick-mortar masonry, a heterogeneous solid, into a homogeneous equivalent Representative Volume Element (RVE).
The nonlinear behavior of the material, including aspects of strain, microcracks and rupture, increase the level of complexity in numerical analysis.Basic and important works about these aspects have been published by Page 8 , Lourenço 9 , Lee et al. 10 , and others.Numerical models based on continuous damage theory provide a good representation of the structural behavior and allow for the prediction of microcrack growth and gradual loss of stiffness of the material in response to increased mechanical loads.Other important works include those of Chow and Wang 11,12 , Maier et al. 13 , Papa 14 and Pegon and Anthoine 15 , all of which used the Finite Element Method (FEM).
However, there are a few studies involving the use of formulations of the Boundary Element Method (BEM) to simulate masonry structures by nonlinear analysis.The applications used by Alessandri and Brebbia 16 and Rashed et al. 17 are BEM formulations in the plane stress state, similar to the one used in this work.The former study simulates masonry walls subjected to horizontal loads, and considers homogeneous materials in elastic-plastic behavior with infinite compressive strength and zero tensile strength.The latter study evaluates the failure modes of a masonry wall under concentrated vertical loading.The same incremental-interactive procedure is employed in an inverse analysis to find unknown values of initial stresses when the failure modes and corresponding stress-strain correlations are known.However, these two studies have not been followed by sequels.It should be noted that the cited works only simulate walls according to the application of a formulation of the numerical method, without providing further information or a discussion about the numerical parameters and their correlations with the physico-elastic properties of the materials.
The current work, therefore, presents the application of a nonlinear boundary element formulation with a continuous damage model, as has been proposed in previous works 18,19 .

Material Properties
In this work, the parameters of a damage model are obtained directly from the material's properties.The ABNT (Brazilian Association of Technical Standards) NBR 7170 [20] standard classifies clay brick in three categories: A, B and C, while the NBR 13281 [21] standard classifies mortar into classes P1 to P6, according to its compressive strength.

Clay brick
The NBR 7170/1983 [20] standard defines solid brick as "brick that has all its sides filled with material, which may have a recessed part containing the manufacturer's logo on one of the sides with a larger area."It also specifies that it should be made of extruded or pressed clay, fired at a temperature that allows the end product to meet the conditions established in the standard.
Common clay bricks are classified into three categories: A, B and C, according to a minimum compressive strength of 1.5 MPa, 2.5 MPa and 4.0 MPa, respectively, as indicated in Table 1.
Solid clay bricks have nominal dimensions of 19 cm length × 9 cm width × 5.7 cm thickness, or 19 cm length × 9 cm width × 9 cm thickness, and the maximum manufacturing tolerance for common bricks in the three dimensions is 3 mm larger or smaller.

Mortar
The NBR 13281/2005 [21] standard defines mortar as: "a homogeneous mixture of fine aggregate, inorganic binder and water, with or without additives or additions, with adhesive and hardening properties, which can be prepared on-site or at an industrial facility (industrial mortar)."This standard classifies mortar according to its compressive strength (see Table 2).

Damage parameter
Considering a damaged solid with the S ef section as a portion of the total area (S ef ≤ S) which effectively resists load F, the difference defines the area of defects: By definition, the local measure of damage, D, can be represented as: By definition, the effective stress, σ ef is When the material remains intact (D=0), i.e., undamaged, then σ=σ ef , which indicates that the stress acting on the material remains unchanged.However, if D is close to 1, σ ef increases indefinitely, causing the material to collapse.
Strain, ε, can be written as: ( ) where, as the damage spreads, Young's modulus varies from E to E ef : ( ) Analogously, if D is close to 1, E ef approaches 0, i.e., the material shows loss of stiffness due to straining.The evolution of the damage directly affects the material's elastic response, which is indicated by its reduced stiffness and strength.
A continuous damage model for concrete formulated by Comi and Perego 22 was used in the present study.This model is represented by a failure criterion of the material in the plane stress state (Figure 2), which is defined by an ellipse (compression function, f c ) and a hyperbola (tension function, f t ) that satisfy Equations 6 and 7: The expressions that relate the parameters of the model to the properties of materials (tensile and compressive stresses) are as follows:

Numerical parameters of the proposed damage model
In this work, the authors created a numerical program that determines the parameters of the damage model for the material properties established by the aforementioned Brazilian standards.Based on Equations 6 and 7, the basic damage parameters for each class of strength are a t , b t , k t (hyperbola) and a c , b c , k c (ellipse).The parameters h t and h c are calculated as a function of the variables of damage D t and D c , respectively 22 .
The tensile strength of the materials is considered to be as follows: • According to Lourenço, Almeida and Barros 23 , a brick's tensile strength is equal to approximately 5% of its compressive strength; and • The tensile strength of mortar is assumed to be equal to 10% of its compressive strength.In the present research, the authors sought to determine the values of the damage parameters D 0c and D 0t that correspond to the material's maximum strength in the hardening process.After compiling Equations 6 and 7 in the computer program with several values for the two damage variables (D c and D t ), the authors determined the maximum expansions for both activation surfaces, finding values of D c = D 0c = 0.555 and D t = D 0t = 0.3, which were used to plot the maximum surfaces.
Figures 3a, b, respectively, illustrate the graphs of the maximum activation surfaces for each strength class of brick and mortar.In these graphs, an elliptical surface delimits the minimum compressive stresses (blue curves, third quadrant), while a hyperbolic surface delimits the maximum tensile stresses (red curves, 1 st , 2 nd and 4 th quadrants) in the composition of a failure criterion related to a plane stress state in principal directions.These graphs show the maximum compressive stresses on the principal stress axis for each class of strength (brick or mortar).
For example, to understand the damage model, Figure 4a depicts the maximum activation surface for class C clay brick, obtained for D c = D 0c = 0.555 and D t = D 0t = 0.3.Note that the maximum strength values are reached on the stress axis at a compressive stress of 4 MPa and at a tensile stress of 0.2 MPa (5% of the compressive strength).
For a better understanding of the behavior of the damage model, Figure 4b illustrates the evolution of the activation surfaces when the damage variable D changes from 0% to 99.9%.
In the case of compressive stress activation surfaces (elliptic surface, Figure 4b), it can be noted that, as expected, the expansion reaches its maximum evolution at a uniaxial stress of 4 MPa, which corresponds to a damage variable of 55.5% (equivalent to the dashed red line, 60%).

Identification MPa
Compressive strength at 28 days (MPa) ≤2.0 1.5 to 3.0 2.5 to 4.5 4.0 to 6.The maximum expansion of the tensile stress activation surfaces occurs at a uniaxial stress of 0.2 MPa (hyperbolic surface), which is equivalent to a damage variable of 20% (solid blue line).Thus, based on the equations in Figure 2 compiled by the numerical program and the assumed values of the variables (D c and D t ), the numerical parameters for each class of strength values were determined as indicated in Table 3.
The efficiency of the formulation of the numerical program in determining the parameters of the properties of masonry materials (clay brick and mortar) was confirmed by its good representation of the structural behavior of the material in the evolution of damage resulting from increasing active stresses.

Boundary element formulation for two-dimensional nonlinear analysis
After suitable adjustments in the computational code, the Boundary Element formulation presented by Botta et al. 24 was reoriented to analyze structural masonry in 2D stress analysis 18,19 , and the parameters of Comi's damage model were implemented in the code.
Following the BEM formulation, after discretization of the boundary and the domain, the algebraic representations of integral equations of displacements (8) and internal stresses (9) results in: and where {U} is the displacement vector and {P} is the vector of boundary surface forces of the problem.The vector {σ 0 } represents the initial stresses in the domain.Matrices or in which {X} and {Y} are the vectors with unknown and imposed boundary values, respectively.Applying the same procedure to the stress Equation 9 and isolating the unknown boundary variables in vector {X}, one has: or Considering the definition of the initial stress as being the difference between the elastic stress and the real stress, and adding { } 0 σ to the two terms, one has: The equilibrium equation can then be rewritten, in terms of strain, as follows: which is a general expression for any nonlinear constituent law.
The algebraic equation must be calculated by means of an incremental-iterative procedure for nonlinear problems in which, for a loading step with an increment of ∆ε n , one has:

Proposed Discretization Procedure
To optimize the numerical procedure, internal points corresponding to the vertices of cells at the boundary are defined inside the cell (Figure 5).Thus, initial stresses are calculated only at internal points, allowing the calculation to be performed numerically or analytically for both domain integrals from Equations 8 and 9 for a resulting non-singular kernel, which is used to calculate the terms of matrices [Q] and [Q'].
The domain integrals are calculated numerically based on the Cauchy principal value sense 24 .To obtain the algebraic representation, linear approximations were used to transform the boundary variables and stresses over the triangular cells.With regard to the homogeneous medium in which the isotropic elastic tensor is defined, different damage parameters were adopted for both component materials, i.e., brick and mortar.It is also necessary to define   The discretization mesh adopted here is more refined in the region where the highest damage is expected.The mesh contains 52 linear boundary elements, 692 triangular cells and 646 internal initial stress points.Figure 9 shows the mesh, in which the cell vertices corresponding to internal points are displaced from their original location in order to highlight them (enclosed region).points inside the cells (discontinuous cells) whose sides coincide with the interfaces of different materials (Figure 5).
This proposal implies a high computational cost due to the large number of internal points created, which considerably increase the size of the [Q] and [Q'] matrices and the calculations of the load functions.Therefore, it should be used to simulate small-scale structures such as prisms or small masonry panels such as those used for shear testing.

Examples
This example is similar to a numeric example presented on Botta et al. 24 in which a rectangular plate is subjected to a uniform longitudinal load distributed over two opposing edges.Herein the plate is considered inhomogeneous, shown in Figure 6, and composed of a brick unit and a half thickness mortar joint.
The brick unit has dimensions of 220.0 × 52.0 mm 2 and the mortar is 5.0 mm thick.The elastic properties of both materials are E = 36000.0MPa and ν = 0.15 (homogeneous material).The discretization mesh (Figure 6) has 52 linear boundary elements, 480 cells and 215 internal points of initial stress.
The numerical data are: Four sets of brick-mortar specimens were considered, two with homogeneous characteristics and two with non-homogeneous characteristics, as described in Table 4.
Figure 7 illustrates the two simulations of the longitudinal stresses along the horizontal axis of symmetry plotted against longitudinal strains in each case.
Figure 7a indicates that the peak stress and the softening phase were determined by the strength of the mortar, i.e., the specimen sets with weak homogeneous material and with weak mortar showed the same stress/strain behavior.Figure 7b shows the same conclusion for the shrinkage process, but the peak stress of mortar and the softening phase could not be reached.Note that the elastic-linear phase occurs with Young´s modulus equal 36,000 N/mm 2 and maximum stress at 1,80 N/mm 2 .

Clay brick wall subjected to a uniformly distributed load at the top and controlled vertical displacement
The example considers a solid clay brick wall (Figure 8a) with dimensions of 757 × 457 mm 2 , supported on each side upon 188 mm of length, built with half-scale bricks of 122 × 37 × 54 mm 3 and 5-mm-thick mortar joints.The load at the top, P, is applied by means of a stiff 381-mm-long steel beam (Lourenço 9 ).

Conclusions
This paper proposed a computational program to determine the numerical parameters of a continuous damage model of the strength properties of masonry materials, i.e., solid clay bricks and mortar.A computer code was created to determine these parameters, and a failure criterion was developed for tensile and compressive surfaces in a plane stress state in principal stress directions.The proposed numerical program also enables visualization of these surfaces and their evolution when the damage variable D changes.The evolution of the activation surfaces offers a very good representation of the material's behavior in terms of loss of strength and stiffness in response to the application of additional stresses.
The results were fed into the numerical program based on the Boundary Element Method to simulate brick walls.The novelty in this computer program is the definition of two different regions, each of them with different damage parameters corresponding to different materials, and the approach that determines discontinuos cells adjacents to boundary or interfaces.
Two examples were presented.Simulations of a brick-mortar specimen and a masonry wall resulted in the expected stress evolution, demonstrating that both computer codes performed well and have a promising potential for simulations involving physically nonlinear analysis based on a damage model.The maximum tensile stresses (σ 1 = 0.56 N/mm 2 , Figure 8c) are located mid-point at the bottom of the wall span (Figure 8b).The present work determined a value of σ x = 0.292 N/mm 2 based on Figure 10a (yellow area).
In the most damaged region, the tensile damage parameter D t reached a value of about 0.586, as indicated in Figure 10c.This minor evolution in tensile damage, without evidence of compressive damage (D c remains null due to the low compressive stresses developed along the entire wall), explains the quasi-linear shape of the load vs. vertical controlled displacement curve.Figure 11b shows that the maximum load reached 85kN, as a result of the higher tensile damage parameter D t and the vertical displacement of 0.50 mm.This load level, which is still close to the experimental collapse load (109 kN) but well below the peak stress simulated numerically, is represented by a critical stage after which uncontrolled compressive damage increases on the inner sides of the supports (Figure 8c).The post-peak can only be simulated by arc-length constraint loading, which has not yet been done by the authors.

Figure 2 .
Figure 2. Constitutive equations of the damage model and linear elastic domain in the principal stress plane 22 .

Figure 3 .
Figure 3. a) Activation surfaces of brick strength classes (stresses in MPa); b) Activation surfaces of mortar strength classes (stresses in MPa).

Figure 4 .
Figure 4. Class C brick: a) Linear elastic envelope; b) Evolution of activation envelopes.
[H], [H´], [G] and [G´] are obtained from integrations over boundary elements and matrices [Q] and [Q´] are obtained from integrations over cells.After imposing boundary conditions, the displacement Equation 8 results in:

Figure 5 .
Figure 5. Internal points for initial stress calculation in cell vertices at boundaries or different material interfaces.

Figure 9 .
Figure 9. Discretization adopted in this work.Discontinuous cells adjacent to boundaries and interfaces; refined mesh in most damaged area (detail); (Paraview 3.2.1 Visualizer).
The numerical results presented by the referenced work are shown in Figure8c.The principal stresses have minimum values concentrated at the ends of the steel beam and the inner side of the supports (between σ 3 = -14.00N/mm 2 and σ 3 = -37.92N/mm 2 ).The results obtained in the present work are depicted in Figure10.The stresses are σ x = -0.430N/mm 2 and σ y = -0.878N/mm 2 at the inner sides of the supports and σ y = -1.29 N/mm 2 at the ends of the loaded region.The similarities of these results confirm the good representation presented by the BEM formulation in this simulation.

Table 2 .
Classification of mortar according to the ABNT NBR 13281/2005 standard.

Table 3 .
Constitutive damage parameters of Comi's model.

Table 4 .
Physical characteristics of the specimens.