Static and dynamic analysis of reinforced concrete shells

Static and dynamic analysis of reinforced concrete shells

Jorge L. P. Tamayo^I,*; Inácio B. Morsch^II; Armando M. Awruch^II

^ICEMACOM, Computational and Applied Mechanical Center, Engineering School, Federal University of Rio Grande do Sul, Av. Osvaldo Aranha, 99-3º Floor, 90035-190, Porto Alegre, RS, Brazil

^IIPPGEC, Department of Civil Engineering, Engineering School, Federal University of Rio Grande do Sul, Av. Osvaldo Aranha, 99-3º Floor, 90035-190, Porto Alegre, RS, Brazil

ABSTRACT

The objective of this work is to provide a reliable numerical model using the finite element method (FEM) for the static and dynamic analysis of reinforced concrete (RC) shells. For this purpose two independently computer programs based on plasticity and viscoplasticity theories are developed. The well known degenerated shell element is used for the static analysis up to failure load, while 3D brick elements are used for the dynamic application. The implicit Newmark scheme with predictor and corrector phases is used for time integration of the nonlinear system of equations. Two benchmark examples analyzed by others are solved with the present numerical model and results are compared with those obtained by other authors. The present numerical model is able to reproduce the path failure, collapse loads and failure mechanism within an acceptable level of accuracy.

Keywords: Reinforced Concrete (RC) Shells, Finite Element Method (FEM)

1 INTRODUCTION

Nonlinear analysis of reinforced concrete shells is an important subject nowadays. Thorough investigation of capacity and safety aspects, concrete structures require to establish that the entire structural response up to collapse. For this purpose, mathematical models for predicting the behavior of concrete in static and dynamic loading are presented. Based on plasticity and viscoplasticity theories, such predictions are possible. The finite element method is the natural choice for spatial discretization of complex structural systems. The basic equilibrium and incremental equations for this method are easily obtained using the principle of virtual work. In general, it may be stated that practical problems require insight and understanding in order to develop efficient methods of analysis. This paper discusses the development of two different computer programs for the analysis of reinforced concrete structures under static and dynamic loading. The first program is suitable for the static analysis of reinforced concrete shells using the theory of plasticity with the degenerated 9-node shell finite element proposed by Figueiras [4]. The second program handles the dynamic loading case using the theory of viscoplasticity using the 20-node brick finite element as presented by Gomes [5]. In Fig. 1, the two finite elements used in this work are shown. The brick element degenerates into the shell element when its upper and lower nodes joint together into its middle plane, being both essentially very similar. Also in the same figure, the natural coordinates systems are shown.

2 FINITE ELEMENT FORMULATION

2.1 Finite element formulation and elasto-plastic constitutive model of the degenerated shell element

The 9-node shell element is used here to represent the concrete shell structure where the reinforcing bars are modeled using the smeared layer approach. The displacement field within the element is defined in terms of quadratic shape functions and displacement values at the nodes. Each nodal point has three degrees of freedom u, v and w, displacement components along the cartersian coordinates x, y and z, respectively and two local degrees of freedoms α and β, which are rotations depending on the local coordinate system defined by the nodal vectors {V}₁, {V}₂ and {V}₃ as illustrated in the right side of Fig. 1. The degrees of freedoms u and v are associated to membrane actions, while w, α and β are associated to bending actions. Therefore, for each element the displacement vector is expressed in the following manner:

In addition, several concrete and steel layers are defined through thickness in order to capture properties variations due to nonlinearities. A plane stress assumption coupled with out of plane shear stresses is considered for each layer. Then, the strain components, in the local system of the element are given by:

where,

or

where {V}_1k and {V}_2k define the nodal coordinate system at each node k of the element (see Fig.1) and matrices [B]_k, [θ] and [C]_k are given explicitly in appendix A. In addition, [B]_s0 is the usual strain-displacement element matrix (e.g. see Figueiras [4]). For the degenerate shell elements employed in this work a specific total Lagragian formulation is adopted for modeling geometrical nonlinear behavior in which large deflections and moderate rotations are accounted for. Within this approach, the current 2nd Piola Kirchhoff stress tensor components and the Green-Lagrange strain tensor components are referred to the original geometric configuration and the displacement field gives the current configuration of the system with respect to its initial position. Refering variables to the original configuration is advantageous for degenerate shell elements, since the computationally expensive transfer of quantities between local and global axes need to be performed only once. The strain-displacement matrix [B]_s0 is calculated once during the nonlinear process ant its corresponding nonlinear part [B]_s1 to be described next is updated using the current displacement by a simple matrix product. Then, the nonlinear strain components vector is given by:

