Engineering Structures 115 (2016) 155–164

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Adaptive mesh refinement FEM for seismic damage evolution in concrete-based structures Bin Sun, Zhaoxia Li ⇑ Department of Engineering Mechanics, Jiangsu Key Laboratory of Engineering Mechanics, Southeast University, Nanjing 210096, China

a r t i c l e

i n f o

Article history: Received 30 June 2014 Revised 11 February 2016 Accepted 16 February 2016

Keywords: Damage Concrete-based structures Multi-grid FEM Reinforced concrete column

a b s t r a c t In order to study the seismic damage mechanisms of concrete-based structures, a adaptive multi-grid method is developed to simulate the process from evolving damage to structural failure in concrete structures under seismic loading. The theory of multi-grid FEM coupled evolving damage is developed on the basis of the improved variational principle, in which the elements in each sub-domain with different grid sizes are under the different state of damage. Meanwhile the FE elements can be refined adaptively with the accumulation of its damage for lower cost, without user intervention in the computation. As a case study of the method, the failure process of reinforced concrete column under seismic loading is simulated, the results show that, the simulated results fit well with the experiment, and the developed method can be used to reveal the seismic damage mechanism of concrete-based structures by considering the process from material damage to structural failure; the adaptive multi-grid method can give satisfactory result with lower cost by compared with experiment. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Due to its good qualities and low cost, concrete is widely used in civil engineering structures [1–3], e.g. concrete buildings and bridges, which exist extensively in the world. In case of the concrete structures failure during the service period, such as suffering earthquake, it will inevitably lead to disaster results. So in order to avoid catastrophe, study on seismic damage of the concrete-based structures to assess their seismic performance has become an important issue nowadays [4–6]. One of the characteristics of the failure of structure, especially for the structure of brittle material, is the sample-specificity [7–9], namely there are significant disparities at failure from sample to sample, even though they are initially identical at the macro-scale. It is because that the effect of material damage on global behavior of structure can be amplified strongly due to sensitive trans-scale coupling [9]. Therefore, in order to study well the failure mechanism of structure, it is necessary to simulate the coupling process of material damage up to structure failure. In the failure process of structures, the initiation and growth of material damage is an important phase. Damage is a physical process which leads to the progressive deterioration of structures due to initiation and growth of micro-cracks or micro-voids leading to ⇑ Corresponding author. E-mail address: [email protected] (Z. Li). http://dx.doi.org/10.1016/j.engstruct.2016.02.021 0141-0296/Ó 2016 Elsevier Ltd. All rights reserved.

the formation of local failure and structural failure for structures [10,11]. Therefore, in order to study on the seismic failure mechanism of concrete structures, the material damage of concrete should first be described. Since the beginning of the last century, the damage evolution processes in solids were attempted to be described [12]. At the beginning of the study, only few simple problems of damage can be solved without effective numerical methods. With the development of numerical methods especially for FE method (FEM), more and more damage evolution processes can be simulated numerically. Ghosh et al. have developed a voronoi cell FE method (VCFEM) to study the damage evolution in composites as well as heterogeneous aluminum alloys [13–17]. Fish delineated the damage phenomena in composite materials using mathematical homogenization method [18,19]. Könke has proposed a simulation method of damage evolution for ductile materials, which divides complete damage range into the micro- and macro-damage range [12,20]. Feng developed an adaptive mesh refinement FEM algorithm to demonstrate the damage evolution of brittle materials [9]. The methods for describing damage evolution can be mainly classified into two classes. The first case implements meso-structure of material in computing model to describe the damage evolution [13–19], while the second case uses a damage variable to describe the deterioration process of materials based on continuum damage mechanics (CDM) ignoring the meso-structure of material [9,12,20]. For the first case, the focus of attention is the damage evolution of materials and the scales

156

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

of the research objects are significantly smaller than the structural dimensions under the limitation of computing power. For the large scale engineering structures, CDM is an effective numerical tool [21]. Although this may decrease the accuracy due to ignoring the realistic meso-structure and micro-damage (e.g., micro-cracks and micro-voids), it brings the convenience of calculation. For concrete structures, CDM is also an effective numerical tool with wide applications to describe the damage evolution of concrete [22–24]. Most previous numerical methods use a damage index to evaluate the damage state and predict the structural failure of concrete structures based on the macroscopic response of structure ignoring the coupling process from material damage to local damage and failure and eventually to global structural failure [25–27]. Although some theories and computational frameworks were attempted to be developed [28,29], the validity of these frameworks had not been proved by experiment. In order to obtain an accurate estimation of damage state of concrete structures under seismic loading, it is necessary to develop a new method to simulate the process from evolving damage to structural failure, which is the basic aim of the paper. As a case of study of the proposed method, simulation of the process from evolving damage to failure of a concrete column under cyclic loading is given and compared with experiment. 2. FEM formulation of the coupling multi-grid interface and evolving damage

tensor and strain tensor in Xk , Dk is damage variable in Xk , Sr is the i . boundary surface subjected to the action of external surface force p

P12 ¼ P23 ¼ P31 ¼

e ð1AÞ 1 p e exp BðAee Þ ½ p

where ters,

0

e 6 ep e > ep

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ P3 jei jþei ; i¼1 hei i hei i ¼ 2

@ X2 \@ X3

ð1Þ

X1

@ X3 \@ X1

@ X1 \@ X2

ð6Þ

#

de 1 ij

@ X1 \@ X2

Pk þ P12 þ P23 þ P31

@ X3 \@ X1

where Pk is the functional of sub-domain Xk (k ¼ 1; 2; 3), P12 is the functional of interface between the boundary of X1 and the boundary of X2 ð@ X1 \ @ X2 Þ, P23 is the functional of interface between the boundary of X2 and the boundary of X3 ð@ X2 \ @ X3 Þ, P31 is the functional of interface between the boundary of X3 and the boundary of X1 ð@ X1 \ @ X3 Þ. The coupling of the sub-domains with different grid sizes can be obtained using the Lagrange multipliers. And the different terms can be rewritten as:

Z

i du1i dS p

1 k31 i dui dS ¼ 0

1 2 dk12 i ðui ui ÞdS ¼ 0

@e

2 ij

Z

@ X1 \@ X2

ð8Þ

#

Z

de2ij F i du2i dX2

Sr \X2

Z 2 k12 i dui dS þ

@ X2 \@ X3

i du2i dS p

2 k23 i dui dS ¼ 0

@ X2 \@ X3

2 3 dk23 i ðui ui ÞdS ¼ 0

Equilibrium in X3 :

" @Aðe3ij ; D3 Þ

@e

3 ij

Z

@ X2 \@ X3

#

Z

de3ij F i du3i dX3 Z

3 k23 i dui dS þ

@ X3 \@ X1

Sr \X3

i du3i dS p

3 k31 i dui dS ¼ 0

@ X3 \@ X1

3 1 dk31 i ðui ui ÞdS ¼ 0

Aðeij ; DÞ ¼

1 1 rij eij ¼ aijkl eij ekl ð1 DÞ 2 2

ð13Þ

where rij ¼ aijkl ekl ð1 DÞ is the Cauchy stress tensor [21], aijkl is the initial elastic stiffness tensor. According to the concept of effective ~ ij can be expressed as [21]: stress, the effective stress tensor r

r~ ij ¼

rij

1D

ð14Þ

In Eqs. (7), (9), (11), we use the relations:

where damage is considered to be related to strain.

ekij are respectively the displacement

ð12Þ

The body of concrete can be approximately considered as elastic body [9]. And the strain energy density coupling damage of elastic body can be written as [21]:

distributed body force, uki and

Sr \Xk

ð3Þ

ð11Þ

Equilibrium on @ X3 \ @ X1 :

where Aðekij ; Dk Þ is the strain energy density, F i ði ¼ 1; 2; 3Þ is

Xk

i uki dS ðk ¼ 1; 2; 3Þ p

ð9Þ

ð10Þ

1 2 @Aðeij ; DÞ 1 @D ~ ij eij ¼ rij r @ eij 2 @ eij

Pk ¼

ð7Þ

Equilibrium on @ X2 \ @ X3 :

Z

Z

ð2Þ

Sr \X1

Z 1 k12 i dui dS

" @Aðe2ij ; D2 Þ

X2

Z

Z

dX1

F i du1i

Equilibrium on @ X1 \ @ X2 :

Z

k¼1

i Aðekij ; Dk Þ F i uki dXk

3 1 k31 i ðui ui ÞdS

@ e1ij

Z

þ

X3

3 X

h

ð5Þ

Z

