Mesolevel simulation of reinforced concrete structures under

Mesolevel simulation of reinforced concrete structures under impact loadings Gianluca Cusatis Rensselaer Polytechnic Institute, Troy (NY), USA Daniele Pelessone ES3, Solana Beach (CA), USA ABSTRACT: Reinforced concrete structures are often subjected to extreme dynamic loading conditions due to direct impact. This paper deals with the numerical simulation of the effect of impact loadings on reinforced concrete structural members such as columns, beams and slabs. Concrete is modeled by a three-dimensional mesolevel lattice model, which can predict accurately concrete behavior in a wide range of loading conditions. The adopted constitutive law simulates fracture, friction and cohesion at the meso-level and takes into account the strain rate dependence of concrete behavior. A mesh of plastic beam elements simulates the effect of steel reinforcement. Since the nodes of the reinforcement do not coincide, in general, with the nodes of the lattice system ad hoc coupling algorithms are proposed. The mesolevel model for concrete and the concreterebar coupling algorithm are implemented in the framework of the object oriented dynamic finite element code MARS. Finally, numerical simulations of the behavior of reinforced concrete structural members subjected to impact loading are performed and the numerical results are compared with experimental data gathered from the literature. 1 INTRODUCTION Nowadays, the assessment of the performance and the vulnerability of reinforced concrete structures subjected to transient dynamic loading conditions is of paramount importance. Typical examples of these loading conditions include transportation structures subjected to vehicle crash impact, marine and offshore structures exposed to ice impact, protective structures subjected to projectile or aircraft impact and structures sustaining shock and impact loads during explosions or earthquakes. This paper presents preliminary results obtained by the Authors in simulating the effect of impacts (in particular projectile penetration) on reinforced concrete slabs. Concrete is modeled by a recently developed three-dimensional mesolevel lattice model, which can predict accurately concrete behavior in a wide range of loading conditions, including uniaxial, biaxial, and triaxial compression, hydrostatic compression, uniaxial tension, biaxial tension (Cusatis et al. 2003b), mode I and mixed mode fracture (Cusatis et al. 2005). The formulation of the constitutive behavior at the mesolevel allows simulating fracture, friction and cohesion at the meso-level and takes into account the strain rate dependence of concrete behavior. This dependence, 1 which is crucial for simulating correctly the effect of impact loadings, is due to two different physical mechanisms: 1) the dependence of the fracture process on the rate of crack opening, and 2) the viscoelastic deformation of cement paste. The first mechanism is described by the activation energy theory applied to the ruptures that occur along a crack, and the second by the microprestress-solidification theory of concrete creep (Baˇ ant et al. 2004). The steel reinforcement is z simulated by a mesh of beam elements which interact with the lattice mesh through a coupling algorithm based on the master-slave approach. Numerical simulations of the behavior of reinforced concrete slab subjected to the impact of a steel projective are carried out by using the object oriented dynamic finite element code MARS (Pelessone 2005). The numerical results are compared with some experimental data gathered from the literature. ADOPTED CONSTITUTIVE MODEL FOR CONCRETE The model used in the present study, called Confinement-Shear Lattice model (CSL model), simulates the concrete meso-structure by taking into account only coarse aggregates. A lattice mesh con2 necting the aggregate centers reproduces the interaction between adjacent aggregate pieces through the embedding matrix. The constitutive law assumed for each lattice element models stiffness, strength and inelastic behavior, and includes the effect of the lower level structure that is not directly introduced in the formulation (cement paste, small aggregate particles, cement-aggregate interface). In the following we present a review of the model formulation, whose detailed description appears in Cusatis et al. (2003a), Cusatis et al. (2003b), and Cusatis et al. (2005). We use the following three-step algorithm to obtain the geometry of the meso-structure for a given concrete specimen. 1) A try-and-reject random procedure places each aggregate particle throughout the specimen volume. New randomly generated positions (coordinate triplets of particle centers) are accepted if they do not overlap with previously placed particles and volume boundaries. In this step zero-radius particles are also generated on boundary surfaces. 2) Delaunay tetrahedralization (Delaunay 1934; Barber et al. 1996) of the generated points defines the lattice mesh. Each tetrahedron edge represents a connecting strut between adjacent aggregate pieces. 3) A domain tessellation (dual complex of the Delaunay tetrahedralization) identifies the contact area through which the interaction forces are transmitted from one aggregate to the adjacent ones. The contact areas of the connecting struts are, in general, not planar (Fig. 1a). For sake of simplicity, the constitutive law is imposed on the projection of this area on a plane that is orthogonal to the connection (Fig. 1b) and that intersects the point located at midway of the strut counterpart belonging to the matrix (point O in Fig. 1c). The application point C of the interaction forces (contact point, Fig. 1b) is the center of mass of the projected area and, in general, does not lie along the axis connecting the particle centers. Rigid body kinematics describes the displacement field along the connection and the displacement jump uC at the contact point, divided by the length L of the connection, defines the strain vector. The components of the strain vector in a local system of reference, characterized by the unit vectors n, l, and m, are the normal and shear strains: εN = nT uC lT uC mT uC ; εL = ; εM = L L L (1) Figure 1: a) Contact area (potential crack surface) between two adjacent aggregates; b) projected contact area; c) geometry and degrees of freedom of a connecting strut; d) stresses acting on the contact area and related nodal forces. The unit vector n is orthogonal to the contact area and the unit vectors l and m are mutually orthogonal and lie in the contact plane. The normal stress σN and the shear stresses σL and σM (depicted in Fig. 1d) are computed through damage-like constitutive relations, which read σN = (1 − D)EεN 2 (2) defined as follows √ εN σN α tan ω = √ = σT αεT (6) The strain εT = ε2 + ε2 and the stress σT = M L 2 2 σM + σL are the total shear strain and stress, respectively. The function σ0 (ω) delimits the elastic domain as shown in Fig. 2a. The evolution of the boundary is governed by the initial slope K(ω) (Fig. 2b) which is defined as  nc Kc 1 − ω+π/2 ω ≤ ω0 ω0 +π/2 (7) K(ω) = nt Kt f (λ) ω−π/2 − 1 ω > ω0 ω0 −π/2 Parameters Kc and nc govern the nonlinear compressive and shear-compressive behavior (ω ≤ ω0 ) at the meso-level and are considered as material constants. For tensile dominated stress states the post-peak slope is made sensitive to the transverse confinement of the connection through the function f (λ), which reads 1 f (λ) = (8) 1 + −λ/λ0 where λ is the confinement strain (Cusatis et al. 2005). The characteristic strain parameter λ0 governs the sensitivity to the confining strain. If no confinement is present the meso-level fracturing behavior depends on the meso-level fracture energies Gt for mode I and Gs for mode II. In order to preserve the correct energy dissipation, the integral of the stress-strain curve for ω = 0 (pure shear) and for ω = π/2 (pure tension) must be equal to Gs and Gt , respectively, divided by the length of the connection. Since the stress-strain relation is exponential, we have −K(0) = Ks = 2αE Lcr /L − 1 s 2E Lcr /L − 1 t (9) (10) Figure 2: a) Elastic domain at the meso-level. b) Initial slope function. and σi = α(1 − D)Eεi (i = M, L) (3) where E is the elastic stiffness of the link, D is the damage parameter, α is the ration between the tangential and the normal stiffness of the link. We can write E = L/(La /Ea + Lm /Em ) if we assume that the aggregates and the matrix, represented by segments La = R1 + R2 and Lm , respectively (Fig. 1c), act in series. The elastic parameters Ea and Em characterize the elastic properties of aggregates and matrix. The damage parameter is defined as D = 1 − σ/Eε, 2 2 2 in which σ = [σN + (σM + σL )/α]1/2 is the effective 2 2 stress, and ε = [εN + α(εM + ε2 )]1/2 is the effective L strain. The evolution of the effective stress σ as function of the effective strain ε is governed by the equations σ = E ε, 0 ≤ σ ≤ σb (ε, ω) ˙ ˙ (4) −K(π/2) = Kt = and nt = ln[Kt /(Kt − Ks )] ln(1 − 2ω0 /π) (11) where σb (ε, ω) is an exponential strain-dependent boundary simulating damage, fracture and plasticity at the meso-level. Its analytical expression is σb (ε, ω) = σ0 (ω) exp K(ω) σ0 (ω) ε− σ0 (ω) E (5) 2 2 where Lcr = 2αEGs /σs and Lcr = 2EGt /σt are two s t characteristic lengths of the material. 3 in which • = max{•, 0}. The coupling strain ω is 3 INCORPORATION OF CREEP AND RATE EFFECT Strength of concrete depends upon the rate of loading. One of sources of the rate effect in concrete is the viscoelastic nature of the cement paste which is modeled in this paper through the microprestress-solidification theory (Baˇ ant et al. 2004). In Eq. 4 the hypo-elastic z constitutive relation σ = E ε should then be replaced ˙ ˙ by the differential equations which describe the evolution of an aging Maxwell chain or, equivalently, by the integral equation obtained through the compliance function J(t, t ) following from the theory. However, since only dynamic applications will be considered, it suffices to use just one term of the Maxwell chain, represented by the differential equation σ + σ/τ = E ε ˙ ˙ (12) a) 1.6 1.4 fc 1.2 ref fc dyn Constitutive Law inertia forces excluded Numerical inertia forces included Abrams Atcheley Watstein Dilger Bresler Rasch Hatano 1.0 0.8 -7 -5 10 10-6 10 10-4 10-3 10-2 10-1 1 10 where τ = η/E is the relaxation time. The elastic modulus E and the viscosity η must, once again, be computed taking into account the different properties of aggregates and embedding mortar. Assuming a series coupling between these two constituents, we get again E = L/(La /Ea + Lm /Em ) and η = ηm L/Lm . The viscoelastic properties of mortar Em , ηm can be obtained following Baˇ ant et al. (2000), in which the z Maxwell model is identified by requiring its compliance function to be tangent to the actual compliance function J(t, t ) at the time tch = tload /2 where tload is the duration of the dynamic event. From this condition one obtains Em = 1/(Jch − tch J˙ ) and ηm = ch ˙ where Jch = J(t + tch , t ), J˙ = J(t + tch , t ), ˙ 1/Jch ch t = age at loading. Rate effect due to creep (material visco-elasticity) is, however, insufficient because it cannot model reversal of softening into hardening after a sudden increase of the loading rate (Baˇ ant et al. 1995). This z phenomenon can be modeled only through a rate dependence of the fracture process. By using activation energy theory, Baˇ ant (1995) showed that in conz crete the rate dependence of crack opening can be expressed, at reference temperature (T = T0 ≈23◦ C), −1 as w = k0 sinh[kr Nb (w)−1 (σ − f 0 (w))] where k0 ˙ and kr are constants, f 0 (w) is the cohesive law observed for the rate of loading typical of the static material tests in the laboratory, and Nb (w) is the number of bonds across the cohesive crack that still have load-carrying capacity. Assuming Nb (w) ∝ f 0 (w) and solving the previous equation for σ, one gets σ(w, w) = [1 + k1 asinh(w/k0 )]f 0 (w) ˙ ˙ (13) Strain rate [sec ] -1 b) 4.8 4.2 3.6 ft ref ft 2.4 dyn Constituive Law inertia forces excluded inertia forces included Numerical Heilmann Hatano Takeda Kormeling et Al. Kvirikadze Sneikin Komlos Birkimer 3.0 1.8 1.2 0.6 10-9 10-7 10-5 10-3 10-1 -1 10 Strain rate [sec ] Figure 3: Rate effect on concrete strength. a) Uniaxial unconfined compression. b) Uniaxial tesnion. model can be calibrated by fitting experimental data (Veerheijm 1992) as shown in Fig. 3a for uniaxial compression and Fig. 3b for uniaxial tension. In representing the numerical results the applied peak stress is plotted along with the peak of the constitutive relation obtained neglecting the effect of the inertia forces. Both stresses are normalized with a reference strength in order to compare the results for different types of concrete. It is evident that for high strain rate the effect of the inertia forces cannot be neglected. This effect is much more important for the tensile strength than for the compressive one. This can be understood by writing fcdyn = fccr + f in (ε) ˙ and ftdyn = ftcr + f in (ε) ˙ (15) (16) which represents a vertical scaling of the static cohesive law by a function of the crack opening rate. In the present model, the rate dependence can then be introduced, adhering to similarity with Eq. 13, by scaling the stress-strain boundary 5 by a function of the strain rate ε ˙ σb = F (ε) σ0 (ω) exp ˙ K(ω) σ0 (ω) ε− σ0 (ω) E (14) where F (ε) = [1 + c1 asinh(ε/c0 )]. The presented ˙ ˙ 4 where fcdyn , ftdyn are the peaks of the applied stress, fccr , ftcr are the peaks of the stress obtained neglecting the inertia forces and f in (ε) is the inertia force abso˙ lute value which is the same in compression and in tension if the same specimen is used. If we normalize the Eqs 15 and 16 by fccr and ftcr respectively, assuming ftcr ≈ fccr /10, we get fcdyn f in (ε) ˙ =1+ cr cr fc fc and ftdyn f in (ε) ˙ f in (ε) ˙ =1+ ≈ 1 + 10 cr cr cr ft ft fc (17) a) (18) b) We see then that the percentage increment on the tensile strength due to the inertia forces is ten times the percentage increment on the compressive strength. 4 STEEL-CONCRETE BOND SIMULATION The behavior of reinforced concrete structures is greatly influenced by the behavior of the reinforcement and by its coupling with the surrounding concrete matrix. The most common type of reinforcement used in practice has the form of steel rebars. In structural dynamic simulations, rebars are effectively modeled using beam finite elements with strainrate dependent elasto-plastic constitutive equations for the steel material. Various algorithms have been developed for tying rebars elements to lattice model, as discussed in Pelessone (2005). The rebar-concrete coupling algorithm employed in the simulations described below is the simplest of them and it provides good coupling while simplifying the model generation. The lattice system and the beam finite element mesh simulating the reinforcement (rebars) are independently generated. The particles (aggregate + surrounding mortar) crossed by the rebars are then constrained to the rebar nodes through either a penalty method or a master-slave method. The two constraint algorithms have their own advantages and disadvantages. The penalty method is conceptually simpler and easier to implement. In addition, it also allows, contrarily to the master-slave approach, modeling nonlinear bond slippage at the concrete-rebar interface. The downside of the penalty method is that the rigid constraint holding in the linear elastic regime can be only approximately enforced by setting the penalty constants to very high values. In turn, these constants, which represent the stiffness of the concrete-rebar bond interface may lead to numerical instabilities. The numerical simulations presented in the next sections are carried out by adopting the master-slave approach. 5 NUMERICAL EXAMPLES This section deals with the numerical simulation of a projectile penetration in plain and reinforced concrete slabs carried out by using the object oriented explicit code MARS (Pelessone 2005). The numerical results are qualitatively compared to the experimental data by 5 c) d) Figure 4: a) Assemblage of aggregates. b) Beam element mesh simulating the slab reinforcement. c) Aggregate particles constrained to the rebars. d) Concrete slab and vertical projectile. Hanchak et al. (1992). The experiments are relevant to the impact of an ogive-nose steel projectile (25.4-mm caliber with total length of 143.7 mm) at velocities ranging from 300 to 1000 m/s. The impact was orthogonal to the target surface. The target consisted in prismatic concrete slabs with a depth of 178 mm reinforced by three layers of square grids of reinforcing steel bars with diameter of 5.69 mm. One layer was placed in the middle of the slab and the other two layers at 12.7 mm (0.5 in) from the surfaces. Each layer consists of a set of seven rebars in one direction and another set of seven rebars in the orthogonal direction. Fig. 4a shows the particle assemblage, consisting of 38448 particles, used to simulate the concrete slabs. Figure 5: Average velocity evolution: a) Vertical impact; b) Oblique impact on plain concrete slab; c) Oblique impact on reinforced concrete slab. Particle diameters range from 4.44 mm (0.175 in) to 6.68 mm (0.263 in) and occupy approximately 40% of the original volume of the concrete specimen. The finite element mesh modeling the rebar layers is shown in Fig. 4b. In the model, the crossing rebars are offset by one diameter like in the actual test specimen and therefore touch each other but do not overlap. The rebars are modeled using beam finite elements with three integration points over the cross-section at the 6 midpoint of the beam. Edge-edge contact conditions prevent rebars from penetrating each other. The edgeedge contact algorithm includes static and dynamic friction forces. Rebars finite elements are treated as cylindrical solids and contact begins when two cylindrical surfaces from two different elements start penetrating each other. A penalty function formulations is used to define the normal contact force as a function of penetration. The rebar-concrete interaction is implemented by constraining the slave particles that intersect the volume occupied by the cylindrical rebars to move with the closest master rebar elements. Obviously, a particle can only be slave to a single rebar element. Fig. 4c shows the 2266 particles that are linked to the rebars in the model. The steel projectile is modeled using 896 8-node hexahedral solid elements (Fig. 4d). Initially, the classical Flanagan-Belytschko formulation (Belytschko and Liu 2000) with a single integration point and hourglass stabilization was tried. However, because of the very large contact forces exerted on the bullet during impact, there was a tendency for the hourglass modes to grow excessively and induce unreasonable deformations in the bullet that would cause the simulation to stop. For this reason, we employed an 8integration point formulation which turned out to be much more stable and made it possible to complete the simulations. The steel bullet interacts with the reinforced concrete slab through two types of contact conditions: 1) particle-face contacts, where the discrete particles interact with the external faces of the bullet finite element mesh, and 2) node-rebar contacts, where the nodes of the bullet interact with the cylindrical surfaces of the rebars. The reinforced concrete target is assumed to be free in space with no initial velocity. The bullet is fired towards the target at different velocities and impact angles. Fig. 5a shows the evolution of the mean velocity of the projectile during the penetration of a projectile fired orthogonally to the slab surface. In the graph three curves are reported. One curve is relevant to the penetration through a plain concrete slab. In this case the projectile completely penetrate the slab (Fig. 6a) with a residual velocity of 359 m/s. The other two curves are relevant to the penetration of a reinforced slab. In one case (Curve labeled 1) the projectile penetrates the slab without hitting the reinforcement rebar grid (Fig. 6b) and the residual velocity is very close (345 m/s) to the preceding case of the plain concrete slab. This result is in agreement with the aforementioned experimental results, in which the measured residual velocity (relevant to a reinforced slab) was slightly higher than 350 m/s. In the second case (Curve labeled 2) the projectile is aimed at the intersection of two rebars and encounters much more resistance. As a result, the projectile is slowed signif- icantly (Fig. 6c) and its residual velocity (27 m/s at 1 ms) approaches zero. a) a) b) b) c) c) Figure 6: Deformed shape at 0.75 s during the vertical impact. a) Plain concrete target. b) Target hit not in correspondence of the reinforcement c) Target hit in correspondence of the reinforcement. Figs. 5b,c show the evolution of the vertical and horizontal components of the mean velocity of a projectile fired obliquely against a plain concrete slab and a reinforced slab, respectively. The inclination of the projectile is 40 ◦ respect to slab surface (Fig. 7a). In both cases the residual velocity approaches zero and the projectile does not go through the slab (Figs. 7b, c). In the case of reinforced concrete, however, the arrest of the projectile is more sudden and the damage caused by the impact is less severe. 6 SUMMARY AND CONCLUSIONS In this paper a previously developed constitutive model, the so-called Confinement-Shear Lattice (CSL) model has been enhanced through the introduction of the effect of creep and rate sensitivity of fracture processes. The new formulation has been calibrated by simulating the effect of strain rate on compressive and tensile strength. Further, the model has been implemented in the object oriented computer code MARS, which can handle effectively very com7 Figure 7: Deformed shape at 1 s during the oblique impact. a) Initial configuration. b) Plain concrete target. c) Reinforced target plex contact situations such those arising during penetration simulation. The CSL model has been coupled with elastic-plastic beam finite elements which simulate the effect of steel rebars. Numerical simulations of the impact of a steel projectile into plain concrete and reinforced slabs have been carried out and the results are very realistic and in agreement with experimental data gathered from the literature. The main conclusion that can be drawn from the present study is that the CSL model formulation fits well the highly complex computational framework needed for the simulation of impact, penetration and fragmentation of reinforced concrete structures. Of course, further work is needed to effectively calibrate and validate the rate-dependent part of the formulation and the steel-concrete bond behavior. 7 ACKNOWLEDGMENTS Financial support by the Engineering Research and Development Center (ERDC) is gratefully acknowledged. REFERENCES Barber, C., D. Dobkin, and H. Huhdanpaa (1996). The Quickhull algorithm for convex hulls. ACM Trans Math Softw 22(4), 469–483. http://www.qhull.org. Baˇ ant, Z., L. Cedolin, and G. Cusatis (2004). z Temperature effect on concrete creep modeled by microprestress-solidification theory. J of Engrg Mech, ASCE 103(6), 691–699. Baˇ ant, Z. P. (1995). Creep and damage in conz crete (J. Skanly and S. Mindess ed.)., Materials science of concrete IV, pp. 355–389. Westerville, OH: Am. Ceramic. Soc. Baˇ ant, Z. P., F. C. Caner, I. Carol, M. D. Adley, z and S. A. Akers (2000). Fracturing rate effect and creep in microplane model for dynamics. J. Engrg. Mech., ASCE 126(9), 962–970. Baˇ ant, Z. P., W.-H. Gu, and K. T. Faber (1995). z Softening reversal and other effects of a change in loading rate on fracture of concrete. ACI Mat. J. 92, 3–9. Belytschko, T. and W. K. Liu (2000). Nonlinear finite elements for continua and structures. Chichester, England: WILEY. Cusatis, G., Z. Baˇ ant, and L. Cedolin (2003a). z Confinement-shear lattice model for concrete damage in tension and compression: I. theory. J of Engrg Mech, ASCE 129(12), 1439–48. Cusatis, G., Z. Baˇ ant, and L. Cedolin (2003b). z Confinement-shear lattice model for concrete damage in tension and compression: Ii. computation and validation. J of Engrg Mech, ASCE 129(12), 1449–58. Cusatis, G., Z. Baˇ ant, and L. Cedolin (2005). z Confinement-shear lattice csl model for fracture propagation in concrete. Comput. Methods Appl. Mech. Engrg.. In press. Delaunay, B. (1934). Sur la sph` re vide. Bull Acad e Sci USSR(VII), Classe Sci Mat Nat, 793–800. Hanchak, S., F. M.J., E. Young, and J. Ehrgott (1992). Perforation of concrete slabas with 48 mpa (7 ksi) and 140 mpa (20 ksi) unconfined compressive strengths. Int J Impact Eng 12(1), 1–7. Pelessone, D. (2005, September). Discrete particle method. Technical report, ES3 Engineering and Softaware System Solutions, Inc., Solana Beach, CA 92075. Final report submitted to Enginnering Research and Development Center (ERDC). Contarct DACA42-02-C-0054. Veerheijm, J. (1992). Concrete under impact tensile loading and lateral compression. Ph. D. 8 thesis, Delft University of Technology, Prins Maurits laboratorium TNO.