or

where [S]_k and [G]_k are matrices given explicitly in appendix A and contain the contribution of each nodal variable to the local derivates ∂w'/∂x' and ∂w'/∂y'. In addition, several concrete and steel layers are defined through thickness. Perfect adherence is considered between concrete and steel reinforcement, so Eq. (1) is also valid for steel layers. The stress components of the Piola-Kirchhoff stress tensor for concrete and steel materials in the element local coordinate system are determined by the following expression:

where [D] is the material constitutive matrix in the local coordinate system for a given integration point. The nodal equivalent forces {P}ⁱ and the stiffness element matrix [K]ⁱ_s, at a given iteration i, are obtained applying the principle of virtual works and are given by the following expressions:

with,

where [D]ⁱ_et is the uncracked, cracked or elasto-plastic tangential constitutive matrix for the concrete material and the elastic or elasto-plastic constitutive matrix for the steel reinforcement, being V the element volume. The second term in the right hand side of Eq. (8) represents the initial stress or geometric stiffness matrix where stress components matrix [σ] is defined in appendix A. A minimum number of eight concrete layers are found to be suitable for the integration of the above expressions besides a gaussian rule of 2x2 for the middle plane of each layer. Concrete in compression is modeled using the associated theory of plasticity; a modified Drucker-Prager yield criterion, which was proposed by Figueiras [4], is used in this work. Due to nonlinear hardening behavior, this yield criterion defines an initial yield surface at an effective stress equal to σ₀ = 0.3 f_c (which is the beginning of the plastic deformation) and a limit surface separating a nonlinear state from a perfect elasto-plastic one, as it is shown in Fig. 2. The yield criterion is defined as:

where σ₀ is the effective stress. In addition, the associated flow rule is defined as:

with the flow vector given by:

In Eq. (8), {dε} contains the components of the total strain, is a component of the plastic strain tensor, [D]_e is the elastic constitutive matrix and H' is the hardening parameter established as the slope of the uniaxial curve which defines the hardening rule. This curve known as "Madrid parabola" is defined by the following expression:

where E_c is the elastic modulus, ε_o represents the total strain at cracking compression stress f_c and _p is the effective plastic deformation. The elasto-plastic constitutive relation is expressed in the following differential form:

where [D]_ep is the elasto-plastic constitutive matrix. Finally, the crushing is given by:

where ε_u represents the ultimate deformation extrapolated from experimental test.

On the other hand, the response of concrete under tensile stresses is assumed to be linear elastic until the fracture surface is reached and then, its behavior is characterized by an orthotropic material. The cracking is governed by a maximum stress criterion. Cracks are assumed to occur in planes perpendicular to the direction of the maximum tensile stress as soon as this stress reaches the specified concrete tensile strength f_t. After cracking has occurred the elastic modulus and Poisson's ratio are assumed to be zero in the perpendicular direction to the cracked plane, and a reduced shear modulus is employed. Due to bond effects, cracked concrete carries, between cracks, a certain amount of tensile force normal to the cracked plane. This effect is considered through a relationship between the strain and the stress normal to the cracking plane direction, as shown in Fig. 3, where ε_ct is the strain associated to f_t and ε_tm is the maximum strain for ω between 0.5-0.7 and the normal stress σ_j is determined from a known value of strain ε_j. Further details of the constitutive model for concrete can be found in Figueiras [4]. The steel reinforcement is modeled as an uniaxial elasto-plastic material with a constant elastic steel E_s and a tangential modulus E_s' according to the bilinear stress-strain relation shown in the right side of Fig. 3.

2.2 Finite element formulation and elasto-viscoplastic constitutive model of the 3D brick element

