Non-linear analysis of reinforced concrete plates by the boundary element method and continuum damage mechanics

In this work, non-linear analyses of reinforced concrete plates are performed by using a BEM formulation, based on Kirchhoff’s theory, which has already proved to be a robust technique to deal with plate problems. The non-linear behavior of concrete is modeled by the Mazars model, which is based on continuum damage mechanics (CDM) and has an easy parametric identification, while the reinforcement is governed by the uniaxial elasto-plastic model with constant hardening. Initially, the different types of damage models and their parametric identification are discussed and the Mazars model is presented. Then the BEM formulation is discussed, which is based on the initial moment technique, where the remaining domain integrals are evaluated by approaching the initial moment field over internal cells. The stress distribution over the plate thickness is obtained by using a Gauss scheme in which the adopted criterion is verified in each Gauss point defined along the plate. Then, the internal values of moments are approached by numerical integrals along the plate thickness. Finally, some numerical examples of reinforced concrete plates are analyzed, where the potentiality of the Mazars model is verified despite the difficulties of the parametric identification.


Introduction Introduction Introduction Introduction
The boundary element method (BEM) is already a well-established numerical technique to deal with an enormous number of complex engineering problems.BEM is particularly suitable to evaluate effort concentrations due to point loads.Moreover, BEM is a far better numerical technique to evaluate shear efforts when compared with other methods.BEM has the advantage of dealing with deflections, slopes, moments and shear forces, approaching them using the same order polynomials.Plate bending BEM analysis has also been successfully extended to non-linear problems, for which explicit techniques usually based on initial moments or curvatures have been adopted.Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 The first works discussing the use of boundary element direct formulation, in conjunction with Kirchhoff's theory, are due to Bezine in 1978and Tottenhan in 1979. Venturini and Paiva (1993) have shown several alternatives to define the algebraic system of equations by placing the collocations either along the boundary or at particular outside points.
All the works reported above deal with the linear plate bending problem.For elasto-plastic plate bending analysis, assuming Kirchhoff's theory and using BEM, we have to mention Moshaiov and Vorus (1986), in which a particular incremental scheme has been proposed.Fernandes and Venturini (2002;2007) have studied this problem as well, introducing approximated models to deal with concrete slabs.Non-linear boundary element models have been also developed in the context of moderately thick plates by taking Reissner's theory into consideration, e.g.Karan and Telles (1992).
In the present paper, new numerical examples using the non-linear BEM formulation to analyze reinforced concrete slabs proposed in Fernandes and Venturini (2002) are discussed.This work intends to show the potentialities of that numerical model (BEM and CDM) applied on new situations where the parametric identification of the damage model is discussed.Although the formulation can admit both the elasto-plastic and damage models for the concrete, in this work only the Mazars (PITUBA; PROENÇA, 2005) continuum damage model will be considered.For the steel bars used as the reinforcement, an elasto-plastic criterion with hardening has been assumed.Integral representations of displacements will be derived from Betti's reciprocal work theorem, taking into account an initial moment field which will be assumed as the corrector vector of the non-linear incremental-iterative procedure.The actual internal forces are computed numerically by performing properly the integral over the plate thickness, using a convenient number of Gauss' stations.The chosen criterion is verified at each Gauss' point leading to the actual and initial distribution of stresses along the plate cross section.Finally, two numerical examples are taken to illustrate the accuracy and stability of the presented formulation.

Basic equations Basic equations Basic equations Basic equations
The equilibrium of an infinitesimal thin plate element, in absence of distributed external moments, gives the following differential equation written in terms of internal forces, where: m ij are bending and twisting moments; q i represents shear forces; g is distributed load acting in the direction perpendicular to plate middle surface.
The plate domain, defined over the middle surface, is denoted by Ω , while its boundary is represented by Γ , along which the following generalized boundary conditions may be assumed: . The generalized displacements u i are deflections, w, and rotations, ∂w ∂n -1 , while generalized tractions are given by normal bending moments, M n , and effective boundary shear forces, V n ( ); (n, s) is the local co-ordinate system, being n and s then normal and tangential directions, respectively).
For non-linear problems, the total strain components are divided into two parts. where: in which ij δ is the Kronecker delta, υ is Poisson's ratio and the classical plate stiffness D F is given by: [ ]

