**TECHNICAL PAPER**

**Numerical modelling for the evaluation of progressive damage to plain concrete structures**

**S A Sadrnejad**

**ABSTRACT**

**Keywords:** crack type, micro-planes, progressive failure, plain concrete structures

**INTRODUCTION**

Surveys through the continuum macroscopic constitutive models of concrete show that micro-planes-based models perform exceptionally well. In these models, instead of presenting constitutive relations in the shape of stress and strain tensors, the stress and strain vectors, which are in turn the projections of stress or strain tensors, are used. This method not only provides a more physical conceptual base, but also leads to more simple mathematical formulation. On the other hand, invariant-based continuum macroscopic models lose some of the important features of material behaviour because they are basically not able to capture and store the data properties in the different directions around a material point, whereas the micro-planes models inherently include the directional characteristics of a material point.

Many of the available 3D commercial analysing codes are not equipped with the proper damage concrete material model, so are not able to perform precisely crack/damage analysis, analysis of crack orientations and analysis of progressive failure of concrete structures.

The basic idea, namely that the constitutive material behaviour as a relation between strain and stress tensors can be "assembled" from the behaviour of material in the planes with different orientations within the material (such as slip planes, micro-planes, particle contacts, etc), can be traced back to the failure envelopes of Mohr (1900) and the "slip theory of plasticity" of Taylor (1938) who was the first to implement this theory for modelling the behaviour of polycrystalline metals. Taylor's idea was formulated in detail by Batdorf and Budiansky (1949). This theory was soon recognised as a realistic constitutive model for plastic-hardening metals. This idea was used in arguments about the physical origin of strain hardening, and was shown to allow easy modelling of anisotropy, as well as the vertex effects for loading increments to the side of a radial path in stress space. All the formulations considered that only inelastic shear strains ("slips"), with no inelastic normal strains, took place on what is now called the "micro-planes". The theory was also adapted to anisotropic rocks and soils under the name "multi-laminate model" (Zienkiewicz & Pande 1977; Pande & Sharma 1981; Sadrnejad & Pande 1989; Sadrnejad 1992).

The static constraint formulation was used extensively and referred to as "slip theory" for metals or "multi-laminate theory" for anisotropic rocks until the first application of this theory by Bažant and Gambarova in 1984 to continuum damage mechanics and cohesive-frictional materials; they referred to this theory as "micro-plane" theory.

The object of this study is the application of a proposed micro-plane, damage-based model of plane concrete (Labibzadeh & Sadrnejad 2008) in a 3D finite element code to show its abilities in crack/damage analysis and analysis of progressive failure of concrete structures such as concrete double curvature arch dams. The proposed code is not only able to predict the crack line, but also allows the determination of which combination of loading conditions occurs on damaged micro-planes. The validity of the proposed code is investigated by means of a few tests and examples of concrete structures.

]]>**MICRO-PLANE FORMULATION WITH KINEMATIC CONSTRAINT**

The orientation of a micro-plane is characterised by the unit normal **n** of components *n _{i}* (indices

*i*and

*j*refer to the components in Cartesian co-ordinates

*x*). In the formulation with a kinematic constraint, which makes it possible to describe the softening behaviour of plane concrete in a stable manner, the strain vector on the micro-plane (Figure 1) is the projection of the macroscopic strain tensor

_{i}*ε*

*. So the components of this vector are*

_{ij}*ε*

*=*

_{Ni}*n*

_{i}*ε*

*that is*

_{Ni}

where repeated indices imply summation over *i* = 1, 2, 3. The mean normal strain, called the volumetric strain *ε** _{v}* and the deviatoric strain

*ε*

*on the micro-plane can also be introduced, defined as follows:*

_{D}*ε*

*and*

_{v}*ε*

*is useful when the effect of the hydrostatic pressure on a number of cohesive frictional materials, such as concrete, needs to be captured.*