" @Aðe1ij ; D1 Þ

Z

ei is the principal strain :

In order to obtain satisfactory result with lower cost, computational domain with different grid sizes can be adaptively generated and changed with the damage evolution in the developed method. As shown in Fig. 1, there are three sub-domains with different grid sizes, which respectively are X1 ; X2 and X3 : And the sum of functionals of all sub-domains and interfaces P can be written as:

Z

2 3 k23 i ðui ui ÞdS

u1i ; u2i ; u3i are respectively the displacement tensor on the boundary of X1 ; X2 and X3 : Using the variational principle dP ¼ 0; the weak formulation of the equations can be expressed as: Equilibrium in X1 :

2.2. Variational formulation

P¼

ð4Þ

Z

X2 ; k23 i is Lagrange multipliers tensor which connect the X2 and X3 ; k31 is Lagrange multipliers tensor which connect the X3 and X1 ; i

Z

ep is the strain threshold of damage, A; B are model parame-

e¼

1 2 k12 i ðui ui ÞdS

where k12 i is Lagrange multipliers tensor which connect the X1 and

Based on CDM, the famous Mazars concrete damage model is chosen for describing the damage of concrete in this paper. And the Mazars concrete damage model can be written as follows [30]:

D¼

@ X1 \@ X2

Equilibrium in X2 :

2.1. Damage model for concrete

