Static and Dynamic Analysis of Reinforced Concrete Shells

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.


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 9node 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. 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 rota- 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:   [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 [ ] 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: Latin American Journal of Solids and Structures 10(2013) 1109 -1134 with, 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 = 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 σ 0 is the effective stress. In addition, the associated flow rule is defined as: with the flow vector given by: 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: [ ] 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. 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: Latin American Journal of Solids and Structures 10(2013) 1109 -1134 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 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 [ ] D degradation, will be described in terms of first and second invariants only. Both surfaces are expressed in the following manner: where I 1 and J 2 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 w p f 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 α 1 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: Figure 4 Basic concept of the model taken from Gomes [5] The parameters β o and β 1 are determined from experimental test and β o f 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 γ (  ε e ) is related to the elastic strain rate through an exponential function of an effective elastic strain rate.
where a o and a 1 are parameters which must be determined experimentally. The effective elastic strain is defined as: because deviatoric strains cause most damage to concrete and ε e eff is equal to the uniaxial elastic strain for a uniaxial stress state, being J 2e the second invariant of the deviatoric elastic tensor.

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: 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.
If required, form the effective stiffness matrix using the expression Otherwise use a previously calculated K * 5 Factorize, forward reduction and back substitute as required to solve If Δd i and/or ψ i do not satisfy the convergence conditions then set i = i + 1 and go to step 3, Otherwise continue 8 Set For use in the next time step. Also set n = n+1, form Cv n+1 + B T σ n+1 (d n+1 )dΩ ∫ and begin Next time step 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.

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  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].    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 2 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 2 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 2 and 3000 Kg/cm 2 ) have a considerable effect in the response at the impact zone. In this work, an elastic modulus of 27000 Kg/cm 2 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 2 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].

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 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.