_{D}To characterise the shear strains on the micro-plane (Figure 1), we need to define two co-ordinate directions, *M* and *L*, given by two orthogonal unit co-ordinate vectors **m** and **l** of components *m _{i}* and

*l*lying on the micro-plane. To minimise the directional bias of

_{i}**m**and

**l**among micro-planes, one of the unit vectors

**m**and

**l**tangential to the plane is considered to be horizontal (parallel to

*x - y*plane).

The magnitudes of the shear strain components on the micro-plane in the direction of **m** and **l** are *ε** _{M}* =

*m*(

_{i}*ε*

*),*

_{ij}n_{j}*ε*=

_{L}*l*(

_{i}*ε*

*). Because of the symmetry of tensor*

_{ij}n_{j}*ε*, the shear strain components may be written as follows (e.g. Bažant et al 1996):

_{ij}in which the following symmetry tensors were introduced:

Once the strain components on each micro-plane have been obtained, the stress components are updated through micro-plane constitutive laws, which can be expressed in algebraic or differential forms.

In the kinematic constraint micro-plane models, the stress components on the micro-planes are equal to the projections of the macroscopic stress tensor *σ*_{ij} only in some particular cases when the micro-plane constitutive laws are specifically prescribed in a manner such that this condition can be satisfied. This happens for example in the case of elastic laws at the micro-plane level, defined with elastic constants chosen so that the overall macroscopic behaviour is the usual elastic behaviour (Carol & Bažant 1997; see also Carol et al 1992, 2001). In general, the stress components determined independently on the various micro-planes will not be related to one another in such a manner that they can be considered as projections of a macroscopic stress tensor. Thus the static equivalence or equilibrium between the micro-level stress components and macro-level stress tensor must be enforced by other means. This can be accomplished by application of the principle of virtual work, yielding:

where Ω is the surface of a unit hemisphere; *σ*_{v} and *σ*_{D} are the volumetric and deviatoric parts of the normal stress component; and *σ*_{L }and *σ*_{M} are shear stress components on the micro-planes respectively. Eq (5) is based on the equality of the virtual work inside a unit sphere and on its surface, rigorously justified by Bažant et al (1996).

This numerical integration technique for evaluation of the integral statement in Eq (5) yields:

*N _{m}* is the number of integration points on the hemisphere. Based on this formulation, the following macroscopic constitutive matrix in the proposed model is obtained as follows:

where *E* and *v* are the elastic modulus and Poisson's coefficient.

**FORMULATION OF NEW ANISOTROPY DAMAGE FUNCTION**

The deviatoric part of the constitutive matrix is computed from superposition of its counterparts on the micro-planes. Such counterparts in turn are calculated according to the types of damage that occur on each plane, depending on its specific loading conditions. This damage is evaluated on the basis of the five separate damage functions. These five loading conditions are:

hydrostatic compression

hydrostatic extension

pure shear

shear + compression

shear + extension

A specific damage function is assigned according to the authoritative laboratory test results available in the literature. Then, for each state of plane loading, one of the five introduced damage functions is computed with respect to the history of micro-stress/strain components. The five damage functions are as follows:

]]>Parameters *a* to *k* in the above relations are computed according to the laboratory test results obtained for each specific concrete. In Eq (9), *ε** _{eq}* is an average strain, and in the other relations it is the magnitude of the projected deviatoric strain vector on each micro-plane. To obtain the parameters

*a*to

*k*, uniaxial compression tests are used to determine the planes' shear-compression; therefore, parameters

*d, e, f, g, h, i, j*and

*k*are found by numerical trial and error. Adjusting the stress path in the triaxial test to provide a plane with pure tension or pure shear can lead to the determination of

*a, b*and

*c*.

In this formulation, for simplicity, we consider just two basic material parameters, namely the elastic modulus and Poisson's ratio. Accordingly, a simple isotropic linear elasticity governs the intact material and all non-linearity, including damage effects, is accounted for through damage functions on the sampling planes. Figure 3 shows the computational sequence used in the proposed model.

