Abstract
The formulations of the advanced concrete RHT model adopted in AUTODYN are investigated and numerical studies are conducted to study the RHT model’s actual performances under various loading conditions. It is found that using of default values in the RHT model is not able to simulate the realistic behavior of concrete under various loading conditions. Thus modified parameters in the RHT model are proposed to better capture the realistic behavior of concrete under such loading conditions. Furthermore, numerical simulation of normal concrete slabs and multilayer concrete slabs subjected to blast loading is conducted using AUTODYN with both the default and modified RHT parameters. Experimental readings from field blast tests are used to validate the numerical model developed. It is shown that the results from numerical simulations using the modified RHT parameters and the measurements from the field blast test agree well in terms of damage pattern, crater diameter, and acceleration. Hence, it can be concluded that the RHT model with modified parameters can capture the mechanical behavior of concrete structures well. The validated model can be further used to conduct a parametric study on the influence of key parameters (i.e., compressive strength, fracture energy, and thickness) on blast resistance of concrete structure.
1. Introduction
Around the world, one of the major construction materials used for both structural and infrastructural elements is concrete. At present, concerns about the safe living environment and the efficient protective structures and infrastructures are arising. The performance of concrete structures subjected to severe loading such as blast has been extensively studied in the last few decades [1–5]. Concrete is a brittle material which consists of cement paste, aggregates, and admixture (optional). The brittle behavior of concrete and other geomaterials, that is, rock and soil, presents obviously different strengths in compression and tension. Concrete also shows pressure hardening and strain hardening behavior under static loading, as well as strain rate hardening in tension and compression under dynamic loading before it fails. When concrete begins to fail, it gradually loses its loading carrying capacity, and this phenomenon is also called the strain softening [6].
Recently, dedicated research has been devoted to the development of reliable methods and algorithms to accurately analyze behaviors of structures and infrastructures subjected to dynamic loadings. Some techniques such as the explicit numerical analysis codes AUTODYN [7] and LSDYNA [8] are available for the modeling of interactive behavior between blast wave and concrete structures. In AUTODYN, different numerical techniques, such as pure Lagrangian and Lagrangian-Eulerian interaction algorithms, are available to analyze interactive response between the blast wave and concrete structure under blast loading [3, 9, 10].
There are many factors that will influence the reliability of numerical simulation. Among these factors, material model plays a key role because it reproduces the essential physical mechanisms of the material under various loading conditions. There are a number of material models for concrete developed in recent years, such as RHT model [1, 11], the Concrete Damaged Plasticity model (also known as MAT 72R3 model in LSDYNA) [11, 12], and HJC concrete model [13, 14]. These robust material models are capable of capturing varying concrete material behaviors under different loading conditions. When subjected to dynamic loading such as blast loading or high impact loading, concrete shows highly nonlinear response. Besides, due to the general complexity of the constitutive models, the determination of the parameters (i.e., residual strength, failure strain, and failure criteria in model) also plays an important role in achieving the actual performance of the concrete materials. This requires sufficient understanding of the modeling formulation and the associated considerations.
In this paper, the formation of the RHT model in hydrocode AUTODYN is carefully evaluated and examined through the numerical tests by generating stress-strain relationships of concrete under various stress conditions. Some potential issues of the simulations using the present form of the RHT model are highlighted and discussed. Within the framework of its current formulation, modified RHT parameters are recommended which are able to reproduce more realistic postpeak softening behavior of concrete under tension and compression. Furthermore, the RHT model with both the default and the modified parameter sets is applied to simulate a series of physical tests on a concrete slab and a multilayer concrete slab under blast loading. The simulation results are compared with the experimental observations. The numerical simulation results using modified RHT model are also compared with those using the MAT 72R3 model in LSDYNA. Conclusions are then drawn on the basis of both the material model exploration and the simulations of concrete slab subjected to the blast loading.
2. RHT Concrete Model in AUTODYN
The RHT concrete model is an advanced plasticity model for brittle materials developed by Riedel et al. [1] from Ernst Mach Institute. This model takes into account some robust features, such as pressure hardening, strain hardening, strain rate hardening, strain softening, and third invariant dependence. In stress space, three pressure-dependent surfaces, elastic limit surface , failure surface , and residual surface (as shown in Figure 1), are implemented to model material hardening and softening responses properly. The key formation of these three surfaces is discussed in the following sections, and the detailed description of this model can be further referred to in [1, 10, 11].
2.1. Elastic Limit Surface and Prepeak Loading Surface
The elastic limit surface dominates the elastic stresses. With the increase in stresses, strain hardening takes place until the failure surface is reached (as to in Figure 1). surface is obtained by scaling down from surface, which can be expressed using the following equation:where the scaling factor is the ratio of elastic tensile stress () or compressive stress () over the respective ultimate strength ( for ultimate tensile strength or for ultimate compressive strength) along a radial path. Considering the fact that concrete initiates the inelastic behavior under tension at around 50~80% and under compression at about 30% of the respective maximum strength [10, 15], in (1) varies linearly with pressure between the value related to uniaxial tension and that related to uniaxial compression [10]. The parabolic cap function is used to ensure the consistency between inelastic volumetric and deviatoric stresses [11].
The prepeak loading surface is subsequently defined through interpolation between the elastic and failure surface using the hardening slope. This can be seen in Figure 2 in the case of uniaxial compression, from which the equation can be expressed as where the definitions of and are shown in Figure 2. From the figure, it is shown that is the plastic strain before the failure strength; is the total plastic strain, and it can be determined by the secant modulus between the elastic limit surface and the failure surface.
2.2. Failure Surface
In the RHT model, the failure surface can be expressed as a function of hydrostatic pressure , Lode angle , and strain rate : where represents the compressive meridian, and it can be expressed as
In (4), denotes the uniaxial compressive strength of the material; and are two constant model parameters without dimension; is the hydrotensile stress normalized by the compressive strength ; represents the dynamic increase factor (DIF) as a function of the strain rate . Equation (3) indicates that the failure surface is established via rotation of the compressive meridian around the hydrostatic axis. By multiplying to the compressive meridian, the shear and tension meridians in the stress space can be taken into account in the failure surface. Thus, can be written as [15] in which / (refer to Figure 3). The Lode angle is a function of the second and third deviatoric stress invariants, and it can be obtained by
Figure 3 illustrates a typical shape of the deviatoric section plane of a strength surface. It should be noticed that for the concrete material the deviatoric section typically transits from a triangular shape at a low pressure (brittle material) to a circular shape at a high pressure (ductile materials) [15].
2.3. Residual Surface and Postfailure Surface
In order to represent the strength of the completely crushed concrete material, an independent residual strength surface is used in RHT model and can be expressed as [11]where and are two constant parameters without dimension, and they can be determined from curve fitting of the experimental data. These two parameters ( and ) determine the residual strength of concrete material, which will be shown in the following section.
After failure is initiated, a damage model is used for strain softening which considers the gradual loss of the load carrying capacity of concrete after reaching its ultimate tensile or compressive strength. This postfailure surface can be achieved by linear interpolation from the failure surface to the residual surface (as to in Figure 1 and (8)) and incorporation of a damage factor D (see (9)). One haswhere the damage parameter is defined as [10] in which is the normalized pressure variable and is the plastic strain increment. From the equation, it can be observed that the damage factor is dependent on the normalized pressure , the shape parameters and , and minimum failure strain (in Section 3, is classified into two types: one is for minimum tensile failure strain and the other is for minimum compressive failure strain). A proper selection of parameters , , and is crucial in order to obtain a reasonable postpeak softening behavior of the concrete material. The default values of these three parameters in RHT model are found to be problematic in some cases [9, 10], and hence modifications will be proposed in the following section.
3. Numerical Test
To test the actual crack softening behavior of the RHT model in tension and compression, single element tests are conducted to evaluate stress-strain relationship of the element under various stress conditions. The detailed single element simulation approach can be found in [10].
3.1. Tension Softening
In the RHT model, the hydrotensile failure criteria are the default tension failure criteria which means if the value of the tensile stress in an element falls below a specified limit, tensile failure is assumed to occur. At the same time, the tensile failure strain is predefined in the hydrotensile failure criteria (the default tensile failure strain ), which indicates that the failure tensile stress and failure strain are not changed with varying element sizes, and hence the value of the fracture energy is greatly mesh-dependent. These phenomena can be clearly explained by (10) in AUTODYN [7]:where is the tensile failure strain, is the tensile failure stress, is the fracture energy, and is the characteristic length of the element. In AUTODYN 3D, is the diameter of a sphere having the same volume as the three-dimensional element [7].
In (10), it is obvious that, for given tensile failure strain and tensile failure stress, the fracture energy of a bigger element size could be larger than that of a smaller element size. However, it is well known that the fracture energy is an important parameter for the simulation of tensile softening behavior of brittle material, and hence a fracture energy value is usually specified in the numerical model. In current study, the crack softening failure criteria [7, 17] are employed to model the tensile softening behavior with consideration of the fracture energy. The crack softening failure criteria adopt a linear softening law, in which, after the peak tensile stress , a constant strength degradation rate occurs with respect to the cracking strain. Hence the tensile failure stress () and the fracture energy () are given as input parameters in the crack softening failure criteria. During the cracking process, the discrete crack width () is smeared out over a certain distance, which is normally equal to the characteristic length of the element . The crack strain can also be calculated by (10).
Figure 4 reports the simulation results of the single element under various stress conditions adopting the default hydrotensile failure criteria and the crack softening failure criteria. After comparison between two failure criteria, it is found that using the default hydrotensile failure criteria tends to show an unrealistically descending process in the strain softening range for tension. For example, for a uniaxial tension test shown in Figure 4(a), the final failure strain is on the order of 0.01, which is more than one order of magnitude higher than the generally observed tensile failure strain from physical tests [15]. However, when adopting the crack softening failure criteria, the RHT model presents more realistic performance in the tensile failure strain under various stress conditions as shown in Figure 4. This modification is also more suitable for simulation of damage pattern of concrete materials under blast loading which will be shown in another application example in Section 4.
(a)
(b)
(c)
3.2. Compression Softening
Shear cracks usually appear after the yield limit is reached. Beyond the yield limit, the concrete exhibits softening behavior with yield stress and stiffness degradations. Such plastic behavior of concrete, to be more specific, the compression softening behavior, has to be considered into the constitutive model. Figure 5(a) shows that the concrete adopting the default RHT parameters possesses a higher residual strength of 21 MPa in the postpeak domain which is unreasonable in the low pressure regime. To make the concrete model more realistic, two steps have been taken to obtain reasonably the residual strength and the failure strain. Firstly, the residual yield stress should be changed [5]. From the formula of the residual surface, it is found that both and control the final residual strength. After a few trials, the values for and of 1.1 and 0.9, respectively, are proposed to replace the default values of 1.6 and 0.61, respectively. With such modification, the residual yield stress decreases significantly, and this improves the concrete material behavior in the low pressure regime. It is also observed that the compressive failure strain of concrete using the default RHT model parameters is around 0.01 as shown in Figure 4(a), which is much larger than the realistic value. Thus, the following step is used to adjust the compressive failure strain and the damage parameter . Based on (9), it is found that compressive failure strain affects the total failure strain under the lower pressure (≤). Therefore, is adjusted to be 0.001 which gives a total strain value of around 0.1% in uniaxial compression and is more reasonable compared to the experimental results [6]. The damage parameter controls the rate of accumulated damage. For the failure behavior under compression, the parameter controls the shape of the softening branch after the peak strength. In this study, is changed to 0.02 and remains unchanged based on numerical tests.
(a)
(b)
(c)
Figures 5(a) and 5(b) present the stress-strain curves of a single element in uniaxial and biaxial compressive tests based on numerical simulation using modified RHT parameters. It is observed that the material with the modified RHT parameters can properly reproduce the behavior of the concrete under both uniaxial and biaxial compressions. Figure 5(c) shows the stress-strain curve of the triaxial compression tests with different confinement pressures based on numerical simulation adopting the modified RHT parameters. It is observed that the peak strength is enhanced with the increase of the confinement. This represents a realistic behavior of concrete in triaxial compression as observed in experimental data [15].
4. Validation of the Modified RHT Model
In order to further investigate the applicability of the RHT model with modified parameters in the analysis of concrete structures under severe dynamic loading, the RHT model is used to simulate full scale field blast tests on a normal concrete slab and a multilayer concrete slab. The simulation results are then compared with the experimental observations. A brief introduction of the field blast tests of the normal concrete slab and the multilayer concrete slab is given below and more details can be found in [16].
4.1. Field Blast Test Configuration and Instrumentation
Two samples are tested in the field blast test, that is, one normal concrete slab sample and one multilayer concrete slab sample. The details of the cross sections of these two slabs are given in Table 1. The multilayer concrete slab is made of Asphalt Concrete (AC), High Strength Concrete (HSC), and Engineered Cementitious Composites (ECC) from the top to bottom. Each slab is subjected to one blast detonation. A bomb equivalent to 7.3 kg TNT charge is placed with its center of gravity about 170 mm above the center of the slab surface. These two samples are cast at site with the dimension of 2.8 m in both length and width and 0.275 m in thickness.
Figure 6 shows the completed normal concrete slab and multilayer concrete slab with four anchors each. The anchors are used to fix the slab to the ground. Standard tests of material properties are conducted for each material cast, and the results are summarized in Table 2. The geogrid used to reinforce the AC layer in this study is the Polyfelt Microgrid MG-100 with bidirectional tensile strength of 100 kN/m and it has an aperture size of 7 mm.
(a) Completed normal concrete slab
(b) Completed multilayer concrete slab
Various instruments are installed onto the slabs to measure the responses of the pavement slabs during the blast. Figure 7 shows the instrumentation installed on the slab. Four accelerometers are installed at the middle of each side of the slab to measure the vertical (V1 and V2 in Figure 7) and horizontal acceleration (H1 and H2 in Figure 7). The accelerometers are mounted on steel frames casted together with the slab. The measurement results of the field blast tests will be discussed and compared with the numerical simulation results in the following section.
4.2. Validation 1: Normal Concrete Slab Subjected to Blast Load
4.2.1. Numerical Model
A 3D numerical model is established to model the normal concrete slab subjected to blast loading as shown in Figure 8. A quarter of the slab is modeled considering symmetry. In the model, the initial detonation and blast wave propagation are modeled using an axis-symmetric 2D model. Before the blast wave reaches the slab, the blast wave distribution is remapped into a 3D model. The foundation soil has a dimension of 825 mm in thickness and 2800 mm in length and width. Eight-node hexahedron Lagrange element with one-point gauss integration is used to represent the concrete slab and foundation soil. The interaction between the slab and foundation soil is considered using a contact algorithm. In AUTODYN, the penalty-based contact is used to detect the penetration between the “slave” surface and “master” surface. When a penetration is found, a force proportional to the penetration depth is then applied to resist the penetration. Thus the interface force could be calculated based on the elastic spring theory [7].
The air and blast waves are modeled using Euler elements. The blast loading acting on the concrete slab can be achieved by performing the interaction between the Lagrange elements and the Eulerian elements. The Lagrange elements apply geometric constraints on the Eulerian elements, while the Eulerian elements provide a pressure boundary to the Lagrange elements. In order to make Lagrangian and Eulerian elements interact effectively, it is important for the Eulerian element size to be equal to or smaller than the Lagrangian element size. Hence, mesh convergence analysis using the sizes of 9 mm, 8 mm, 7 mm, and 6 mm and 5 mm, 4 mm, and 3 mm for Eulerian elements is performed in current study. The size of Lagrangian elements is fixed to be 10 mm for the convergence study. Note that 3 mm is the smallest possible size for Euler elements using the computer resource available. During the convergence study, the pressure and impulse exerted on the center part of the slab are monitored and compared. After the convergence study, it is found that the size of 4 mm for Eulerian elements and 10 mm for the Lagrangian elements in current model is preferred considering the optimum CPU time and computer resource.
The boundary condition for air is outflow which allows the air to flow out. For the foundation soil, transmit condition is employed to simulate the semi-infinite space. The anchors on the slab are simulated as fixed points in the 3D model. The RHT model with modified parameters is employed to simulate the concrete slab subjected to blast loading, and an elastic-plastic Drucker-Prager model [15] is employed to simulate the foundation soil. The parameters of the air, foundation soil, and normal concrete are given in Tables 3–5, respectively.
4.2.2. Results and Discussion
Figures 9(a) and 9(b) show the damage patterns of the concrete slab in the numerical model using default parameters and those using the modified RHT parameters in comparison with those in the field blast test as shown in Figure 9(c). Much smaller craters and less severe damage cracks are observed in the numerical simulation using the default RHT parameters than those observed after the field test. Such findings indicate that using the default RHT parameters could largely underestimate the damage pattern of the concrete slab under blast loading in terms of the cracks and the crater diameter. The concrete slab simulated using default RHT parameters appears to behave in a substantially ductile manner. This echoes the findings in the single element tests using the RHT model with the default parameter values. When adopting the modified RHT parameters, the damage and failure models of the concrete slab in the numerical analysis are comparable with the experimental results in terms of the crater diameter and the crack pattern as shown in Figures 9(b) and 9(c).
(a) Damage pattern with default RHT model
(b) Damage pattern with modified RHT model
(c) Physical damage of normal concrete slab after blast event
The damage pattern of the normal concrete slab simulated using the modified RHT parameters is then compared with that using MAT 72R3 model in LSDYNA [8]. The detailed numerical model using the MAT 72R3 material model in LSDYNA can be referred to in [16]. Numerical simulation results using the MAT 72R3 model are given in Figure 10. Figures 9(b) and 10 indicate that the damage patterns of the concrete slab simulated using both material models are almost the same. However, the crater diameter predicted by the MAT 72R3 model is smaller than that by the modified RHT model and that observed after the field blast test. This is possible that different blast simulation method and different mesh sizes are adopted for the numerical model in LSDYNA. Considering that there are 20 parameters that need to be determined in the MAT 72R3 model which makes it time consuming to establish a reliable numerical model, the RHT model with modified parameters is more suitable for the simulation of concrete structures under blast loading for quick damaged assessment.
Table 6 reports the results of vertical and horizontal acceleration of the concrete slab obtained from both the field measurement and numerical model. It is found that the variation of the vertical acceleration between the values from the field test and numerical simulation is around 6.8%, and the numerical simulation predicts a higher vertical acceleration than the value measured in the field blast test. However, in view of the inherent variation in the field blast tests, a prediction of 6.8% deviation of the numerical simulation results from field blast test results can be considered as acceptable. For the horizontal acceleration, two different readings were obtained from the accelerometers because the charge is not perfectly placed at the center above the slab. It can be observed from the table that horizontal acceleration predicted in the numerical simulation is right the average of the two measured acceleration values.
4.3. Validation 2: Multilayer Concrete Slab Subjected to Blast Load
4.3.1. Numerical Model
The dimension and mesh size of multilayer concrete slab in numerical simulation are the same as those in the normal concrete slab. The parameters used for the air and foundation soil are the same as those given in Tables 1 and 2. The numerical model for the multilayer concrete slab is shown in Figure 11. The thickness of each material is consistent with that used in the field blast test as stated in Section 4.1, that is, 100 mm for the ECC layer, 100 mm for the HSC layer, and 75 mm for the reinforced asphalt layer. The boundary condition and the location of detonation are the same as those in the case for normal concrete slab model, which are detailed in Section 4.2.
Laboratory tests [18] indicate that the geogrid enhances tensile strength of a material if the geogrid is added. Thus, in current numerical model, the asphalt layer is modeled with greater stiffness and tension strength in order to take into account the inclusion of the geogrid. The parameters used for the asphalt material are given in Table 7. For High Strength Concrete (HSC), the approach of the RHT model with modified parameters is used to simulate the material with consideration of high compressive strength. The RHT model with modified parameters is also employed to simulate the ECC material with consideration of its high compressive strength and superior ductile property. The parameters of the HSC and ECC material are given in Tables 8 and 9, respectively.
4.3.2. Results and Discussion
The results of damage patterns of the multilayer concrete slab after the blast in the numerical simulations are shown in Figures 12(a) and 12(b). Figures 12(c) and 12(d) present the damage patterns after the field test. In the field test, the asphalt layer breaks into a few fragments after the blast, while in the numerical simulation, only a centralized crater is produced and no fragment is observed. This could be because the Drucker-Prager model may not be suitable for the simulation of fragmentation. Further, it should be noted that the AC layer is mainly destroyed by the combination of the high temperature and blast pressure as observed in the field blast test. An advanced material model for asphalt with consideration of the temperature effect should be explored in future study.
(a) Damage contour of asphalt layer
(b) Damage contour of HSC layer with modified RHT model
(c) Physical damage of asphalt layer
(d) Physical damage of HSC layer
After removing the asphalt layer, the damage pattern of the HSC layer in the field test is shown in Figure 12(d). The damage pattern of the HSC layers in the numerical simulation is illustrated in Figure 12(b). The diameter of the crater is around 0.7 m in the field test, and it is 0.75 m in the numerical simulation. These two results are close to each other. The damage pattern of the HSC layer simulated using the modified RHT parameters is also compared with that using MAT 72R3 model in LSDYNA as shown in Figure 13 [16]. It is observed that the results from the simulations using both material models are comparable in terms of the damage pattern and crater diameter. Hence, the modified RHT model can represent the real behavior of the concrete slab under blast loading.
Table 10 reports the results of the vertical and horizontal acceleration of the multilayer concrete slab obtained from both the field measurement and numerical simulation. The vertical acceleration from the numerical simulation is quite close to that measured after the field test with a small deviation of 5.3%. The two horizontal acceleration readings on two sides are not the same because the charge is not perfectly located at the center above the slab.
5. Conclusion
This paper intensively studies the equation of the RHT model in the hydrocode AUTODYN. It is found that using the default RHT parameters such as hydrotensile failure criteria, residual strength parameters, and failure strain cannot predict the realistic behavior of concrete material under various loading conditions. Through numerical studies on behavior of concrete material in various loading conditions, the key parameters such as tensile failure criteria, residual strength parameters, and failure strain parameter controlling the behavior of concrete under uniaxial, biaxial, and triaxial tension/compression are discussed, and the modified RHT parameters are proposed to improve the model. Furthermore, numerical simulations of a concrete slab and a multilayer concrete slab under blast loading are conducted using AUTODYN with both the default and modified RHT parameters. The numerical simulations are validated by field blast test. Based on the numerical simulations, the incorporation of the modified RHT parameters reproduces more realistic structural behavior under blast loading than that of the default RHT parameters. The results in numerical simulations using RHT model with modified parameters agree well with the corresponding field test results including the damage pattern, crater diameter, and acceleration readings. The numerical model using the RHT model with modified parameters is then compared with that of MAT 72R3 in LSDYNA for the normal concrete slab and multilayer concrete slab. It is observed that the predicted damage pattern and crater diameter using the RHT model with the modified parameter and using the MAT 72R3 model show a better comparison with the experimental result. Considering that 20 parameters are needed in MAT 72R3 model which makes it time consuming to establish the numerical model, the RHT model with modified parameters is preferred for the simulation of the concrete structure under blast loading for quick damaged assessment. Hence, it can be concluded that the current numerical model using AUTODYN with modified RHT model could represent more realistic behavior of concrete structures under blast loading. The numerical model developed in this paper can be further used to conduct a parametric study on the influence of key parameters (i.e., compressive strength, fracture energy, and thickness) on blast resistance of concrete structure.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
Part of this research is sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry of China (Grant no. E262021514), Youth Teacher Training Scheme from Shanghai Education Committee (Grant no. ZZGCD15053), and foundation from Shanghai University of Engineering Science (Grant no. E10501140170).