The 20-node isoparametric quadratic brick element is used here to represent the concrete shell structure where the reinforcement bars are modeled also using the smeared layer approach. The displacement field within the element is defined in terms of the shape functions and displacement values at the nodes. Each nodal point has three degrees of freedom u, v and w along the cartersian coordinates x, y and z, respectively. Therefore, for each element the displacement vector is expressed in the following manner:

The strain components vector is defined by:

or

where N_k is the shape function of node k and [B]_b is the strain-displacement matrix. The stresses and strains are related by the following expression:

where [D] is the material constitutive matrix in the global system. Equivalent nodal forces, at a given iteration i, are expressed in the following manner:

The stiffness matrix for a concrete element of volume V can be expressed as:

where |J| is the determinant of the jacobian matrix and ξ, η and ζ represent the local natural coordinates at each integration point within the element. As it was mentioned earlier, the steel bars are assumed to be distributed over the concrete element in any direction in a layer with uniaxial properties. A composite concrete-reinforcement constitutive relation is used in this case. To derive such a relation, perfect bond is assumed between concrete and steel. Crushing criteria is similar to that presented in the previous section while the cracking monitoring algorithm presented in Gomes [5] is used.

Earlier developments and studies suggest that a concrete model intended for transient analysis should be rate and history dependent. Apparently, elasto-plastic and elasto-viscoplastic models in their basic formulations do not meet these requirements, but they can be modified to incorporate such effects. However, Cotosvos and Pavlovic [3] have given a detailed explanation for this apparent strength gain at high loading rates. They attribute this fact due to the effect of the inertia loads that reduce the rate of cracking and consequently lead to an increase of the load-carrying capacity of the structure. Therefore, the use of a constitutive model which is dependent on the strain loading rate is questionable. As the objective of this work is not to discuss this issue in detail, the wide acceptance rate and history dependent load model presented by Gomes [5] is still used.

The present model is a strain-rate sensitive elasto-viscoplastic model with progressive degradation of the strength and it introduces two main differences when compared with the traditional elasto-viscoplastic model. Firstly, the fluidity parameter is not constant and it is assumed to be dependent on the elastic strain rate and secondly, a variable strength limit surface is introduced to monitor the damage caused by the viscoplastic flow. If the stress point reaches the strength limit surface, then the degradation of the yield surface is initiated. The yield surface F_D, defining the onset of viscoplastic behavior, and the strength limit surface F_f, defining the initiation of material degradation, will be described in terms of first and second invariants only. Both surfaces are expressed in the following manner:

where I₁ and J₂ are the first and the deviatoric second stress invariants, respectively. The constants c and β are evaluated from experimental test and are equal to 0.1775 and 1.355, respectively. During inelastic straining both surfaces change, depending on the amount of accumulated damage, which is expressed as the viscoplastic energy density w_p and it is given by:

where t_f and are, respectively, the time and viscoplastic energy density when the strength limit is reached. While the stress path remains inside the yield surface, the behavior of concrete is linearly elastic, no viscoplastic straining is developed, and the yield and failure surfaces remain stationary. When the stress path is outside the yield surface, inelastic straining occurs, and the values of σ_o and σ_f vary. Expansion of the yield surface due to hardening is not considered here. However, σ_f decreases with the increase of damage and the strength limit surface shrinks. The strength limit surface is only a monitoring device to define where and when the failure occurs. When the stress path reaches the strength limit surface, degradation of the material is initiated. After that, the strength limit surface is not longer considered and the yield surface begins to shrink according to the post-failure dissipated energy density k. All those paths described before are well depicted in Fig. 4.

An exponential function will be used to describe the post-failure behavior. Therefore, the function σ_o(w_p, k) is defined by the following expression:

where α₁ defines the limit for elastic behavior (typically between 0.3-0.4) and α_c models the degradation after failure. The parameter f'_c is the static compressive strength of concrete. The failure stress will be assumed to be a linear function of the viscoplastic energy density and the function σ_f(w_p) is defined by the following expression:

The parameters β_o and β₁ are determined from experimental test and β_of'_c is the compressive strength obtained with infinite load rates and no inelastic strains. The rate of viscoplastic straining depends on the rate of elastic strain and on the position of the yield surface. It is written as:

The fluidity parameter is related to the elastic strain rate through an exponential function of an effective elastic strain rate.

where a_o and a₁ are parameters which must be determined experimentally. The effective elastic strain is defined as:

because deviatoric strains cause most damage to concrete and is equal to the uniaxial elastic strain for a uniaxial stress state, being J_2e the second invariant of the deviatoric elastic tensor.

3 NUMERICAL ALGORITHM

In order to introduce the implicit numerical algorithm for the solution of the dynamic equation, it is necessary to describe the predictor and corrector form of the Newmark scheme for the integration of the semi-discrete system of equations governing nonlinear transient dynamic problems. Typically at time station t_n+1 these equations take the following form:

where [M] and [C] are the mass and damping matrices while {a}_n+1, {v}_n+1 and {d}_n+1 are the acceleration, velocity and displacement vectors. The tangential stiffness matrix [K]_et is related to the internal forces in the following manner:

with

In the Newmark scheme the displacement and velocity at time t_n+1 can be expressed in the following form:

with

Note that {d}_n, {v}_n and {a}_n are the approximations to d(t_n), (t_n) and (t_n) and β and γ are free parameters which control the accuracy and stability of the method. {}_n+1 and {}_n+1 are the predictor values and {d}_n+1 and {v}_n+1 are the corrector values. Initially the displacements {d}₀ and velocities {v}₀ are provided and the acceleration {a}₀ is obtained from the following expression:

By using Eq. (30) to Eq. (36), an effective static problem is formed which is solved using a Newton Raphson type scheme. This algorithm is summarized in Table 1.

The above algorithm is also useful for static problems by simply omitting all the inertial terms such as the mass and damping matrices and also the velocity and acceleration vectors. In this case the time step is a fictitious time which is used as an incremental multiplier.

4 NUMERICAL APPLICATIONS

4.1 Geometrically nonlinear reinforced parabolic cylindrical shell experimentally tested by Hedgren and Billington (1967)

A parabolic cylindrical shell with variable thickness, subjected to uniformly distributed pressure load P, which was tested by Hedgren and Billington (see e.g. Figueiras [4]), is analyzed with the present finite element code, which includes a total Lagrangian scheme to analyze geometrically nonlinear problems. The shell structure was tested with end support diaphragms and free edges. A finite element mesh of 36 heterosis finite elements is used to model one quarter of the structure due to symmetry considerations. Material properties are listed in Table 2 and the layout plant and transversal cross section of the shell are shown in Fig. 5. For a detailed description of the layer pattern zones of the reinforcement, the reader is referred to the work of Figueiras [4]. Because different amount of reinforcement with different directions are present in various zones of the shell, this example is considered to be a very challenger case. Then, several steel layers need to be defined for each finite element to take into account space steel variation.

In Fig. 6, the load-vertical deflection curve at midspan of the free edge is compared with the results obtained by Figueiras [4] using a nonlinear geometrical model. As it can be observed good agreement is found at each load level. It is important to explain that all the results presented here are normalized with respect to the results obtained using a reference load P_r = 0.358 N/cm². The present numerical model predicts a failure load factor equal to 4.25 while the ultimate experimental load factor was 4.35. The vertical deflections of the transverse section at midspan of the shell are presented in Fig. 7 for load factors of 0.8, 2.4, 3.2 and 4.0. The horizontal displacements at support section are shown in Fig. 8 for load factors 0.4, 2.4 and 4.2. The transverse variation of stress resultants N_x and M_y at midspan are shown in Fig. 9 and Fig. 10 for load factors of 0.4 and 3.2, respectively.

Fig. 11 shows the stresses in the main reinforcement along the free edge for three different load levels. The full line corresponds to the bottom reinforcement and the dashed line refers to the reinforcement top layer. From the difference between the stresses in the top and bottom layer, the moment acting on the free edge can be estimated. It should be noted that the ultimate load has been found before the main reinforcement has reached its ultimate strength. All the previous results shown good agreement with the results reported by Figueiras [4].

4.2 Aircraft impact on nuclear power plant (1964)