Et
, with E being Young's modulus.It is important to note that in the present work constitutive models for the concrete will not be considered elasto-plastic, but only a damage model.
Although, in any case the plastic moments ( represent, in fact, the corrector moments that must be considered in the incremental-iterative procedure in order to achieve plate equilibrium. Acta Scientiarum. Technology Maringá, v. 31, n. 1, p. 25-33, 2009 By differentiating equation ( 3) and then replacing it into equation (1a), one can write shear forces in terms of deflection derivatives, as follows: Replacing equation ( 1) into (1) and then using equation (4), one can achieve the equilibrium final expression, for plates with constant stiffness D F , written in terms of displacements: iijj w,

(
) where: w w iijj 2 , ∇ = , stands for the Laplacian operator here applied on deflections w.

Integral representations Integral representations Integral representations Integral representations
The integral representations will be derived from the reciprocity relation, which can be written in terms of moments and curvatures, as follow: In equation ( 7) ) ( , q w ij and ) (q m e ij represent, respectively, curvatures and elastic moments related to the actual plate bending problem valid over the domain Ω, while w * (s,q) is the well known fundamental solution, i.e., the deflection at the field point q due to a unit load applied at a source point s.Differentiating (7) twice by parts and considering (3) one can derive the following deflection integral representation: where: S is the source point that can be placed anywhere, Ω g is the sub-region where the load g is applied; the free term C(S) is given by: C(S)=1 for an internal point, C(S)=0 for outside points and C(S)=β C /2π for boundary points, being β C the internal angle at S.
In equation ( 8), w, ∂w/∂n, M nn , V n and R C are actual plate boundary values, deflections, rotations, normal bending moment, effective shear forces and corner reactions, respectively.The values indicated by * are the corresponding ones derived from the fundamental solution.To obtain the curvature integral representations, one has to differentiate twice equation (8).Then, bending and twisting moment integral representations are obtained by simply applying the definition given in equation ( 4), resulting into: where: g ijkl is a free term resulting of differentiating a singular integral term, being given by [ ] The shear force integral representation can be obtained by differentiating once the curvature equation to apply the definition given in (5).

Algebraic equations Algebraic equations Algebraic equations Algebraic equations
To transform the integral representations (8), and (9) into algebraic expressions the plate boundary has been discretized into geometrically linear elements over which the boundary values have been approximated by quadratic shape functions.Besides, the plastic moments, ) q ( m p ij , are approximated by adopting linear shape functions over domain triangular cells. After selecting an appropriate number of algebraic equations, one can assemble a convenient set of equations to solve the problem in terms of boundary values.As it has already been observed in previous study (FERNANDES;VENTURINI, 2002), it is better to use only deflection representations to define the final set of equations, taking the collocation points either along the boundary or outside the body.We tested algebraic equations achieved by either selecting only outside points or outside and boundary points.Both schemes gave accurate results.Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 Writing equation ( 8) into its discretized form for a convenient number of collocation points, one finds the BEM classical matrix equation, Vectors U and P, in equation ( 10), contain the generalized displacement and traction nodal values; vector T represents the loading, and M p contains plastic (or corrector) moments at internal and boundary nodes.H and G are the standard square matrices achieved by integrating all boundary elements, while E is a matrix obtained by performing the integrals over all cells.
The moment integral representation (9) can also be written into its algebraic form: with: I being the identity matrix.All matrices appearing in equation ( 11) are similar to the ones in equation ( 10), obtained by using the corresponding kernel exhibited in equation ( 9).Equations ( 10) and ( 11) can be conveniently arranged to express the solution in terms of boundary values, moments and shear forces, as follows: In equations ( 12) L, N and N' give the elastic solution due to the prescribed loads acting along the boundary or over the domain, while M p effects are represented by matrices R, S and S'.

Stratified model Stratified model Stratified model Stratified model
To capture different states of stresses that occur along the thickness, due to the different behavior of the concrete in tension and compression, the criterion is verified at particular points defined along the thickness.Then, at any cross section the actual moment increments will be computed by integrating the stress increment distribution along the thickness using a Gauss scheme as follow: σ is the concrete stress distribution along the thickness, sk ij σ represents the steel stress of a layer placed at x 3sk , while A sk represents the steel bar cross section, being N s the reinforced layer number and δ ij the Kronecker delta, IG ξ is the Gauss point homogenous co-ordinate.
Equation ( 13) is used to compute the internal moment tensor, for a specific plate cross section.This value is computed assuming the plate's middle surface to define the neutral axis, therefore to deal with a symmetric stress distribution diagram.For concrete slabs, one must take into account that the material behaves differently when in tension or in compression.In addition, reinforcement is usually non-symmetrically distributed.Thus, the neutral axis is no longer defined by the middle surface.To continue the analysis in the context of simple bending, one needs to define the new position of the neutral axis enforcing the stress resultant to be zero (FERNANDES;VENTURINI, 2002).After that, by conveniently using equation ( 13) one can compute the actual and plastic moments or their increments, M ∆ and p M ∆ , to be adopted in the non-linear algorithm described in section 7.For damage analysis, the same procedure is followed, replacing plastic values by damaged values.

Remarks about damage constitutive models Remarks about damage constitutive models Remarks about damage constitutive models Remarks about damage constitutive models
The Continuum Damage Mechanics (CDM) is a tool for the simulation of the material deterioration in equivalent continuous media due exclusively to microcracking process.A material can be simulated as a continuous medium and the influence of the internal changes caused by the microcracks are considered through scalar or tensor damage variables, which implies in a reduction of several rigidity components, where the damaged material can keep its isotropic properties or to become anisotropic.In the isotropic models, the damage does not affect the number of either the symmetry directions and the initial symmetry planes of the material, i.e., if the medium is initially isotropic or anisotropic with some degrees, those characteristics are preserved during the damage process.
On the other hand, the anisotropic models have the ability to change the number of both the symmetry directions and the initial symmetry planes of the material.In the last years, many constitutive models have been proposed in order to take into account anisotropic characteristic of the medium (BRÜNIG, 2004;PITUBA, 2007).
In order to represent the mathematical formulation of media with anisotropy induced by damage, it is considered that the damaged rigidity Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 constitutive tensor can be composed by two additive parts: the first one represents the undamaged material and second one is dependent of the actual damage state.When that last part is written like a summation of the scalar valued functions of the damage scalar variables that multiplies the fourthorder constitutive tensors, the resulting model is socalled scalar damage model.For instance, simple additive forms can be expressed as follows: where: M and M are fourth-order tensors and λ[D] is a scalar valued function of the damage scalar variables.
For initially isotropic materials, it can preserve the isotopic characteristic or create an anisotropic effect on the constitutive tensors depending how the tensors M and M were defined.
In the other hand, the general forms to the fourth-order damage tensor D can be proposed in order to take into account the anisotropy induced by damage.In Pituba (2007) the definition of the damage tensor follows a so-called scalar form expressed as: D = f j (D i ) M j , where f j (D i ) are scalar valued functions of the damage scalar variables D i and M j are anisotropic tensors.
Finally, constitutive models with orthotropy and bimodularity induced by damage can be proposed in order to evaluate the potentialities in the analysis of reinforced concrete structures.However, it must be observed that the parametric identification could become unfeasible at the practical point of view.The parametric identification is essential for using a constitutive model.The values to be determined can present reasonable sensibility mainly in the cases where the parameters number is reduced, which is the case of the proposed model.The adoption of a small number of parameters leads to an easier identification.However the capture of the characteristics of the stress-strain experimental curve of the material becomes more complicated.On the other hand, to get a correct identification and a good adjustment of the model is indispensable reliable experimental results.
Therefore, in this work, one has chosen an isotropic damage model due to the formulation simplicity and its easy parametric identification.In this context, the damage constitutive model proposed by Mazars is a well-known model whose efficiency has been already tested in various situations.

Constitutive models Constitutive models Constitutive models Constitutive models
As already mentioned, for the reinforcement the steel behavior is governed by a uniaxial elasto-plastic curve exhibiting hardening effects.For the concrete we have considered the well-known damage model proposed by Mazars (PITUBA;PROENÇA, 2005).In this model, the damage is represented by the scalar variable D (with 0 ≤ D ≤ 1), whose evolution occurs when the equivalent extension deformation ε ~ is bigger than a reference value.The plastic deformations evidenced experimentally are not considered.
The equivalent extension deformation is given by: where: i ε is a principal deformation component, being The damage activation occurs when = ε ~εd0 , being ε d0 the deformation referred to the maximum stress of a uniaxial tension test.Thus the criterion is given by: Considering the principles of thermodynamics, the damage evolution can be expressed by: where: is written in terms of ε ~ and defined continuous and positive.
As the concrete behaves differently in tension and compression, the damage variable D is obtained by combining properly the variables D T and D C , related to tension and compression, respectively, as follows: Now, it is necessary to identify the evolution laws for the scalar damage variables.In a general way, the evolution laws can be obtained by: Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 -direct identification from experimental tests obtained at laboratories; -definition of a dissipation potential written in terms of the associated variables, that includes states that can be achieved without any additional energy dissipation.In this case, the evolution laws are obtained from the variation of the dissipation potential.
An inconvenience of the first strategy is its limited use.One can lead to evolution laws that have their application restricted to loading paths combinations similar to those of the experimental tests.The second strategy has a more theoretical stamp, though with some advantages such as: application to general situations of loading paths, constitutive models with formal structure similar to models derived from other theories, and the possibility of maintenance of the constitutive tensor symmetry.
In this work, in order to consider a more practical employment of the numerical model, one has followed the first strategy.Therefore, the evolution laws proposed for the scalar damage variables are resulting from fittings on experimental results (PITUBA; PROENÇA, 2005): In equations ( 19) A T and B T are parameters related to uniaxial tension tests while A C and B C .are obtained from uniaxial compression tests.The parametric identification results of constitutive models obtained from experimental tests in the concrete specimens are presented in Pituba and Proença (2005).
To compute the α T and α C values defined in (18), we have to obtain, initially, the deformations ε T and ε C associated, respectively, to tension and compression states as follows: where: I is the identity tensor, E the elastic modulus of a non-damaged material, where: ε is given by: Finally, the constitutive relation can be expressed in terms of the actual deformation tensor as follows:

Solution technique Solution technique Solution technique Solution technique
As the plate behavior was assumed to be nonlinear, an incremental-iterative procedure is required.Thus, for the iteration j of the load increment i, the elastic moment increments e M ∆ , for all boundary and domain nodes, are given by: In equation ( 23), β i is the load factor, N and S are defined in (12) and ∆M P is the plastic moment increment computed during the previous iteration.
To compute the actual moment field taking into account particular criteria for each concrete layer and for the reinforcement, we follow the scheme defined in section 5. Initially, using Hooke's law, the corresponding curvature increment is evaluated from the elastic moment increment.Then, at each Gauss point and at the reinforcement positions, strain and stress elastic increments are computed.As superposition is admitted, those values can be accumulated at proper vectors.After verifying the adopted constitutive model for all Gauss points and reinforcements, we are able to compute the three inplane components N x , N y and N xy .The actual strain and stress distribution along the plate thickness is Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 obtained by enforcing these components to be zero.Finally, the actual moment increment (∆M) vector is evaluated by using equation ( 13).The plastic (or corrector) moment increments are then computed at all boundary and internal nodes: The corrector value { } 1 is then applied to the next iteration as elastic moment increments (equation 23b).The iterative procedure ends when the convergence criterion is satisfied according to a tolerance previously defined.At the increment end, boundary values and shear forces can be evaluated using expressions (12), for which M P represents the accumulated plastic moment vector.

Results and discussion
In this section, two numerical examples of concrete reinforced slabs were chosen to verify the accuracy of the BEM formulation prescribed previously.This model was implemented using the FORTRAN code as well as the direct boundary element method based on collocation points.For both examples a conventional ultimate state, defined by ultimate strain values cu ε = -0.35%and su ε = 1.0%, respectively for concrete and steel materials, was assumed to define the ultimate load.
In the first example, the plate geometry and boundary conditions are particularly defined in order to simulate the behaviour of a 30 cm thick beam (t = 30 cm), as shown in Figure (1a).The single applied load taken to analyze this beam is given by the moment M = 500 kN cm, distributed along the simply supported sides of length 2b = 20 cm.The other two sides of length equal to 2a = 150 cm are free, with no prescribed displacement or rotation.The boundary was discretized into 8 quadratic elements, while 24 cells with the definition of 5 internal points were adopted to approximate initial moments over the plate domain, as shown in Figure (1b).The matrix H and G have been computed by defining boundary and outside collocation points, being the distance of the outside point to the boundary given by: d 1 = 0.1l, where l is the element length.The elastic modulus E c , for the concrete material, was assumed equal to 3,000 kN cm -2 , being its Poisson's ration ν = 0 and the concrete ultimate stress, f c = 3 kN cm -2 .The reinforcement is constant and displayed only in the x 1 direction The steel elastic modulus E s was assumed equal to 27,000 kN cm -2 with the yielding stress f y = 24 kN cm -2 and hardening modulus K s = 14,000 kN cm -2 .The total steel cross section A s =0.5 cm 2 cm -1 was considered for all bars that have been placed at x 3 = 12 cm.The convergence was controlled by a tolerance of 0.0001%.To compute the stress distribution along the plate thickness 12 Gauss points have been considered.
The load was applied in 36 increments.Elastic responses were observed until the load factor β reached 0.1, when the concrete was damaged.Plastic strains appeared in the steel bars as well, but only for β equal to 0.60.The analysis stopped at β = 2.7 (M = 1350 kN cm), when the assumed steel maximum strain of -0.01 was reached.Figure (2) shows the deflection at central point 23 through the incremental procedure.Note that for all increments the computed values for moments were the following ones: M 11 = 500 βkN cm, M 12 = M 22 = 0.0, as expected.It is interesting to note that this same example has been considered in Fernandes and Venturini (2002), adopting an elastoplastic model for the concrete.Although in Fernandes and Venturini (2002) we have also computed internal moments equal to M 11 = 500 βkN cm, M 12 = M 22 = 0.0, the ultimate load was obtained for the applied moment M = 450 kN cm, when the steel bars has reached the maximum strain.One can observe that this ultimate moment is much lower than the one obtained in the present work.
Load Factor -β w (cm) Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 The example was also analyzed adopting a finer mesh with 16 quadratic elements, giving practically the same results obtained with the first discretization.
The second example consists of analyzing a square plate of side length a = 18in and thickness t = 1.75in supported only at the corners.A distributed load g = 200 psi has been applied over a small sub-region at the center of the plate (see Figure 3a), resulting into a total load P = 3.2 kpis.The matrix H and G have been computed by defining only outside collocation points, whose distance to the boundary are given by: d 1 = 0.1l and d 2 = 0.25l, respectively, being l the element length.The plate was uniformly discretized by adopting 32 quadratic boundary elements and 512 triangular internal cells, giving 68 boundary nodes and 225 internal points (Figure 3b).The reinforcement is characterized by the steel elastic modulus E s = 29000000 psi and yielding strength f y = 50000 psi with no hardening.The bars cross sections in both directions are A s = 0,011135"²/" which are placed at x 3 = 0,435".The adopted concrete elastic modulus is E c = 4150000 psi , while its Poisson's ratio and ultimate strength are ν = 0.15 and f c = 5500 psi, respectively.The tolerance of 0.2% was assumed to govern the convergence criterion, while 12 Gauss' points were used to perform the integrals across the slab thickness.
The load has been applied in 26 increments, starting by β = 0.1.The analysis has stopped for β = 0.79 (P = 2.53 Kpis), when the plate equilibrium has not be achieved.The concrete started to damage for P = 1.12 Kpis, while the steel bars had plastic deformations for P = 2.24 Kpis.In Figure 4 is displayed the load-deflection curve for the point 1 (Figure 3a), adopting different limits for the damage variable D. One can observe that the best results have been obtained for D = 0.7.As there are not stress x deformation curves in tension and compression for the concrete, in order to identify the model parameters, one has adopted average values for these parameters and limited the damage variable value.Thus, this example intends to show the capacity of the numerical model (BEM and CDM) to simulate the behavior of reinforced concrete plates.

Conclusion Conclusion Conclusion Conclusion
New examples using the BEM formulation proposed in Fernandes and Venturini (2002) to perform non-linear analysis of reinforced plates have been discussed.Once again the technique has proved to be efficient to model the concrete element stiffness deterioration using damage models, for which stability was always observed within the range of practical applications; for these situations, ultimate loads have been easily obtained without showing any kind of instabilities in the computations.These numerical instabilities can be observed in many situations even for intermediate levels of the damage processes leading to a deformation localization phenomenon.In these cases, it is necessary to use non-local versions of damage models.BEM has already proved to be a Acta Scientiarum.Technology Maringá, v. 31, n. 1, p. 25-33, 2009 suitable numerical tool to deal with plate bending problems.The method is particularly recommended to evaluate concentrations of internal forces or deformations that very often appear in practical problems.Moreover, the same order of errors is expected when computing deflections, slopes, moments and shear forces.Shear forces, for instance, are much better evaluated when compared with other numerical methods.They are not obtained by differentiating approximation function as for other numerical techniques.
The inclusion of a non-local version of a damage model in the numerical model presented here can be used in future works for the development of other numerical model depending of the involved phenomena.
Finally, note that despite of having a limited parametric identification, because there was enough data about the concrete used in the structure, one can observe the potentialities of the damage constitutive model to simulate the concrete behavior.
, positive and negative parts of the deformations ε T and ε C defined in (20);

Figure 1
Figure 1.a) Plate geometry and loading; b) Plate discretization.

Figure 3
Figure 3. a) Plate supported at corners with a central load; b) Plate discretization.