**CORRELATION STUDIES OF THE UNIAXIAL COMPRESSION, TRIAXIAL AND UNIAXIAL TENSION TESTS**

To demonstrate the validity of the proposed model, correlation studies were done between the model's analytical results and experimental evidence from tests regarding the stress-strain response of concrete specimens under different loading conditions. The tests used were the uniaxial compression (UC) test (Figures 4 and 5), the conventional (cyclic) triaxial compression (CTC) test and the uniaxial tension (UT) test. The proposed model should be able to show the activity of each predefined micro-plane up to and even after failure. This specific novel aspect has been evaluated for each test.

]]>

The material parameters used in the above analysis are: E = 25 000 MPa and *v* = 0,20.

**Uniaxial compression (UC) testing**

In Figure 5, the volumetric changes of the concrete specimen under uniaxial compressive loading are shown to compare well with the experimental results of Kupfer and coworkers (1969).

]]> Figure 6 compares the variation in micro-stress. As the figure shows, during the application of the uniaxial compressive load, micro-plane number 11 (see Figure 2) is in compressive stress, whereas micro-planes numbers 9, 10, 12 and 13, which are orientated normal to the load direction on the unit sphere, experience only tensile stress. Compressive stress accompanied by shear strains affects the remaining planes. It is noted that during increments of the uniaxial compressive load, the compressive and shear stress components acting on micro-planes numbers 1 to 8 increase simultaneously, with a greater rise of shear stress at first; however, near to the peak stress () the compressive stress decreases suddenly.

Figure 7 shows the growth of the damage function values of different micro-planes during the uniaxial compression test of concrete obtained with the proposed model. As can be observed from this figure, damage evolves faster on micro-planes numbers 9, 10, 12 and 13 on the unit sphere than on the other planes. This is because there are different modes of loading on those planes. In the uniaxial compression test done by the proposed model, the axial compressive load is applied on the x-axis that is normal to micro-plane number 11 (Figure 2). So on this plane there is only a normal compressive load (mode *I*) which means that no damages could be caused to it. On micro-planes numbers 5, 6, 7 and 8 there is a shear stress combined with the normal compressive load (mode *IV*); this causes less damage than on micro-planes 1, 2, 3 and 4 on which there is the same mode of loading (mode *III*). This is because on micro-planes numbers 1, 2, 3 and 4 the magnitude of the compressive stress component is less than that on micro-planes numbers 5 to 8 (see Figure 6), so damage grows faster. Finally, on micro-planes numbers 9, 10, 12 and 13 there is only normal tension (mode *II*), which causes damage to grow faster than on all the other planes.

]]> If the points on the top and bottom loading plates are assumed not to be horizontally constant, the priority of micro-plane activity in the above test is the other way round. This means that damage on micro-planes numbers 1, 2, 3 and 4 will be greater and cracks will be initiated first on these micro-planes. This phenomenon is in agreement with the model's predictions and is depicted in Figure 8.

In this test, the hydrostatic pressure is first applied to the specimen up to a certain level and then the axial compression is increased while the lateral or confining pressure is held constant. Therefore, in this test, up to a certain level of hydrostatic pressure there should be no shear forces on the micro-planes. This can be seen in Figure 9 which shows the evolution of the micro-stress components on the different micro-planes during the CTC test.

]]> Finally, the effect of lateral confining pressures on the compressive cylindrical strength of concrete specimens simulated by the proposed model is compared with the experimental data of Ansari and Li (1998) in Figure 10. The stress-strain response of the concrete cylindrical specimen under axial tension load is depicted in Figure 11.

**Cyclic triaxial compression (CTC) testing**

In Figure 12, the predicted response of the model under the cyclic compression test is compared with the experimental results of Sinha et al (1964).

In addition, Figure 13 shows the behaviour of concrete simulated by the proposed model under complete cyclic triaxial loading.

]]>

**Finite element code**