The horizontal impact of an aircraft on the shield building of a nuclear power plant is analyzed (see Cervera and Hinton [2]). The geometry, the loading function and the reinforcement are specified in Fig. 12. The built-in reinforced concrete shell is composite of cylindrical and spherical parts of constant thickness. The reinforcement placed circumferentially and meridionally on the interior and exterior surfaces consist of bars of 40 mm diameter, spaced at 8 cm. The material properties are shown in Table 3. The impact is assumed to occur horizontally and is analyzed following an uncoupled procedure in which its effect onto the nuclear power plant is considered through the application of a short-time load (See References [2],[6],[8],[9] and [10]) as shown in Fig. 12. The location of the area of impact of 28 m² is also shown in Fig. 12. The load history is also indicated in the same figure and it is noted that the load has a maximum value of 9000 ton.

Since the loading and geometry of the shell are symmetric, only one half of the structure is modeled. A mesh of 54 solid elements is used in the analysis, with a local refinement in the vicinity of the impact load where a rectangular area of 14 m² is defined to apply the distributed load (see Fig. 13). The implicit Newmark scheme with β = 0.25 and γ = 0.5 is used to integrate in time with a time step Δt = 0.00475. In Table 4, a summary of results for the linear and nonlinear horizontal peak displacement at the impact zone found in different references is presented. It is observed that there exists a range of variation for this value even for the linear case. The variation in the nonlinear case is generally attributed to the use of different cracking strain values between 0.0015 and 0.0020, being some of them not reported in their original references. It is important to indicate that the peak horizontal displacement does not occur at the dome-cylinder junction (point A). It occurs some distance below point A, near the centre of the impact area. It is possible that some authors reported their results considering that point A is located at the centre of the impact area. Also, it was verified that small variations in the value of the elastic modulus of concrete (between 25000 Kg/cm² and 3000 Kg/cm²) have a considerable effect in the response at the impact zone. In this work, an elastic modulus of 27000 Kg/cm² was chosen because this value yielded similar results to those published by Cervera and Hinton [2] as shown later. An additional analysis was carried out using an elastic modulus equal to 30000 Kg/cm² and considering that point A is located at the centre of the impact area. Obtained results (not shown here) are in agreement with those obtained by Abbas et al. [1], whose response greatly differs from that obtained by Cervera and Hinton [2] for a cracking value 0.0015.

Horizontal displacements at points A, B and C are plotted as functions of time in Fig. 14, Fig. 15 and Fig 16, respectively. Two different values of cracking strain are considered here (0.0015 and 0.00185). In order to compared similar profiles of displacements (see Fig. 14), the response obtained using a cracking strain value 0.00185 is compared directly to that obtained by Cervera and Hinton [2] using a cracking value 0.0020 because it yields approximately the same peak displacement at time 0.274s. It is important to emphasize that some differences are expected due to the slightly different mesh used and because of the great sensibility of the response at the impact zone due to cracking. This means that a small variation in the value of the cracking strain could yield significant changes in displacements. Others facts to take into account in the final response are the sensibility established by the type of the integration rule used (an eighteen point selective integration rule was used) and the value of the covering of the reinforcement mesh, which is considered to be 10 cm in this work. Results obtained at point B are also expected to be slightly different because the finite element mesh presents a small circular opening at the upper part of the dome. In general, the profile patterns and magnitude of displacements at the three different locations are in good agreement with those presented by Cervera and Hinton [2].

5 CONCLUSIONS

In this work, the usual elasto-plastic and elasto-viscoplastic algorithms are used to predict the dynamic and static behavior of reinforced concrete shells. Validation of the present algorithms and codes are provided by modeling two usual benchmark examples found in the technical literature about this topic. These two examples are considered to be complex because of the nonlinearities involved and in the case of the static example also due to the complex reinforcement arrangement. For the dynamic analysis of the nuclear power plant, different results (even for the linear case) are reported in the literature by various authors. Then, some explanations are given to justify these differences such as that the final response is very sensitive to the way in which cracking develops within the impact zone, the location of the peak displacement and the value of the elastic modulus of concrete. Several tests have shown that mesh refinement near the impact zone, time step value for the dynamic analysis and covering value of the reinforcement mesh have little influence in the final results.