(

Z

eij ¼ ðui;j þ uj;i Þ

ð15Þ ð16Þ

157

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

Ω2

Ω1

Using the developed method to simulate the damage and failure process of structures

Ω3

Ω=Ω1 +Ω 2 +Ω3

Ω=Ω1 (a) Initial FE model consists only of coarse grid before the operation of the method

(b) FE model consists of multi-grid after the operation of the method Fig. 1. Adaptive split criterion.

2.3. Multi-grid FEM equations coupling multi-grid interface and damage The displacement uki ðk ¼ 1; 2; 3Þ in Xk ; and the Lagrange multipliers 23 31 k12 i ; ki ; ki

are interpolated from nodal values using shape functions as:

k

uk ¼ N uk

ðk ¼ 1; 2; 3Þ

ð17Þ

k12 ¼ N12 k K12

ð18Þ

N23 k K23 N31 k K31

ð19Þ

23

k

¼

k31 ¼

ð20Þ

Eqs. (7)–(12) can be rewritten as:

where L is the operator matrix transforming displacements into strains, which the has the similar form to that in the undamaged strain–displacement for small deformation problems, I is identity matrix, E is initial elastic stiffness matrix for undamaged material, is the prescribed surface F is the prescribed body force matrix, p force matrix, the matrix Mk ðk ¼ 1; 2; 3Þ is defined as:

Mkij ¼

@Dk k e @ ekij ij

The Newton–Raphson iterative method is chosen to solve the non-linear Eqs. (21)–(26), which requires the linearization of each equation, according to the linearized form:

Z

R X1 ¼

n o T T ðLN1 Þ ðI D1 ÞELN1 ðM1 LN1 Þ ELN1 u1 dX1 X1 Z Z dS N1T FdX1 N1T p X1 Sr \X1 Z Z N1T N12 N1T N31 þ k K12 dS k K31 dS ¼ 0 @ X1 \@ X2

Z RX1 \X2 ¼ R X2

@ X1 \@ X2

@ X3 \@ X1

N12T N1 u1 N2 u2 dS ¼ 0 k

n o T T ¼ X2 ðLN2 Þ ðI D2 ÞELN2 ðM2 LN2 Þ ELN2 u2 dX2 R R dS X2 N2T FdX2 Sr \X2 N2T p R R 2T 12 @ X1 \@ X2 N Nk K12 dS þ @ X2 \@ X3 N2T N23 k K23 dS ¼ 0

ð27Þ

n Rnþ1 ðÞ ¼ R ðÞ þ

@RnðÞ @bi

Dbi

ð : X1 ; X1 \ X2 ; X2 ; X2 \ X3 ; X3 ; X3 \ X1 Þ ð28Þ

ð21Þ

where ðnÞ denotes the value at nth iteration, bi are the independent variables. Using the Newton–Raphson iterative method, Eqs. (21)– (26) can be expressed as:

2 ð22Þ

R

ð23Þ

K1

6 1T 6K 6 12 6 60 6 6 60 6 60 4 K1T 31

K112

0

0

0

K131

0

K2T 12 0

0

0

K212

K2

K223

0

0

0

K2T 23

0

K3T 23 0

0

0

K323 K3

0

0

0

K3T 31

K331 0

@ X2 \@ X3

2 3 N23T k ðN u2 N u3 ÞdS ¼ 0

ð24Þ

Z

R X3

n o T T ¼ ðLN3 Þ ðI D3 ÞELN3 ðM3 LN3 Þ ELN3 u3 dX3 X3 Z Z dS N3T FdX3 N3T p X3 Sr \X3 Z Z N3T N23 N3T N31 k K23 dS þ k K31 dS ¼ 0 @ X2 \@ X3

@ X3 \@ X1

Z RX3 \X1 ¼

@ X3 \@ X1

7 7 7 7 7 7 7 7 7 7 5

3 3ðnÞ 2 RX1 Du1 ðnÞ 6 DK 7 6 RX \X 7 6 12 7 6 1 27 7 7 6 6 6 Du 2 7 6 RX2 7 7 7 6 6 7 6 DK 7 ¼ 6 R X2 \X3 7 6 23 7 6 7 7 6 6 4 Du 3 5 4 RX3 5 2

DK31

RX3 \X1 ð29Þ

Z RX2 \X3 ¼

3ðnÞ

N31T k

N3 u3 N1 u1 dS ¼ 0

where

Kk ¼

Z n o T T ðLNk Þ ðI Dk ÞELNk ðMk LNk Þ ELNk dXk ðk ¼ 1;2;3Þ ð30Þ Xk

Z Kk12 ¼

@ X1 \@ X2

NkT N12 k dS ðk ¼ 1;2;3Þ

ð31Þ

NkT N23 k dS ðk ¼ 1;2;3Þ

ð32Þ

NkT N31 k dS ðk ¼ 1;2;3Þ

ð33Þ

Z ð25Þ

Kk23 ¼

ð26Þ

Kk31 ¼

@ X2 \@ X3

Z

@ X3 \@ X1

158

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

Repeat the iterative process until nth the iteration until the convergence criterion is met:

2 R 2 3 R 3ðnÞ 1T 1T RX dS X1 N FdX1 þ Sr \X1 N p 1 6 6 7 6 6 RX1 \X2 7 7 0 7 6 7 6 , 7 6 R N2T FdX þ R 6 R dS 7 7 N2T p X2 6 X2 6 2 7 Sr \X2 7 6 6 7 6 tol 6 RX2 \X3 7 7 6 0 7 6 R 6 7 7 R 6 4 R 7 3T 3T 5 dS 5 X3 4 X3 N FdX3 þ Sr \X3 N p RX3 \X1 0 2

2

ð34Þ where tol is tolerance of equilibrium.

adaptively refined into 8 elements consisting of X3 , once its damage value reaches D2c . Third, a global damage index of structures is developed to evaluate the global seismic damage state of the whole structures due to material damage evolution and local failure. And the damage index reflects the degradation of the structural behavior. When structure loses bearing capacity, it is considered to be failure. For example, with accumulation of the material damage and local failure area, the force of the load point is close to 0 even if the displacement loading acting on the load point is very large, while the displacement of the load point is very large even if the force loading acting on the load point is very small. So the global damage index Dglobal is established based on the load–displacement curve of structures:

Dglobal ¼ 1 3. Computational strategy for adaptive mesh refinement FEM coupling evolving damage 3.1. Description of the damage and failure process In order to study well damage and failure mechanism of concrete structures, it is necessary to simulate the process from damage evolution to structural failure. Depending on the purpose, there are three main processes to be simulated in the developed method. First, the process of material damage evolution should be delineated in the developed method. Here famous Mazars concrete damage model based on CDM is chosen and coupled in the stress–strain relationship to describe the initiation and growth of micro-cracks or micro-voids using a damage variable. Second, the process from the material damage to local failure of structure should be delineated in the developed method. When the damage of element reaches 1, the element is considered to be failure and deleted. Meanwhile, in order to obtain the relatively accurate initiation and growth position of the local failure in structure and enhance the accuracy of finite element solutions of problems with high gradients due to damage, FE meshes consisting of a lot of elements are necessary when using traditional methods for the whole small-scale simulation. This often leads to consume a lot of computation time and even make the analysis impossible. Therefore, in this paper a 3D adaptive mesh refinement procedure is developed and implemented in the method for sufficient precision and lower cost. By using the method coupled with the procedure, the elements can be refined to smaller scales adaptively, once their damage values reach preset value. As shown in Fig. 1, One element in X1 can be adaptively refined into 8 elements consisting of X2 , once its damage value reaches D1c . One element in X2 can be

ΔlN

Δl0 F0

(a) Initial undamaged stage

FN Dl N

F0 Dl 0

ð35Þ

where F N ; DlN are respectively the load and displacement of the load point of the structure at Nth loading, F 0 ; Dl0 are respectively the initial load and displacement of the load point, when the global structure is in elastic stage, as shown in Fig. 2. And a procedure for monitoring Dglobal is implemented in the method. Once the Dglobal reaches 1, the concrete column is considered to be failure. 3.2. The basic procedure of the developed method As described the strategy for the simulation of the process from material damage to structural failure in Section 3.1, the flow chart of the developed method is given in Fig. 3, which chooses route 1 after the position A. 3.3. Technique for overcoming ‘‘island effect” After the operation of the algorithm, which choose route 1 after the position A in Fig. 3, some initial coarse elements are refined to smaller scale elements and some elements are failure, whose damage reach the preset critical value Dc, as shown in Fig. 4a. Note that in Fig. 4a there appear some connect regions consisting of elements suspending in the air, which does not agree with the fact due to gravity. The cause of this condition is that the damage of all elements around the connect region reach 1. The connect region suspending in the air look like computational debris which are needed to be eliminated automatically in the developed method. We define the formation of in the connect region suspending in the air as ‘‘island effect”, and the element in the region as ‘‘island element”. To eliminate the island elements, a procedure is developed and implemented in the method to overcome island effect. And the flow chart of the procedure is in the dashed boxes

Damage

FN

(b) Damaged stage at Nth loading Fig. 2. Definition of global damage index.

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

159

Fig. 3. Flow chart of the method.

of Fig. 3. After the operation of the algorithm, which choose route 2 after the position A in Fig. 3, the island elements can be automatically deleted, as shown in Fig. 4b. 4. Numerical analyses on seismic damage evolution of a reinforced concrete column

the numerical instability. And the ultimate strength of steel reinforcement bar can be founded in the experimental reports of the reinforced concrete structure under seismic loading [31]. Based on this way, steel reinforcement bar can be considered in the simulation of damage evolution of the reinforced concrete column. 4.2. Experimental validation

4.1. Description of steel reinforcement bar in the model of the reinforced concrete column Since steel reinforcement bar is also very important in the reinforced concrete structures except concrete, steel reinforcement bar should be considered in the simulation. And in order to consider the effect of steel reinforcement bar in the reinforced concrete column, which provides tensile strength, discrete steel reinforcement elements are coupled in the FE model as shown in Fig. 5b. Once steel reinforcement bar element reaches to the ultimate strength, the element is considered to be failure for preventing

Before using the developed method, parameters for material property of concrete should first be determined. And here the model parameters can be obtained as following way. Based on the given standard intervals for the model parameters of concrete in Ref. [32], adjusting macro-model parameters of concrete makes that the macroscopic analysis is in good accordance with experimental results, such as stress–strain curves. Based on this way, which is similar to the Ref. [1], a good set of model parameters can be obtained and the details can be found in our previous work [33].

160

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

Island elements

(b) Multi-grid of concrete which choose route 2 in the method

(a) Multi-grid of concrete which choose route 1 in the method

Fig. 4. Calculation result made up of concrete elements after the operation of the algorithm.

(b) Initial FE model

(a) Test setup of concrete column Fig. 5. Numerical example of RC column.

mm 60

Horizontal displacement loading

KN

Horizontal force loading

30 20 10 0 -10 -20

40 20 0 -20 -40

-30 0

0.02

0.04

0.06

0.08

0.1

0.12

N / Nf

0.14

0.165

-60 0.165

0.3

0.4

0.5

0.6

0.7

0.8

0.9

N / Nf

(b) Horizontal displacement loading

(a) Horizontal force loading Fig. 6. Load spectrum.

1

161

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

(a1) Calculation result at the cycle ratio

(b1) Calculation result at the cycle ratio

N / N f =29/38

N / N f =33/38

(c1) Calculation result at the cycle ratio N / N f =1

Damage

(a2) Experimental result at the cycle ratio

(b2) Experimental result at the cycle ratio

N / N f =29/38

N / N f =33/38

(c2) Experimental result at the cycle ratio N / N f =1

Fig. 7. Comparison between calculation and experimental results.

(a1) Undamaged state at N / N f =5/38

(a2) Initiation of damage at N / N f =6/19

(a3) Accumulation of damage at N / N f =23/38

(a5) Accumulation of (a4) Initiation of damage and local local failure area failure area at at N / N f =13/19

(a6) Failure of concrete column at N / N f =1

N / N f =31/38

Damage

Fig. 8. Simulation of the trans-scale failure process.

162

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164 1 0.9

In order to study the seismic damage mechanisms of concrete structures, the columns under cyclic loading is chosen as the numerical example, since its collapse experiment had been carried out by Tsinghua University [31]. As shown in Fig. 5. The dimension of the concrete column is 200 mm 200 mm 850 mm. The load spectrum is shown in Fig. 6. And N is the current number of loading cycles, the failure number of loading cycles N f is 38. As shown in Fig. 7, the failure evolution pattern predicted by simulation using the developed method agrees well with the experimental result by comparison with experiment. And the location of the initiation and growth of local failure area predicted by simulation is verified through experiment.

Calculation result Calculation result Experimental result Experimental result

0.8 0.7

Dglobal

0.6 0.5 0.4 0.3 0.2 0.1 0

4.3. Numerical simulation of the damage and failure process 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

N / Nf Fig. 9. Global damage index of the concrete column under the cyclic loading.

(a1) With scale-1 at N / N f =29/38

(b1) With scale-2 at N / N f =29/38

(c1) With scale-3 at N / N f =29/38

(a2) With scale-1 at N / N f =33/38

(b2) With scale-2 at N / N f =33/38

(c2) With scale-3 at N / N f =33/38

Compared with the previous numerical methods, one of the main features of the developed method is that the process from evolving material damage to structural failure can be simulated concurrently in a coupled manner. And the damage and failure pro-

(a3) With scale-1 at N / N f =1

(b3) With scale-2 at N / N f =1

(c3) With scale-3 at N / N f =1

Damage Fig. 10. Failure evolution patterns predicted by simulation with different grid scales.

163

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

235.1

25

106496

12.5

33.3 1.9

11.8

(a) Computational time

1664

13312

6.25

13298

(b) Total number of element

6.25

(c) Grid size of serious damaged area

Fig. 11. Comparison between the single scale simulation and adaptive multi-grid simulation.

cess of the concrete column under the cyclic loading is simulated using the developed method, as shown in Fig. 8. The process from undamaged state to initiation of material damage, to initiation and growth of local failure area due to accumulation of material damage, and eventually to structural failure, can be simulated using the developed method. Meanwhile the global damage index of the concrete column under the cyclic loading can be obtained from the monitoring data in the numerical method and experimental data based on Eq. (35), as shown in Fig. 9. For such a complex process, it is unavoidable that there exist error between the simulation and experiment. There are many causes. For example, the highly complex damage mechanisms of concrete structures cannot be completely contained in the numerical method. And the used famous concrete damage model cannot completely reflect the behavior of material damage. Similarly, it may also be due to the measuring error under the limitation of the experimental condition, or the real constraints for experimental condition cannot be completely reflected in simulation. When the global damage index reaches 1, the concrete column is considered to be failure. 5. Computational efficiency of the adaptive multi-grid method From recently researches, it shows that the calculated results of damage analysis using FEM considerably depend upon FE mesh [34,35]. In order to enhance the accuracy of finite element solutions of problems with high gradients due to damage, the refinement grid is required in the damaged area. But if using traditional methods for the whole small-scale simulation, this often leads to consume a lot of computation time and even make the analysis impossible. So a new numerical method coupled with 3D adaptive mesh refinement procedure is developed here to overcome the difficulty. 5.1. Numerical analysis by using the method in single scale In order to validate the efficiency of adaptive multi-grid method, we compare the computing cost of single scale simulations with different grid sizes and the adaptive multi-grid simulation, as shown in Fig. 10. The scales of grid with different colors in the singe scale simulation is given in Fig. 10, are the same as the scales of elements with the same color in the adaptive multi-grid simulation, as shown in Fig. 4. The element size of scale-1 is twice the size of scale-2, and the element size of scale-2 is twice the size of scale-3. And the simulation of damage and failure process using the method in single scale is given Fig. 10. As shown in Fig. 10, FE mesh size plays a role in simulating the precise initiation location of the local failure area and the pattern

for local failure area growth. The finer FE mesh, the more disordered details of local failure area initiation and growth can be simulated, including the specific position of initiation and the path of the growth, but it cost too much computational cost. 5.2. Comparison of efficiency As shown in Figs. 7 and 10, the adaptive multi-grid method can give satisfactory prediction with sufficient precision and lower cost by comparison with experiment and traditional pure small-scale numerical method, e.g., as shown in Fig. 11, the total numbers of the concrete elements of the two simulations are 13,298 and 106,496, respectively, and the computational time of the two simulations are 33.3 h and 235.1 h, respectively. Since the simulation of the process from material damage to failure of the reinforced concrete column is considerably nonlinear here, the computational time seems long here. But compared with the traditional single whole small-scale simulation method, the developed adaptive method indeed can save computational time and memory, and the developed method’s efficiency is verified. Therefore, the adaptive multi-grid FEM provides a numerical tool to simulate the process from material damage evolution to structural failure of the concrete-based structures with sufficient precision and lower cost. And the adaptive multi-grid method is much more economical in memory and time consumed for computation. 6. Conclusions The following specific conclusions can be obtained from the present study: (1) The theory and strategy for the adaptive mesh refinement FEM is proposed to simulate the process from material damage evolution to structural failure. (2) In the developed multi-grid FEM, FE elements can be refined adaptively with the damage accumulation without user intervention in computation. (3) Compared with the traditional method in single scale, the developed adaptive multi-grid method provides an excellent tool to obtain satisfactory result with sufficient precision and lower cost. (4) The developed method has been successfully applied to simulate the process from evolving material damage to global failure of a reinforced concrete column under cyclic loading. And the simulated results fit well with the experimental data. It is shown that the developed method can be used to reveal the damage and failure mechanisms of concrete-

164

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

based structures by considering the process from material damage in concrete of stress concentration zone to local failure in vulnerable area and eventually to structural failure.

Acknowledgements The work described in this paper was substantially supported by A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (No. CE02-2-46), Scientific Research Foundation of Graduate School of Southeast University (No. YBJJ1440), and National Natural Science Foundation of China (No. 11072060). The authors are very grateful to the reviewer for carefully reading the paper and for his/her comments and suggestions which have improved the paper quality. References [1] Wriggers P, Moftah SO. Mesoscale models for concrete: homogenisation and damage behaviour. Finite Elem Anal Des 2006;42:623–36. [2] Zhou XQ, Hao H. Mesoscale modelling of concrete tensile failure mechanism at high strain rates. Comput Struct 2008;86:2013–26. [3] Ghosh A, Chaudhuri P. Computational modeling of fracture in concrete using a meshfree meso–macro-multiscale method. Comput Mater Sci 2013;69: 204–15. [4] Kim TH, Lee KM, Chung YS, Shin HM. Seismic damage assessment of reinforced concrete bridge columns. Eng Struct 2005;27(4):576–92. [5] Del Gaudio C, Ricci P, Verderame GM, Manfredi G. Development and urbanscale application of a simplified method for seismic fragility assessment of RC buildings. Eng Struct 2015;91:40–57. [6] Bai JL, Ou JP. Earthquake-resistant design of buckling-restrained braced RC moment frames using performance-based plastic design method. Eng Struct 2016;107:66–79. [7] Xia MF, Song ZQ, Xu JB, Zhao KH, Bai YL. Sample-specific behavior in failure models of disordered media. Commun Theor Phys 1996;25:49–54. [8] Xia MF, Wei YJ, Ke FJ, Bai YL. Critical sensitivity and trans-scale fluctuations in catastrophic rupture. Pure Appl Geophys 2002;159:2491–509. [9] Feng R, Xia MF, Ke FJ, Bai YL. Adaptive mesh refinement FEM for damage evolution of heterogeneous brittle media. Modell Simul Mater Sci Eng 2005;13:771–82. [10] Lin J, Liu Y, Deng TA. A review on damage mechanisms, models and calibration methods under various deformation conditions. Int J Damage Mech 2005;14:299–319. [11] Arun CO, Rao BN, Srinivasan SM. Stochastic mesh free method for elastoplastic damage analysis. Comput Methods Appl Mech Eng 2010;199: 2590–606. [12] Könke C. Damage evolution in ductile materials: from micro- to macrodamage. Comput Mech 1995;15:497–510. [13] Ghosh S, Lee K, Moorthy S. Multiple scale analysis of heterogeneous elastic structures using homogenization theory and Voronoi cell FE method. Int J Solids Struct 1995;32(1):27–62.

[14] Moorthy S, Ghosh S. A Voronoi cell FE model for particle cracking in composite materials. Comput Methods Appl Mech Eng 1988;151:377–400. [15] Moorthy S, Ghosh S. Adaptivity and convergence in the Voronoi cell FE model for analyzing heterogeneous materials. Comput Methods Appl Mech Eng 2000;185:37–74. [16] Ghosh S, Lee K, Raghavan P. A multi-level computational model for multiscale damage analysis in composite and porous materials. Int J Solids Struct 2000;38 (14):2335–85. [17] Ghosh S, Paquet D. Adaptive multi-level model for multi-scale analysis of ductile fracture in heterogeneous aluminum alloys. Mech Mater 2013;65: 12–34. [18] Fish J, Yu Q, Shek KJ. Computational damage mechanics for composite materials based on mathematical homogenization. Int J Numer Meth Eng 1999;45:1657–79. [19] Fish J, Shek KJ. Multiscale analysis of composite materials and structures. Compos Sci Technol 2000;60:2547–56. [20] Könke C. Coupling of micro- and macro-damage models for the simulation of damage evolution in ductile materials. Comput Struct 1997;64:643–53. [21] Lemaitre J, Desmorat R. Engineering damage mechanics. New York: Springer; 2005. [22] Scotta R, Vitaliani R, Saetta A, Oñate E, Hanganu A. A scalar damage model with a shear retention factor for the analysis of reinforced concrete structures: theory and validation. Comput Struct 2001;79:737–55. [23] Arruda MRT, Manuel L, Castro S. A new hybrid-mixed stress model for the analysis of concrete structures using damage mechanics. Comput Struct 2013;125:23–44. [24] Paredes JA, Barbat AH, Oller S. A compression–tension concrete damage model, applied to a wind turbine reinforced concrete tower. Eng Struct 2011;33 (12):3559–69. [25] Erduran E, Yakut A. Component damage functions for reinforced concrete frame structures. Eng Struct 2007;29:2242–53. [26] Kim TH, Lee KM, Chung YS, Shin HM. Seismic damage assessment of reinforced concrete bridge columns. Eng Struct 2005;27:576–92. [27] Kamaris GS, Hatzigeorgiou GD, Beskos DE. A new damage index for plane steel frames exhibiting strength and stiffness degradation under seismic motion. Eng Struct 2013;46:727–36. [28] Fish J, Yu A. Multiscale damage modelling for composite materials: theory and computational framework. Int J Numer Meth Eng 2001;52:161–91. [29] Markovic D, Niekamp R, Ibrahimbegovic A, Matthies HG, Taylor RL. Multi-scale modeling of heterogeneous structures with inelastic constitutive behaviour: Part I – physical and mathematical aspects. Eng Comput 2005;22:664–83. [30] Mazars J, Pyaudier-Cabot G. Continuous damage theory-application to concrete. J Eng Mech 1989;115(2):345–65. [31] Lu XZ, Ye LP, Pan P, Tang DY, Qian JR. Pseudo-static collapse experiments and numerical prediction competition of RC Frame structure II: key elements experiment. Build Struct 2012;42(11):23–6 (in Chinese). [32] Lemaitre J. Handbook of materials behavior models: failures of materials. San Diego: Academic Press; 2001. [33] Sun B, Li ZX. Adaptive image-based method for integrated multi-scale modeling of damage evolution in heterogeneous concrete. Comput Struct 2015;152:66–81. [34] Labergère C, Rassineux AA, Saanouni KK. 2D adaptive mesh methodology for the simulation of metal forming processes with damage. Int J Mater Form 2011;4:317–28. [35] Toi Y, Hasegawa K. Element-size independent, elasto-plastic damage analysis of framed structures using the adaptively shifted integration technique. Comput Struct 2011;89:2162–8.

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Adaptive mesh refinement FEM for seismic damage evolution in concrete-based structures Bin Sun, Zhaoxia Li ⇑ Department of Engineering Mechanics, Jiangsu Key Laboratory of Engineering Mechanics, Southeast University, Nanjing 210096, China

a r t i c l e

i n f o

Article history: Received 30 June 2014 Revised 11 February 2016 Accepted 16 February 2016

Keywords: Damage Concrete-based structures Multi-grid FEM Reinforced concrete column

a b s t r a c t In order to study the seismic damage mechanisms of concrete-based structures, a adaptive multi-grid method is developed to simulate the process from evolving damage to structural failure in concrete structures under seismic loading. The theory of multi-grid FEM coupled evolving damage is developed on the basis of the improved variational principle, in which the elements in each sub-domain with different grid sizes are under the different state of damage. Meanwhile the FE elements can be refined adaptively with the accumulation of its damage for lower cost, without user intervention in the computation. As a case study of the method, the failure process of reinforced concrete column under seismic loading is simulated, the results show that, the simulated results fit well with the experiment, and the developed method can be used to reveal the seismic damage mechanism of concrete-based structures by considering the process from material damage to structural failure; the adaptive multi-grid method can give satisfactory result with lower cost by compared with experiment. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Due to its good qualities and low cost, concrete is widely used in civil engineering structures [1–3], e.g. concrete buildings and bridges, which exist extensively in the world. In case of the concrete structures failure during the service period, such as suffering earthquake, it will inevitably lead to disaster results. So in order to avoid catastrophe, study on seismic damage of the concrete-based structures to assess their seismic performance has become an important issue nowadays [4–6]. One of the characteristics of the failure of structure, especially for the structure of brittle material, is the sample-specificity [7–9], namely there are significant disparities at failure from sample to sample, even though they are initially identical at the macro-scale. It is because that the effect of material damage on global behavior of structure can be amplified strongly due to sensitive trans-scale coupling [9]. Therefore, in order to study well the failure mechanism of structure, it is necessary to simulate the coupling process of material damage up to structure failure. In the failure process of structures, the initiation and growth of material damage is an important phase. Damage is a physical process which leads to the progressive deterioration of structures due to initiation and growth of micro-cracks or micro-voids leading to ⇑ Corresponding author. E-mail address: [email protected] (Z. Li). http://dx.doi.org/10.1016/j.engstruct.2016.02.021 0141-0296/Ó 2016 Elsevier Ltd. All rights reserved.

the formation of local failure and structural failure for structures [10,11]. Therefore, in order to study on the seismic failure mechanism of concrete structures, the material damage of concrete should first be described. Since the beginning of the last century, the damage evolution processes in solids were attempted to be described [12]. At the beginning of the study, only few simple problems of damage can be solved without effective numerical methods. With the development of numerical methods especially for FE method (FEM), more and more damage evolution processes can be simulated numerically. Ghosh et al. have developed a voronoi cell FE method (VCFEM) to study the damage evolution in composites as well as heterogeneous aluminum alloys [13–17]. Fish delineated the damage phenomena in composite materials using mathematical homogenization method [18,19]. Könke has proposed a simulation method of damage evolution for ductile materials, which divides complete damage range into the micro- and macro-damage range [12,20]. Feng developed an adaptive mesh refinement FEM algorithm to demonstrate the damage evolution of brittle materials [9]. The methods for describing damage evolution can be mainly classified into two classes. The first case implements meso-structure of material in computing model to describe the damage evolution [13–19], while the second case uses a damage variable to describe the deterioration process of materials based on continuum damage mechanics (CDM) ignoring the meso-structure of material [9,12,20]. For the first case, the focus of attention is the damage evolution of materials and the scales

156

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

of the research objects are significantly smaller than the structural dimensions under the limitation of computing power. For the large scale engineering structures, CDM is an effective numerical tool [21]. Although this may decrease the accuracy due to ignoring the realistic meso-structure and micro-damage (e.g., micro-cracks and micro-voids), it brings the convenience of calculation. For concrete structures, CDM is also an effective numerical tool with wide applications to describe the damage evolution of concrete [22–24]. Most previous numerical methods use a damage index to evaluate the damage state and predict the structural failure of concrete structures based on the macroscopic response of structure ignoring the coupling process from material damage to local damage and failure and eventually to global structural failure [25–27]. Although some theories and computational frameworks were attempted to be developed [28,29], the validity of these frameworks had not been proved by experiment. In order to obtain an accurate estimation of damage state of concrete structures under seismic loading, it is necessary to develop a new method to simulate the process from evolving damage to structural failure, which is the basic aim of the paper. As a case of study of the proposed method, simulation of the process from evolving damage to failure of a concrete column under cyclic loading is given and compared with experiment. 2. FEM formulation of the coupling multi-grid interface and evolving damage

tensor and strain tensor in Xk , Dk is damage variable in Xk , Sr is the i . boundary surface subjected to the action of external surface force p

P12 ¼ P23 ¼ P31 ¼

e ð1AÞ 1 p e exp BðAee Þ ½ p

where ters,

0

e 6 ep e > ep

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ P3 jei jþei ; i¼1 hei i hei i ¼ 2

@ X2 \@ X3

ð1Þ

X1

@ X3 \@ X1

@ X1 \@ X2

ð6Þ

#

de 1 ij

@ X1 \@ X2

Pk þ P12 þ P23 þ P31

@ X3 \@ X1

where Pk is the functional of sub-domain Xk (k ¼ 1; 2; 3), P12 is the functional of interface between the boundary of X1 and the boundary of X2 ð@ X1 \ @ X2 Þ, P23 is the functional of interface between the boundary of X2 and the boundary of X3 ð@ X2 \ @ X3 Þ, P31 is the functional of interface between the boundary of X3 and the boundary of X1 ð@ X1 \ @ X3 Þ. The coupling of the sub-domains with different grid sizes can be obtained using the Lagrange multipliers. And the different terms can be rewritten as:

Z

i du1i dS p

1 k31 i dui dS ¼ 0

1 2 dk12 i ðui ui ÞdS ¼ 0

@e

2 ij

Z

@ X1 \@ X2

ð8Þ

#

Z

de2ij F i du2i dX2

Sr \X2

Z 2 k12 i dui dS þ

@ X2 \@ X3

i du2i dS p

2 k23 i dui dS ¼ 0

@ X2 \@ X3

2 3 dk23 i ðui ui ÞdS ¼ 0

Equilibrium in X3 :

" @Aðe3ij ; D3 Þ

@e

3 ij

Z

@ X2 \@ X3

#

Z

de3ij F i du3i dX3 Z

3 k23 i dui dS þ

@ X3 \@ X1

Sr \X3

i du3i dS p

3 k31 i dui dS ¼ 0

@ X3 \@ X1

3 1 dk31 i ðui ui ÞdS ¼ 0

Aðeij ; DÞ ¼

1 1 rij eij ¼ aijkl eij ekl ð1 DÞ 2 2

ð13Þ

where rij ¼ aijkl ekl ð1 DÞ is the Cauchy stress tensor [21], aijkl is the initial elastic stiffness tensor. According to the concept of effective ~ ij can be expressed as [21]: stress, the effective stress tensor r

r~ ij ¼

rij

1D

ð14Þ

In Eqs. (7), (9), (11), we use the relations:

where damage is considered to be related to strain.

ekij are respectively the displacement

ð12Þ

The body of concrete can be approximately considered as elastic body [9]. And the strain energy density coupling damage of elastic body can be written as [21]:

distributed body force, uki and

Sr \Xk

ð3Þ

ð11Þ

Equilibrium on @ X3 \ @ X1 :

where Aðekij ; Dk Þ is the strain energy density, F i ði ¼ 1; 2; 3Þ is

Xk

i uki dS ðk ¼ 1; 2; 3Þ p

ð9Þ

ð10Þ

1 2 @Aðeij ; DÞ 1 @D ~ ij eij ¼ rij r @ eij 2 @ eij

Pk ¼

ð7Þ

Equilibrium on @ X2 \ @ X3 :

Z

Z

ð2Þ

Sr \X1

Z 1 k12 i dui dS

" @Aðe2ij ; D2 Þ

X2

Z

Z

dX1

F i du1i

Equilibrium on @ X1 \ @ X2 :

Z

k¼1

i Aðekij ; Dk Þ F i uki dXk

3 1 k31 i ðui ui ÞdS

@ e1ij

Z

þ

X3

3 X

h

ð5Þ

Z

" @Aðe1ij ; D1 Þ

Z

ei is the principal strain :

In order to obtain satisfactory result with lower cost, computational domain with different grid sizes can be adaptively generated and changed with the damage evolution in the developed method. As shown in Fig. 1, there are three sub-domains with different grid sizes, which respectively are X1 ; X2 and X3 : And the sum of functionals of all sub-domains and interfaces P can be written as:

Z

2 3 k23 i ðui ui ÞdS

u1i ; u2i ; u3i are respectively the displacement tensor on the boundary of X1 ; X2 and X3 : Using the variational principle dP ¼ 0; the weak formulation of the equations can be expressed as: Equilibrium in X1 :

2.2. Variational formulation

P¼

ð4Þ

Z

X2 ; k23 i is Lagrange multipliers tensor which connect the X2 and X3 ; k31 is Lagrange multipliers tensor which connect the X3 and X1 ; i

Z

ep is the strain threshold of damage, A; B are model parame-

e¼

1 2 k12 i ðui ui ÞdS

where k12 i is Lagrange multipliers tensor which connect the X1 and

Based on CDM, the famous Mazars concrete damage model is chosen for describing the damage of concrete in this paper. And the Mazars concrete damage model can be written as follows [30]:

D¼

@ X1 \@ X2

Equilibrium in X2 :

2.1. Damage model for concrete

(

Z

eij ¼ ðui;j þ uj;i Þ

ð15Þ ð16Þ

157

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

Ω2

Ω1

Using the developed method to simulate the damage and failure process of structures

Ω3

Ω=Ω1 +Ω 2 +Ω3

Ω=Ω1 (a) Initial FE model consists only of coarse grid before the operation of the method

(b) FE model consists of multi-grid after the operation of the method Fig. 1. Adaptive split criterion.

2.3. Multi-grid FEM equations coupling multi-grid interface and damage The displacement uki ðk ¼ 1; 2; 3Þ in Xk ; and the Lagrange multipliers 23 31 k12 i ; ki ; ki

are interpolated from nodal values using shape functions as:

k

uk ¼ N uk

ðk ¼ 1; 2; 3Þ

ð17Þ

k12 ¼ N12 k K12

ð18Þ

N23 k K23 N31 k K31

ð19Þ

23

k

¼

k31 ¼

ð20Þ

Eqs. (7)–(12) can be rewritten as:

where L is the operator matrix transforming displacements into strains, which the has the similar form to that in the undamaged strain–displacement for small deformation problems, I is identity matrix, E is initial elastic stiffness matrix for undamaged material, is the prescribed surface F is the prescribed body force matrix, p force matrix, the matrix Mk ðk ¼ 1; 2; 3Þ is defined as:

Mkij ¼

@Dk k e @ ekij ij

The Newton–Raphson iterative method is chosen to solve the non-linear Eqs. (21)–(26), which requires the linearization of each equation, according to the linearized form:

Z

R X1 ¼

n o T T ðLN1 Þ ðI D1 ÞELN1 ðM1 LN1 Þ ELN1 u1 dX1 X1 Z Z dS N1T FdX1 N1T p X1 Sr \X1 Z Z N1T N12 N1T N31 þ k K12 dS k K31 dS ¼ 0 @ X1 \@ X2

Z RX1 \X2 ¼ R X2

@ X1 \@ X2

@ X3 \@ X1

N12T N1 u1 N2 u2 dS ¼ 0 k

n o T T ¼ X2 ðLN2 Þ ðI D2 ÞELN2 ðM2 LN2 Þ ELN2 u2 dX2 R R dS X2 N2T FdX2 Sr \X2 N2T p R R 2T 12 @ X1 \@ X2 N Nk K12 dS þ @ X2 \@ X3 N2T N23 k K23 dS ¼ 0

ð27Þ

n Rnþ1 ðÞ ¼ R ðÞ þ

@RnðÞ @bi

Dbi

ð : X1 ; X1 \ X2 ; X2 ; X2 \ X3 ; X3 ; X3 \ X1 Þ ð28Þ

ð21Þ

where ðnÞ denotes the value at nth iteration, bi are the independent variables. Using the Newton–Raphson iterative method, Eqs. (21)– (26) can be expressed as:

2 ð22Þ

R

ð23Þ

K1

6 1T 6K 6 12 6 60 6 6 60 6 60 4 K1T 31

K112

0

0

0

K131

0

K2T 12 0

0

0

K212

K2

K223

0

0

0

K2T 23

0

K3T 23 0

0

0

K323 K3

0

0

0

K3T 31

K331 0

@ X2 \@ X3

2 3 N23T k ðN u2 N u3 ÞdS ¼ 0

ð24Þ

Z

R X3

n o T T ¼ ðLN3 Þ ðI D3 ÞELN3 ðM3 LN3 Þ ELN3 u3 dX3 X3 Z Z dS N3T FdX3 N3T p X3 Sr \X3 Z Z N3T N23 N3T N31 k K23 dS þ k K31 dS ¼ 0 @ X2 \@ X3

@ X3 \@ X1

Z RX3 \X1 ¼

@ X3 \@ X1

7 7 7 7 7 7 7 7 7 7 5

3 3ðnÞ 2 RX1 Du1 ðnÞ 6 DK 7 6 RX \X 7 6 12 7 6 1 27 7 7 6 6 6 Du 2 7 6 RX2 7 7 7 6 6 7 6 DK 7 ¼ 6 R X2 \X3 7 6 23 7 6 7 7 6 6 4 Du 3 5 4 RX3 5 2

DK31

RX3 \X1 ð29Þ

Z RX2 \X3 ¼

3ðnÞ

N31T k

N3 u3 N1 u1 dS ¼ 0

where

Kk ¼

Z n o T T ðLNk Þ ðI Dk ÞELNk ðMk LNk Þ ELNk dXk ðk ¼ 1;2;3Þ ð30Þ Xk

Z Kk12 ¼

@ X1 \@ X2

NkT N12 k dS ðk ¼ 1;2;3Þ

ð31Þ

NkT N23 k dS ðk ¼ 1;2;3Þ

ð32Þ

NkT N31 k dS ðk ¼ 1;2;3Þ

ð33Þ

Z ð25Þ

Kk23 ¼

ð26Þ

Kk31 ¼

@ X2 \@ X3

Z

@ X3 \@ X1

158

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

Repeat the iterative process until nth the iteration until the convergence criterion is met:

2 R 2 3 R 3ðnÞ 1T 1T RX dS X1 N FdX1 þ Sr \X1 N p 1 6 6 7 6 6 RX1 \X2 7 7 0 7 6 7 6 , 7 6 R N2T FdX þ R 6 R dS 7 7 N2T p X2 6 X2 6 2 7 Sr \X2 7 6 6 7 6 tol 6 RX2 \X3 7 7 6 0 7 6 R 6 7 7 R 6 4 R 7 3T 3T 5 dS 5 X3 4 X3 N FdX3 þ Sr \X3 N p RX3 \X1 0 2

2

ð34Þ where tol is tolerance of equilibrium.

adaptively refined into 8 elements consisting of X3 , once its damage value reaches D2c . Third, a global damage index of structures is developed to evaluate the global seismic damage state of the whole structures due to material damage evolution and local failure. And the damage index reflects the degradation of the structural behavior. When structure loses bearing capacity, it is considered to be failure. For example, with accumulation of the material damage and local failure area, the force of the load point is close to 0 even if the displacement loading acting on the load point is very large, while the displacement of the load point is very large even if the force loading acting on the load point is very small. So the global damage index Dglobal is established based on the load–displacement curve of structures:

Dglobal ¼ 1 3. Computational strategy for adaptive mesh refinement FEM coupling evolving damage 3.1. Description of the damage and failure process In order to study well damage and failure mechanism of concrete structures, it is necessary to simulate the process from damage evolution to structural failure. Depending on the purpose, there are three main processes to be simulated in the developed method. First, the process of material damage evolution should be delineated in the developed method. Here famous Mazars concrete damage model based on CDM is chosen and coupled in the stress–strain relationship to describe the initiation and growth of micro-cracks or micro-voids using a damage variable. Second, the process from the material damage to local failure of structure should be delineated in the developed method. When the damage of element reaches 1, the element is considered to be failure and deleted. Meanwhile, in order to obtain the relatively accurate initiation and growth position of the local failure in structure and enhance the accuracy of finite element solutions of problems with high gradients due to damage, FE meshes consisting of a lot of elements are necessary when using traditional methods for the whole small-scale simulation. This often leads to consume a lot of computation time and even make the analysis impossible. Therefore, in this paper a 3D adaptive mesh refinement procedure is developed and implemented in the method for sufficient precision and lower cost. By using the method coupled with the procedure, the elements can be refined to smaller scales adaptively, once their damage values reach preset value. As shown in Fig. 1, One element in X1 can be adaptively refined into 8 elements consisting of X2 , once its damage value reaches D1c . One element in X2 can be

ΔlN

Δl0 F0

(a) Initial undamaged stage

FN Dl N

F0 Dl 0

ð35Þ

where F N ; DlN are respectively the load and displacement of the load point of the structure at Nth loading, F 0 ; Dl0 are respectively the initial load and displacement of the load point, when the global structure is in elastic stage, as shown in Fig. 2. And a procedure for monitoring Dglobal is implemented in the method. Once the Dglobal reaches 1, the concrete column is considered to be failure. 3.2. The basic procedure of the developed method As described the strategy for the simulation of the process from material damage to structural failure in Section 3.1, the flow chart of the developed method is given in Fig. 3, which chooses route 1 after the position A. 3.3. Technique for overcoming ‘‘island effect” After the operation of the algorithm, which choose route 1 after the position A in Fig. 3, some initial coarse elements are refined to smaller scale elements and some elements are failure, whose damage reach the preset critical value Dc, as shown in Fig. 4a. Note that in Fig. 4a there appear some connect regions consisting of elements suspending in the air, which does not agree with the fact due to gravity. The cause of this condition is that the damage of all elements around the connect region reach 1. The connect region suspending in the air look like computational debris which are needed to be eliminated automatically in the developed method. We define the formation of in the connect region suspending in the air as ‘‘island effect”, and the element in the region as ‘‘island element”. To eliminate the island elements, a procedure is developed and implemented in the method to overcome island effect. And the flow chart of the procedure is in the dashed boxes

Damage

FN

(b) Damaged stage at Nth loading Fig. 2. Definition of global damage index.

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

159

Fig. 3. Flow chart of the method.

of Fig. 3. After the operation of the algorithm, which choose route 2 after the position A in Fig. 3, the island elements can be automatically deleted, as shown in Fig. 4b. 4. Numerical analyses on seismic damage evolution of a reinforced concrete column

the numerical instability. And the ultimate strength of steel reinforcement bar can be founded in the experimental reports of the reinforced concrete structure under seismic loading [31]. Based on this way, steel reinforcement bar can be considered in the simulation of damage evolution of the reinforced concrete column. 4.2. Experimental validation

4.1. Description of steel reinforcement bar in the model of the reinforced concrete column Since steel reinforcement bar is also very important in the reinforced concrete structures except concrete, steel reinforcement bar should be considered in the simulation. And in order to consider the effect of steel reinforcement bar in the reinforced concrete column, which provides tensile strength, discrete steel reinforcement elements are coupled in the FE model as shown in Fig. 5b. Once steel reinforcement bar element reaches to the ultimate strength, the element is considered to be failure for preventing

Before using the developed method, parameters for material property of concrete should first be determined. And here the model parameters can be obtained as following way. Based on the given standard intervals for the model parameters of concrete in Ref. [32], adjusting macro-model parameters of concrete makes that the macroscopic analysis is in good accordance with experimental results, such as stress–strain curves. Based on this way, which is similar to the Ref. [1], a good set of model parameters can be obtained and the details can be found in our previous work [33].

160

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

Island elements

(b) Multi-grid of concrete which choose route 2 in the method

(a) Multi-grid of concrete which choose route 1 in the method

Fig. 4. Calculation result made up of concrete elements after the operation of the algorithm.

(b) Initial FE model

(a) Test setup of concrete column Fig. 5. Numerical example of RC column.

mm 60

Horizontal displacement loading

KN

Horizontal force loading

30 20 10 0 -10 -20

40 20 0 -20 -40

-30 0

0.02

0.04

0.06

0.08

0.1

0.12

N / Nf

0.14

0.165

-60 0.165

0.3

0.4

0.5

0.6

0.7

0.8

0.9

N / Nf

(b) Horizontal displacement loading

(a) Horizontal force loading Fig. 6. Load spectrum.

1

161

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

(a1) Calculation result at the cycle ratio

(b1) Calculation result at the cycle ratio

N / N f =29/38

N / N f =33/38

(c1) Calculation result at the cycle ratio N / N f =1

Damage

(a2) Experimental result at the cycle ratio

(b2) Experimental result at the cycle ratio

N / N f =29/38

N / N f =33/38

(c2) Experimental result at the cycle ratio N / N f =1

Fig. 7. Comparison between calculation and experimental results.

(a1) Undamaged state at N / N f =5/38

(a2) Initiation of damage at N / N f =6/19

(a3) Accumulation of damage at N / N f =23/38

(a5) Accumulation of (a4) Initiation of damage and local local failure area failure area at at N / N f =13/19

(a6) Failure of concrete column at N / N f =1

N / N f =31/38

Damage

Fig. 8. Simulation of the trans-scale failure process.

162

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164 1 0.9

In order to study the seismic damage mechanisms of concrete structures, the columns under cyclic loading is chosen as the numerical example, since its collapse experiment had been carried out by Tsinghua University [31]. As shown in Fig. 5. The dimension of the concrete column is 200 mm 200 mm 850 mm. The load spectrum is shown in Fig. 6. And N is the current number of loading cycles, the failure number of loading cycles N f is 38. As shown in Fig. 7, the failure evolution pattern predicted by simulation using the developed method agrees well with the experimental result by comparison with experiment. And the location of the initiation and growth of local failure area predicted by simulation is verified through experiment.

Calculation result Calculation result Experimental result Experimental result

0.8 0.7

Dglobal

0.6 0.5 0.4 0.3 0.2 0.1 0

4.3. Numerical simulation of the damage and failure process 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

N / Nf Fig. 9. Global damage index of the concrete column under the cyclic loading.

(a1) With scale-1 at N / N f =29/38

(b1) With scale-2 at N / N f =29/38

(c1) With scale-3 at N / N f =29/38

(a2) With scale-1 at N / N f =33/38

(b2) With scale-2 at N / N f =33/38

(c2) With scale-3 at N / N f =33/38

Compared with the previous numerical methods, one of the main features of the developed method is that the process from evolving material damage to structural failure can be simulated concurrently in a coupled manner. And the damage and failure pro-

(a3) With scale-1 at N / N f =1

(b3) With scale-2 at N / N f =1

(c3) With scale-3 at N / N f =1

Damage Fig. 10. Failure evolution patterns predicted by simulation with different grid scales.

163

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

235.1

25

106496

12.5

33.3 1.9

11.8

(a) Computational time

1664

13312

6.25

13298

(b) Total number of element

6.25

(c) Grid size of serious damaged area

Fig. 11. Comparison between the single scale simulation and adaptive multi-grid simulation.

cess of the concrete column under the cyclic loading is simulated using the developed method, as shown in Fig. 8. The process from undamaged state to initiation of material damage, to initiation and growth of local failure area due to accumulation of material damage, and eventually to structural failure, can be simulated using the developed method. Meanwhile the global damage index of the concrete column under the cyclic loading can be obtained from the monitoring data in the numerical method and experimental data based on Eq. (35), as shown in Fig. 9. For such a complex process, it is unavoidable that there exist error between the simulation and experiment. There are many causes. For example, the highly complex damage mechanisms of concrete structures cannot be completely contained in the numerical method. And the used famous concrete damage model cannot completely reflect the behavior of material damage. Similarly, it may also be due to the measuring error under the limitation of the experimental condition, or the real constraints for experimental condition cannot be completely reflected in simulation. When the global damage index reaches 1, the concrete column is considered to be failure. 5. Computational efficiency of the adaptive multi-grid method From recently researches, it shows that the calculated results of damage analysis using FEM considerably depend upon FE mesh [34,35]. In order to enhance the accuracy of finite element solutions of problems with high gradients due to damage, the refinement grid is required in the damaged area. But if using traditional methods for the whole small-scale simulation, this often leads to consume a lot of computation time and even make the analysis impossible. So a new numerical method coupled with 3D adaptive mesh refinement procedure is developed here to overcome the difficulty. 5.1. Numerical analysis by using the method in single scale In order to validate the efficiency of adaptive multi-grid method, we compare the computing cost of single scale simulations with different grid sizes and the adaptive multi-grid simulation, as shown in Fig. 10. The scales of grid with different colors in the singe scale simulation is given in Fig. 10, are the same as the scales of elements with the same color in the adaptive multi-grid simulation, as shown in Fig. 4. The element size of scale-1 is twice the size of scale-2, and the element size of scale-2 is twice the size of scale-3. And the simulation of damage and failure process using the method in single scale is given Fig. 10. As shown in Fig. 10, FE mesh size plays a role in simulating the precise initiation location of the local failure area and the pattern

for local failure area growth. The finer FE mesh, the more disordered details of local failure area initiation and growth can be simulated, including the specific position of initiation and the path of the growth, but it cost too much computational cost. 5.2. Comparison of efficiency As shown in Figs. 7 and 10, the adaptive multi-grid method can give satisfactory prediction with sufficient precision and lower cost by comparison with experiment and traditional pure small-scale numerical method, e.g., as shown in Fig. 11, the total numbers of the concrete elements of the two simulations are 13,298 and 106,496, respectively, and the computational time of the two simulations are 33.3 h and 235.1 h, respectively. Since the simulation of the process from material damage to failure of the reinforced concrete column is considerably nonlinear here, the computational time seems long here. But compared with the traditional single whole small-scale simulation method, the developed adaptive method indeed can save computational time and memory, and the developed method’s efficiency is verified. Therefore, the adaptive multi-grid FEM provides a numerical tool to simulate the process from material damage evolution to structural failure of the concrete-based structures with sufficient precision and lower cost. And the adaptive multi-grid method is much more economical in memory and time consumed for computation. 6. Conclusions The following specific conclusions can be obtained from the present study: (1) The theory and strategy for the adaptive mesh refinement FEM is proposed to simulate the process from material damage evolution to structural failure. (2) In the developed multi-grid FEM, FE elements can be refined adaptively with the damage accumulation without user intervention in computation. (3) Compared with the traditional method in single scale, the developed adaptive multi-grid method provides an excellent tool to obtain satisfactory result with sufficient precision and lower cost. (4) The developed method has been successfully applied to simulate the process from evolving material damage to global failure of a reinforced concrete column under cyclic loading. And the simulated results fit well with the experimental data. It is shown that the developed method can be used to reveal the damage and failure mechanisms of concrete-

164

B. Sun, Z. Li / Engineering Structures 115 (2016) 155–164

based structures by considering the process from material damage in concrete of stress concentration zone to local failure in vulnerable area and eventually to structural failure.

Acknowledgements The work described in this paper was substantially supported by A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (No. CE02-2-46), Scientific Research Foundation of Graduate School of Southeast University (No. YBJJ1440), and National Natural Science Foundation of China (No. 11072060). The authors are very grateful to the reviewer for carefully reading the paper and for his/her comments and suggestions which have improved the paper quality. References [1] Wriggers P, Moftah SO. Mesoscale models for concrete: homogenisation and damage behaviour. Finite Elem Anal Des 2006;42:623–36. [2] Zhou XQ, Hao H. Mesoscale modelling of concrete tensile failure mechanism at high strain rates. Comput Struct 2008;86:2013–26. [3] Ghosh A, Chaudhuri P. Computational modeling of fracture in concrete using a meshfree meso–macro-multiscale method. Comput Mater Sci 2013;69: 204–15. [4] Kim TH, Lee KM, Chung YS, Shin HM. Seismic damage assessment of reinforced concrete bridge columns. Eng Struct 2005;27(4):576–92. [5] Del Gaudio C, Ricci P, Verderame GM, Manfredi G. Development and urbanscale application of a simplified method for seismic fragility assessment of RC buildings. Eng Struct 2015;91:40–57. [6] Bai JL, Ou JP. Earthquake-resistant design of buckling-restrained braced RC moment frames using performance-based plastic design method. Eng Struct 2016;107:66–79. [7] Xia MF, Song ZQ, Xu JB, Zhao KH, Bai YL. Sample-specific behavior in failure models of disordered media. Commun Theor Phys 1996;25:49–54. [8] Xia MF, Wei YJ, Ke FJ, Bai YL. Critical sensitivity and trans-scale fluctuations in catastrophic rupture. Pure Appl Geophys 2002;159:2491–509. [9] Feng R, Xia MF, Ke FJ, Bai YL. Adaptive mesh refinement FEM for damage evolution of heterogeneous brittle media. Modell Simul Mater Sci Eng 2005;13:771–82. [10] Lin J, Liu Y, Deng TA. A review on damage mechanisms, models and calibration methods under various deformation conditions. Int J Damage Mech 2005;14:299–319. [11] Arun CO, Rao BN, Srinivasan SM. Stochastic mesh free method for elastoplastic damage analysis. Comput Methods Appl Mech Eng 2010;199: 2590–606. [12] Könke C. Damage evolution in ductile materials: from micro- to macrodamage. Comput Mech 1995;15:497–510. [13] Ghosh S, Lee K, Moorthy S. Multiple scale analysis of heterogeneous elastic structures using homogenization theory and Voronoi cell FE method. Int J Solids Struct 1995;32(1):27–62.

[14] Moorthy S, Ghosh S. A Voronoi cell FE model for particle cracking in composite materials. Comput Methods Appl Mech Eng 1988;151:377–400. [15] Moorthy S, Ghosh S. Adaptivity and convergence in the Voronoi cell FE model for analyzing heterogeneous materials. Comput Methods Appl Mech Eng 2000;185:37–74. [16] Ghosh S, Lee K, Raghavan P. A multi-level computational model for multiscale damage analysis in composite and porous materials. Int J Solids Struct 2000;38 (14):2335–85. [17] Ghosh S, Paquet D. Adaptive multi-level model for multi-scale analysis of ductile fracture in heterogeneous aluminum alloys. Mech Mater 2013;65: 12–34. [18] Fish J, Yu Q, Shek KJ. Computational damage mechanics for composite materials based on mathematical homogenization. Int J Numer Meth Eng 1999;45:1657–79. [19] Fish J, Shek KJ. Multiscale analysis of composite materials and structures. Compos Sci Technol 2000;60:2547–56. [20] Könke C. Coupling of micro- and macro-damage models for the simulation of damage evolution in ductile materials. Comput Struct 1997;64:643–53. [21] Lemaitre J, Desmorat R. Engineering damage mechanics. New York: Springer; 2005. [22] Scotta R, Vitaliani R, Saetta A, Oñate E, Hanganu A. A scalar damage model with a shear retention factor for the analysis of reinforced concrete structures: theory and validation. Comput Struct 2001;79:737–55. [23] Arruda MRT, Manuel L, Castro S. A new hybrid-mixed stress model for the analysis of concrete structures using damage mechanics. Comput Struct 2013;125:23–44. [24] Paredes JA, Barbat AH, Oller S. A compression–tension concrete damage model, applied to a wind turbine reinforced concrete tower. Eng Struct 2011;33 (12):3559–69. [25] Erduran E, Yakut A. Component damage functions for reinforced concrete frame structures. Eng Struct 2007;29:2242–53. [26] Kim TH, Lee KM, Chung YS, Shin HM. Seismic damage assessment of reinforced concrete bridge columns. Eng Struct 2005;27:576–92. [27] Kamaris GS, Hatzigeorgiou GD, Beskos DE. A new damage index for plane steel frames exhibiting strength and stiffness degradation under seismic motion. Eng Struct 2013;46:727–36. [28] Fish J, Yu A. Multiscale damage modelling for composite materials: theory and computational framework. Int J Numer Meth Eng 2001;52:161–91. [29] Markovic D, Niekamp R, Ibrahimbegovic A, Matthies HG, Taylor RL. Multi-scale modeling of heterogeneous structures with inelastic constitutive behaviour: Part I – physical and mathematical aspects. Eng Comput 2005;22:664–83. [30] Mazars J, Pyaudier-Cabot G. Continuous damage theory-application to concrete. J Eng Mech 1989;115(2):345–65. [31] Lu XZ, Ye LP, Pan P, Tang DY, Qian JR. Pseudo-static collapse experiments and numerical prediction competition of RC Frame structure II: key elements experiment. Build Struct 2012;42(11):23–6 (in Chinese). [32] Lemaitre J. Handbook of materials behavior models: failures of materials. San Diego: Academic Press; 2001. [33] Sun B, Li ZX. Adaptive image-based method for integrated multi-scale modeling of damage evolution in heterogeneous concrete. Comput Struct 2015;152:66–81. [34] Labergère C, Rassineux AA, Saanouni KK. 2D adaptive mesh methodology for the simulation of metal forming processes with damage. Int J Mater Form 2011;4:317–28. [35] Toi Y, Hasegawa K. Element-size independent, elasto-plastic damage analysis of framed structures using the adaptively shifted integration technique. Comput Struct 2011;89:2162–8.