Three-dimensional (3D) finite element software was developed using the proposed damage model for concrete as an alternative to the non-linear material model. This new software can perform 3D crack/damage analysis of concrete structures under static/dynamic and monotonic/cyclic loadings. The operating sequence of the proposed code for non-linear static analysis is illustrated in Figure 14.

**Three-point bending test**

Figure 15 shows the geometry and finite element meshes for this test. The tensile strength of concrete and the elastic modulus are assumed to be = 2,8 MPa and *E* = 27 500 MPa. For modelling, 428 3D eight-noded brick elements and 958 nodes are used. The path of crack growth simulated with the proposed model is illustrated in Figure 16. As can be seen from this figure, the damaged region of the specimen coincides with experimental observations.

Figure 17 shows a comparison of the load-displacement curve obtained with the proposed model with that found by Bažant and Luzio (2004). Clearly, there is reasonable agreement between the obtained and previously reported results. This example shows the mode *I* fracture of concrete in this 3D condition. Furthermore, because of the specific formulation of damage on micro-planes, the suggested model is able to predict the direction of the crack line at each integration point.

]]>

This can be observed in Figure 18, in which the damaged micro-planes at integration point number 681 are illustrated (orange planes). On these planes there is only tensile strain (load condition *II*) and mode *I* fracture has occurred on them. A comparison between the damage values obtained with the proposed code on various micro-planes at integration point 681 is shown in Figure 19. Clearly, in the case of non-symmetrical conditions, different stress/strain paths appear on the sampling planes and the arrangement of the active planes leads to a different failure configuration from that presented in the next section.

]]>

**Anti-symmetrical four-point shear test**

To show how the model is also able to predict a crack with a curved path, a notched specimen, subjected to an anti-symmetrical four-point shear loading, is considered. The thickness of the beam is 38 mm. A mesh of 428 3D brick elements is used (Figure 20).

The material parameters are *E* = 27 500 MPa and *v* = 0,20. When the anti-symmetrical load is applied, a crack starts from the left-hand corner of the notch and grows upwards to the left-hand side of the loading plate. Figure 21 shows the crack path in the specimen obtained from the proposed model and the actual monitored crack path obtained from experiment. The experimentally obtained crack ends in the left-hand corner of the loading test plate, while the predicted one ends in the right-hand corner of the load plate. There are two possible reasons for this difference. The first is that point loads are applied in the model, whereas plate loading is applied in the test. The second is that there is a small mismatch between the sampling plane orientations and those of the real cracked planes. Increasing the number of sampling planes can sort out this mismatching to some extent. Furthermore, if other element types are chosen, or a denser finite element mesh is used, the predicted crack pattern may be unchanged from the prediction with the coarser mesh of linear brick elements. Obviously, as Figure 21 shows, the proposed model is able to simulate the crack line fairly well. Figure 22 shows the modelled crack evolution in the direct tensile test.

]]>

Mode *II* fracture has occurred on microplanes numbers 5, 6, 7 and 8 (blue planes in Figure 18) due to a combination of tensile and shear stresses (load condition *V*) at the damage integration points.

**ARCH DAM CRACK ANALYSIS**

To demonstrate further the ability of the proposed model, a double curvature arch dam was selected. This dam is in the feasibility stage of consulting studies in Dez-Ab consulting engineering company and is located in the south-west of Iran near Aligoodarz city.

]]>**Dam geometry**

The geometry of the dam, including the plan and the profiles of the crown cantilever and centre lines are shown in Figures 23 and 24. The dam has a height of 220 m and the length of its crown is about 430 m.

]]>

**Material library**

The elastic modulus and Poisson's ratio of plain concrete for the dam are considered to be *E _{c}*= 30 000 MPa and

*v*= 0,20. The mass density of the concrete is

_{c}*p*= 2 400 kg/m

^{3}. The elastic modulus of the rock mass is selected as

*E*

_{c}= 20 000 MPa.