Otherwise, because the elasto-viscoplatic algorithm is explicit, a suitable time step size must be chosen for the stability of the integration of the constitutive relations. This time step size is related to the specific form of the viscoplastic potential used in the flow rule. In addition, the global stability of the procedure also depends on the time step length chosen for the integration of the dynamic equation of motion. These dependencies make the elasto-plastic algorithm more attractive because automatic sub-incrementation can be used instead of selecting a suitable time step size for the integration of the constitutive relations. In this way, the procedure is reduced only to the choice of a suitable time step length for the integration of the dynamic equation of motion. Finally, results obtained here are in good agreement with those numerically obtained by other authors. Further investigation is carried out to implement a strain rate sensitive elasto-plastic model for concrete and also to include other effects such as concrete creep and shrinkage.

Acknowledgments The financial support provided by CAPES and CNPQ is gratefully acknowledged.

Received 08 Sep 2012

In revised form 21 Nov 2012

Appendix A

For a given node k, the following expressions apply:

where [θ] is the transformation matrix, which relates the local coordinate system with the global one. The terms of the jacobian matrix [J] are given by the following expressions:

where the superscripts sup and inf stand for superior and inferior nodes. For the geometrical nonlinear analysis, the following matrices are needed:

with,

The term [G] is a matrix with two rows and a number of columns equal to the number of element nodal variables. The first row contains the contribution of each nodal variable to the local derivate ∂w'/∂x' (corresponding shape function derivates) and the second row contains similar contributions for ∂w'/∂y'.

Then, the terms of matrix [G] are rearranged to for the matrix [S] in the following manner:

The stresses components are also rearranged conveniently in the following manner:

• [1] H. Abbas, D.K. Paul, P.N. Godbole, and G.C. Nayak. Aircraft crash upon outer containment of nuclear power plant. Nuclear Engineering and Design, 160(1):13-50, 1996.
• [2] M. Cervera and E. Hinton. Análisis dinámico en rotura de estructuras laminares y tridimensionales de hormigón armado. Rev. Internacional de Métodos Numéricos para Calculo y Diseño en Ingeniería, 3(1):61-76, 1987.
• [3] D.M. Cotsovos and N. Pavolic. Numerical investigation of concrete subjected to compressive impact loading. Part1: A fundamental explanation for the apparent strength gain at high loading rates. Computers and Structures, 86(1-2):145-163, 2008.
• [4] J.A. Figueiras. Ultimate load analysis of anisotropic and reinforced concrete plates and shells, PhD. Thesis, University of Wales, U.K., 519, 1983.
• [5] H. Gomes, 1997, Análise da confiabilidade de estruturas de concreto armado usando o método dos elementos finitos e processos de simulação, Dissertação de Mestrado, Universidade Federal do Rio Grande do Sul, 1997.
• [6] M.A. Iqbal, S. Rai, M.R. Sadique, and P. Bhargava. Numerical simulation of aircraft crash on nuclear containment structure. Nuclear Engineering and Design, 243:321-335, 2012.
• [7] G.Q. Liu. Nonlinear and transient finite element analysis of general reinforced concrete plates and shells, PhD. Thesis, University of Wales, U.K., 1985.
• [8] M. Kukreja. Damage evaluation of 500 MWe Indian pressurized heavy water reactor nuclear containment for aircraft impact. Nuclear Engineering and Design, 235(17-19):1807-1817, 2005.
• [9] C.M. Madasamy, R.K. Singh, H.S. Kushwaha, S.C. Mahajan, and A. Kakodkar. Nonlinear transient analysis of Indian PHWR reinforced concrete containments under impact load A case study. Transactions of the 13^th International Conference on Structural mechanics in Reactor Technology (SMIRT), H038-1:243-248, 1995.
• [10] A.K. Pandey, R. Kumar, D.K. Paul, and D.N. Trikha. Strain rate model for dynamic analysis of reinforced concrete structures. Journal of Structural Engineering, 132(9):1393-1401, 2006.

Publication Dates

• Publication in this collection
  15 May 2013
• Date of issue
  Nov 2013

History

• Received
  08 Sept 2012
• Accepted
  21 Nov 2012