**Finite element mesh, results and failure mechanism**

The model discretisation of the arch dam used in this analysis is shown in Figure 24. A total of 652 3D brick elements and 1 146 nodes are used for the finite element calculations. The proposed arch dam is analysed by the code in Figure 14 under its own weight and the hydrostatic pressure of the reservoir in the first stage. Damaged portions of the dam are shown in Figure 25 in terms of the damage parameters (damage value *w* 1,0) of the proposed model. It can be observed that the damage or cracks at all designated integration points in Figure 25 occurs on micro-planes numbers 7 and 8 where tensile and shear stress components (mode *II* fracture) exist. This indicates that the combination of tensile and shear stresses (load condition *V*) have the greatest damage effects on the arch dam under its own weight and the pressure of the reservoir at that stage.

The dam's deformations and stresses under its own weight and normal water pressure are shown in Figure 26. The stress paths of the tensile and shear components on micro-planes numbers 7 and 8 due to the initial loading up to the point of failure are shown in Figure 27.

]]>

The six stages of progressive failure are shown in Figure 28. Damage starts in the bottom right-hand location under the dam's own weight and water pressure increases. With the compliance matrices kept constant for the damaged planes, failure extends to the left and upwards to configure a failure mechanism that finally divides the dam body into five pieces (Figure 29).

]]>

**CONCLUSIONS**

A constitutive damage model for concrete predicting the effects of any arbitrary loading conditions was developed using the approach of a theoretical framework of micro-planes and damage functions. This model has some characteristics that distinguish it from similar work done recently in this area of research. Some of those differences are listed below.

A new damage formulation has been employed in the micro-plane model. It was built on the basis of five fundamental force conditions that essentially can occur on each micro-plane. Consequently, any change in six strain/stress components in a

dxdydzelement leads to one of the five stress cases introduced into plane conditions and the use of corresponding damage functions. Therefore, the proposed model is capable of predicting the concrete behaviour under any arbitrary strain/stress path. These five force conditions are: () hydrostatic compression, (I) hydrostatic extension, (II) pure shear, (III) shear + compression, and (IV) shear + extension.VFive damage evolution functions (all functions of equivalent strain) were formulated for any of five force conditions. The equivalent strain for the first two conditions is defined as volumetric strain and for the others it may be considered as the magnitude of the projected vectors of the deviatoric strain tensor on the corresponding micro-planes.

These damage evolution functions were constructed with reference to the experimental test results for concrete specimens under compressive and tensile loading conditions reported by the various researchers.

In our formulation, according to the damage theory, the value of total damage functions in each micro-plane varies between zero and one, the first value being related to the undamaged state (no crack initiation on the micro-plane) and the second value to the totally damaged state (where the crack opening on the micro-plane is greater than the specified critical value).

As soon as any crack initiates and starts to grow on each of 26 micro-planes, the stiffness of the plane orthogonal to the crack line gradually decreases. When the crack opening reaches to its critical value, the components of the stiffness matrix related to the direction of normal to the crack are reduced to zero. If the on-plane loading conditions change in such a way that the crack starts to close and the crack opening then reaches a value smaller than the critical crack opening, these stiffness components increase again and return to their initial values.

]]> It was observed that for a particular micro-plane, the combination of shear and extension strain components (mode) could result in much more damage than the other situations. Thereafter, hydrostatic extension (V), pure shear (II), shear + compression (III) and hydrostatic compression (IV) respectively have a smaller effect on damage evolution.IThe novel micro-plane damage model is able to simulate the behaviour of concrete specimens under both compressive and tensile loadings with few model parameter requirements.

Although the proposed model has excellent features, such as pre-failure configuration of the internal material, prediction of the final failure mechanism, the capability of seeing induced/inherent anisotropy and also any fabric effects on material behaviour, the basis of its formulation is simple and logical, and it confers some physical insights that make it convenient to understand.

**REFERENCES**

Ansari, F & Li, Q B 1998. *Fiber Optic Sensors for Construction Materials and Bridges*. Lancaster, PA: Technomic Publishing Co., p 368. [ Links ]

Batdorf, S B & Budiansky, B 1949. A mathematical theory of plasticity based on the concept of slip. Technical Note 1871, National Advisory Committee for Aeronautics. [ Links ]

Bažant Z P & Luzio G D 2004. Nonlocal microplane model with strain-softening yield limits. *International Journal of Solids and Structures*, 41: 7209-7240. [ Links ]

Bažant, Z, Xiang, Y & Prat, P C 1996. Microplane model for concrete. I. Stress-strain boundaries and finite strain. American Society of Civil Engineers, Journal of Mechanical Engineering, 122(3): 245-254. [ Links ]

Bažant, Z P & Gambarova, PG 1984. Crack shear in concrete: crack band micro plane model. *Journal of Structural Engineering*, ASCE, 110: 2015-2036. [ Links ]

Carol, I, Jirasek, M & Bažant, Z P 2001. A thermodynamically consistent approach to micro plane theory. Part I: Free energy and consistent micro plane stresses. *International Journal of Solids and Structures*, 38: 2921-2931. [ Links ]

Carol, I & Bažant, Z 1997. Damage and plasticity in micro plane theory. *International Journal of Solids & Structures*, 34: 3807-3835. [ Links ]

Carol, I, Bažant, Z P & Prat, P. 1992. New explicit micro-plane model for concrete: theoretical aspects and numerical implementation. *International Journal of Solids and Structures*, 29: 1173-1191. [ Links ]

Kupfer, H B, Hilsdorf, H K & Rusch, H 1969. Behaviour of concrete under biaxial stresses. *Journal of the American Concrete Institute*, 66(8): 656-666. [ Links ]

Labibzadeh, M & Sadrnejad, S A 2008. Mesoscopic damage-based model for plane concrete under static and dynamic loadings. *American Journal of Applied Science*, 3(9): 2011-2019. [ Links ]

Pande, G N & Sharma, K G 1983. Multilaminate model of clays - A numerical evaluation of the influence of rotation of principal stress axes. *International Journal of Numerical and Analytical Methods in Geomechanics*, 7: 397-418. [ Links ]

Sadrnejad, S A & Pande G N 1989. A multi-laminate model for sand. *Proceedings,* 3rd International Symposium on Numerical Models in Geomechanics, NUMOG-III, Niagara Falls, Canada, 8-11 May 1989. [ Links ]

Sadrnezhad, S A 1992. Multi-laminate elasto-plastic model for granular media. *Journal of Engineering*, (Iran), 5(1 & 2): 11. [ Links ]

Sinha, B P, Gerstle, K H & Tulin, L G 1964. Stress-strain relations for concrete under cyclic loading. *Journal of the American Concrete Institute*, 61(2): 195-211. [ Links ]

Taylor, G I 1938. Plastic strain in metals. *Journal of the Institute of Metals*, 62: 307-324. [ Links ]

Zienkiewicz, O C, & Pande, G N 1977. Time-dependent multi-laminate model of rocks. *International Journal of Numerical and Analytical Methods in Geomechanics*, 1: 219-247. [ Links ]

**Contact details:** ]]>
Department of Civil Engineering KN Toosi University of Technology

Tehran Iran

Tel: 98 21 8888 1128 Fax: 98 21 8877 9385

e-Mail: sadrnejad@kntu.ac.ir

SEYED AMIRODIN SADRNEJAD received his BSc from Sharif University of Technology, Teheran, his MSc from University College Cardiff, UK, and his PhD from University College Swansea, UK. He is Professor in Civil Engineering at the K N Toosi University of Technology and vice-principal of the university. He is a consultant to Teheran Municipality on design and construction activities for residential and industrial construction. His fields of interest include earth and concrete dams, numerical methods in geomechanics and structural engineering, finite element methods, hydromechanics and hydraulic structures. |