Evaluation of automatic versus material test-based calibrations of concrete models for ballistic impact simulations

Abstract

Concrete is used for protective structures all over the world. Accurate response estimates to a given threat is vital for designing such structures. Concrete models often require numerous input parameters for which sufficient experimental data can be challenging to obtain. Some models are accompanied by parameter generators which use the unconfined compression strength to extrapolate the remainder of the parameters based on experimental databases. This study investigates simulation of ballistic impact on high-strength concrete with 75 MPa nominal unconfined cylindrical compressive strength. The first objective is to investigate the accuracy parameter generators to produce input data for commonly used concrete material models. The second objective is to establish and evaluate a simplified parameter calibration procedure based on standard material experiments and data from the literature. The results employing parameter generators varied notably between the models while still giving decent ballpark estimates. The parameters obtained from inverse modelling of standardized material tests improved the results significantly. The findings of this study recommend caution when using automatic parameter generators. Although a detailed calibration of these concrete models is complicated, a simplified calibration gives reasonable predictions, making this the advisable approach for designing concrete protective structures.

Introduction

Concrete is the most used construction material in the world. When weight and space are not constraints, concrete is the preferred option for critical infrastructure, nuclear reactors, protective structures and general fortification (Corbett et al., 1996). Because of its ubiquitous use, there are numerous studies on the behaviour of concrete subjected to extreme loadings, for example, vehicle collision (Do et al., 2019; Yonten et al., 2005), explosions (Kristoffersen et al., 2021a; Silva and Lu, 2009) and ballistic impact (Frew et al., 1998; Hanchak et al., 1992). Studies of this kind are in general divided into experimental investigations, analytical/empirical models, and numerical simulations, or any combination thereof. While experimental data is still of vital importance, and analytical/empirical approaches are still used (Warren et al., 2014; Yankelevsky and Feldgun, 2020), more and more studies are nowadays utilizing numerical simulations to study the behaviour of concrete exposed to extreme loadings. Such simulations are now accepted tools for evaluating blast or impact loads against materials and structures (Abdel-Kader, 2018). They are particularly useful where experiments are not feasible, for instance due to economic, spatial or practical restraints. If the simulations are to provide useful insight into the problem at hand, the constitutive relations describing the material behaviour of the concrete must be able to accurately account for all phenomena relevant to the specific load case and structure. In addition, the material parameters and numerical setup should be validated against material tests and/or component level tests so that the model may be used to make assessments for structures of a larger scale (Kristoffersen et al., 2019).

Numerical modelling of concrete is challenging due to the brittleness, inherent heterogeneity, microcracking, and in general, complex behaviour of the material (Lu et al., 2010). Nevertheless, there are good, homogenized models which can give reasonable macroscopic predictions (Cui et al., 2017), and it is models of this kind that is the concern of this study. Such models come with different features, limitations, and levels of complexity (Tu and Lu, 2009; Yonten et al., 2005), and choosing the appropriate material model for the task at hand is not always straightforward. Some models require only a few input parameters, while other more complex models require numerous parameters. Given a thorough calibration of the material input based on appropriate test data, the accuracy should in general improve for any material model. Producing this input for the relevant concrete can be an arduous ordeal, so several concrete models come with generators for the input parameters, or they use empirical relations to generate the required input. Parameter generators typically use the unconfined compressive strength $f_{\mathrm{c}}$ to extrapolate the parameter set from an existing experimental database. This can be a very useful tool for making ballpark estimates but should be used with caution. Different concrete models using the same $f_{\mathrm{c}}$ may elicit different responses for what should be the same concrete (Feldgun and Yankelevsky, 2019). Still, experimental observations of impact against massive semi-infinite concrete targets have shown that $f_{\mathrm{c}}$ is the primary variable for the target’s resistance to deep penetration (Forrestal et al., 1988; Frew et al., 1998; Li et al., 2005; Luk and Forrestal, 1987). For perforation of concrete slabs with a finite thickness, however, $f_{\mathrm{c}}$ has less predictive power. Hanchak et al. (1992) found that tripling the unconfined compressive strength resulted in a very modest increase in perforation resistance for 178 mm thick slabs impacted by 25.4 mm diameter ogive-nose projectiles. Similar results have since been replicated (Børvik et al., 2007), and recent experimental and numerical studies show that $f_{\mathrm{c}}$ alone has limited influence on the ballistic perforation resistance of thin concrete slabs (Kristoffersen et al., 2021b). This study is thus limited to ballistic impact against concrete slabs with a finite thickness where the unconfined compressive strength is one of several governing parameters. Quantities that typically correlate with $f_{\mathrm{c}}$ (for instance, Young’s modulus) may still affect the results to some extent.

Today, some of the most widely used concrete material models which include parameter generators are the Karagozian and Case (K&C) concrete model - Release III (Magellanes et al., 2010; Malvar et al., 1997), the continuous surface cap model (CSCM) (Murray, 2007; Zhao et al., 2021), and the Riedel–Hiermaier–Thoma (RHT) model (Grunwald et al., 2017; Riedel et al., 2009). All three models are readily available as standard models in the explicit finite element (FE) code LS-DYNA (LST, 2020). The Holmquist–Johnson–Cook (HJC) model (Holmquist et al., 1993; Johnson et al., 2018; Rajput et al., 2018) is another widely used model, and a modified version has been included due to its relative ease of calibration (Polanco-Loria et al., 2008). Numerous constitutive concrete models may be employed to simulate the ballistic impact resistance of concrete, for example, the Winfrith Concrete (WC) model (Broadhouse and Neilson, 1987), Concrete-Damage-Plasticity (CDP) (Lubliner et al., 1989) and a damage model based on effective rate-dependency proposed by Pereira et al. (2018). Alternatively, mesh-free methods are extensively used such us smoothed particle hydrodynamics (SPH) (Gingold and Monaghan, 1977), lattice discrete particle method (LDPM) (Cusatis et al., 2011) and discrete element method (DEM) (Antoniou et al., 2018).

The goal of this study is then to investigate the behaviour of commonly used concrete material models when applied to perforation of finite thickness concrete slabs by ogive-nosed steel projectiles. First, the models – those being the K&C concrete model, the CSCM and the RHT model – were used with their default (or generated) parameter sets, and FE simulations of material and component tests were carried out. Next, the material models were calibrated based on a set of standardized material tests, those being the cylinder compression test and the tensile splitting (Brazilian) test. Arguably, triaxial tests, dynamic tests, etc., should be included for a complete characterization of the material, but such tests are rare. The equipment required to complete the entire collection of tests is extensive. Thus, this work aims to investigate to which extent calibrations based on common material tests and some data from the literature can improve the numerical results compared with the default values and assess whether it is worth the effort. Several studies of similar nature are available in the literature (Abdel-Kader, 2018; Cui et al., 2017; Feldgun and Yankelevsky, 2019; Jiang and Zhao, 2015; Markovich et al., 2011; Ren et al., 2017; Tu and Lu, 2009; Wu et al., 2012; Yonten et al., 2005), but this work compares not only the different models with each other, but also the default calibration with a calibration produced on instrumented material tests and validation of the various models against ballistic impact test results. This calibration is conducted for a commercially produced concrete with a nominal unconfined compressive strength $f_{\mathrm{c}}$ of 75 MPa, dubbed in the following C75. It turns out that the calibrations based on the material tests, laborious as they are, do significantly improve the concrete perforation simulation results.

Concrete models

The concrete models under evaluation; Modified Holmquist–Johnson–Cook (MHJC), K&C - Release III, CSCM and RHT are all isotropic and elasto-plastic, and they use a yield surface to separate the elastic and plastic domains. The yield surface follows a different expression for each model, though all are functions of the hydrostatic pressure, the equivalent shear stress and the Lode angle.

The first invariant $I_1$ of the stress tensor $\mathbf{\sigma }$ characterizes the hydrostatic pressure $P$ (defined as positive in compression), that is,

$P=-\frac{I_1}{3}=-\frac{1}{3}\mathrm{t}\mathrm{r}\left(\mathbf{\sigma }\right)$

(1)

The second deviatoric invariant $J_2$ of the stress deviator tensor $\mathbf{s}$ describes the equivalent shear stress $q$ as

$q=\sqrt{3J_2}=\sqrt{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\mathrm{t}\mathrm{r}\left({\mathbf{s}}^{\mathbf{2}}\right)}$

(2)

The Lode angle $\theta$ is a function of the equivalent shear stress $q$ and the third deviatoric invariant $J_3=\left|\mathbf{s}\right|$ . It signifies the stress state position in the deviatoric plane and is given as

$\theta =\frac{1}{3}\arccos \left(\frac{27J_3}{2q^3}\right)$

(3)

Dependency on the stress triaxiality is established by a reduction factor $\mathcal{R}$ as a function of the Lode angle $\theta$ and determines the radial distance of any stress state from the hydrostatic axis. The function $\mathcal{R}$ is defined in the polar range $0°\le \theta \le 60°$ , and extends to all directions using threefold rotational symmetry. Thus, the angles $\theta =60°$ , $180°$ and $300°$ represent the shear strength in triaxial compression (compressive meridian), while $\theta =0°$ , $120°$ and $240°$ give the triaxial extension (tensile meridian). The expression for $\mathcal{R}$ typically permits a smooth transition in the deviatoric plane from a triangular shape with smooth corners at low pressures to a circular shape at high pressures as shown in Figure 1.

Figure 1.

Lode angle and yield stress in the deviatoric plane.

Therefore, the yield surface represents the deviatoric elastic limit in the principal stress space with a certain degree of symmetry around the hydrostatic axis as illustrated in Figure 2. The yield function for all the concrete models in the current study is expressed as a function of $P$ , $q$ and $\theta$ , that is,

${\mathcal{F}}_y=f\left(P,q,\theta \right)$

(4)

Figure 2.

Yield surface in the principal stress space.

Further, the three concrete models K&C, CSCM and RHT include an additional strength limit, the maximum surface. The maximum surface involves plastic hardening after yielding until reaching the limit strength. All four models employ a residual surface for the final strength after irreversible damage.

The volumetric pore collapse evolution differs for each model in terms of hydrostatic pressure-volumetric strain equations. In general, the pressure-porosity constitutive model consists of three phases. These are as follows: (1) fully elastic behaviour, (2) compaction, that is, pores collapse with plastic deformation and (3) solidification, or fully compacted. Thus, irreversible closure of the porosity is observed at the final phase. Then, the volumetric response returns to an elastic behaviour while retaining the irreversible deformation from the second phase.

Likewise, the strain-rate sensitivity is introduced with a different expression for each model. Therefore, the strength of all the models enhances with the increase of the strain rate.

Lastly, all four concrete models include damage accumulation, which has not been evaluated in the present study for brevity. For a more detailed description of the concrete models, the readers are referred to Grunwald et al. (2017) for RHT, Murray (2007) for CSCM, Malvar et al. (1997) and Magellanes et al. (2010) for K&C, and Polanco-Loria et al. (2008) for MHJC.

Modified Holmquist–Johnson–Cook model

The MHJC yield surface is given by

${\sigma }_{\mathrm{e}\mathrm{q}}^{*}=\{\begin{array}{c}B{\left[P^{*}+T^{*}\left(1-D\right)\right]}^{N}F({\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{*})\mathcal{R}\left(\theta ,e\right)\le S_{\max },P^{*}\ge -T^{*}\left(1-D\right)\\ \hspace{0.17em}0,P^{*}<-T^{*}\left(1-D\right)\end{array}$

(5)

in which ${\sigma }_{\mathrm{e}\mathrm{q}}^{*}=\frac{{\sigma }_{\mathrm{e}\mathrm{q}}}{f_{\mathrm{c}}}=\frac{q}{f_{\mathrm{c}}}$ is the normalized equivalent stress and $f_{\mathrm{c}}$ is the quasi-static uniaxial compressive strength. Note that $f_{\mathrm{c}}$ is often denoted $f_{\mathrm{c}}^{\prime }$ in relevant literature, for example, Malvar et al. (1997) and the CEB-FIB mode code (1990). Further, $P^{*}=P/f_{\mathrm{c}}$ is the normalized hydrostatic pressure, $P$ is the actual hydrostatic pressure and $T^{*}=T/f_{\mathrm{c}}$ is the normalized hydrostatic tension. The constant $T$ indicates the maximum tensile pressure $P$ as the initial intersection of the yield surface on the axis $P$ . In addition, $B$ and $N$ are the pressure hardening coefficient and exponent, whilst $D$ is the total damage and $S_{\max }$ is the normalized maximum available strength.

The reduction factor $\mathcal{R}\left(\theta ,e\right)$ describes the dependence of the yield surface on the Lode angle $\theta$ and the shape factor $e$ , and is defined as proposed by Willam and Warnke (1975). The shape of the yield function in the deviatoric plane transitions from a triangular shape with smooth corners at low pressure to a circular shape at high pressure, that is,

$\mathcal{R}\left(\theta ,e\right)=\frac{2\left(1-e^2\right)\cos \theta +\left(2e-1\right)\sqrt{4\left(1-e^2\right){\cos }^2\theta +5e^2-4e}}{4\left(1-e^2\right){\cos }^2\theta +{\left(1-2e\right)}^2}$

(6)

The normalized shape factor $e$ gives the strength ratio of the tensile to compressive meridian with respect to the normalized hydrostatic pressure $P^{*}$ and the reference pressure $P_{\mathrm{r}\mathrm{e}\mathrm{f}}^{*}$ , which is set to 10 internally. The normalized shape factor $e$ is determined as

$e=\{\begin{array}{c}0.65\\ 0.65+\left(1-0.65\right)\left(\frac{P^{*}}{P_{\mathrm{ref}}^{*}}\right)\\ 1.0\end{array}\begin{array}{c},\\ ,\\ ,\end{array}\begin{array}{c}{P}^{*}<0\\ \hspace{0.17em}0\le P^{*}\le \\ \hspace{0.17em}{P}^{*}>P_{\mathrm{r}\mathrm{e}\mathrm{f}}^{*}\end{array}{P}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{*}$

(7)

The pressure-porosity expression of the MHJC model constructs the three phases of the pore collapse evolution based on the crush volumetric strain ( ${\mu }_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}})$ and the solidification volumetric strain $({\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}$ ) according to

$P=\{\begin{array}{c}\begin{array}{cc}{\hspace{0.17em}K}_0\hspace{0.17em}\mu ,& \hspace{0.17em}\mu <{\mu }_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}\end{array}\\ \begin{array}{cc}\hspace{0.17em}{K}_{\mathrm{a}\mathrm{v}}\hspace{0.17em}\left(\mu -{\mu }_{\mathrm{p}}\right),& \hspace{0.17em}{\mu }_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}\le \mu \le {\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}\end{array}\hspace{0.17em}\\ \begin{array}{cc}{K}_1\hspace{0.17em}\overline{\mu }+K_2\hspace{0.17em}{\overline{\mu }}^2+K_3\hspace{0.17em}{\overline{\mu }}^3,& \hspace{0.17em}\mu >{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}\end{array}\hspace{0.17em}\end{array}$

(8)

Here the elastic phase is controlled by the bulk modulus $K_0=\frac{P_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}}{{\mu }_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}}$ , where $P_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}$ is the crush pressure. During the second phase the average bulk modulus $K_{\mathrm{a}\mathrm{v}}$ varies with the plastic volumetric strain ${\mu }_{\mathrm{p}}$ . $K_{\mathrm{a}\mathrm{v}}$ is computed internally by $K_{\mathrm{a}\mathrm{v}}=K_0(1-\frac{{\mu }_{\mathrm{p}}}{{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}})+K_1\frac{{\mu }_{\mathrm{p}}}{{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}\left(1-{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}\right)}$ with bounds $K_0$ and $K_1$ . Here, $K_1$ is the solidified bulk modulus at the end of the compaction phase and corresponds to the solidification pressure $P_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}$ . The solidification phase is governed by a cubic function depending on a modified volumetric strain $\overline{\mu }=\left(\mu -{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}\right)/\left(1+{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}\right)$ , where $K_1$ , $K_2$ and $K_3$ are the solidification function coefficients.

The associated plastic flow rule is employed and the plastic rate of deformation tensor $D_{ij}^{\mathrm{p}}$ is obtained as

$D_{ij}^{\mathrm{p}}={\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}\frac{3}{2}\frac{s_{ij}}{{\sigma }_{\mathrm{e}\mathrm{q}}}$

(9)

The equivalent plastic strain rate ${\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}$ is found by multiplying equation (9) with itself $D_{ij}^{\mathrm{p}}{D}_{ij}^{\mathrm{p}}={\left({\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}\right)}^2{(\frac{3}{2})}^2\frac{s_{ij}{s}_{ij}}{{\sigma }_{\mathrm{e}\mathrm{q}}^2}$ and combining with equation (2), that is,

${\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}=\sqrt{\frac{2}{3}{D}_{ij}^{\mathrm{p}}{D}_{ij}^{\mathrm{p}}}$

(10)

The strain-rate sensitivity is described as follows

$F({\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{*})={(1+{\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{*})}^C$

(11)

in which $C$ is the strain-rate sensitivity exponent, ${\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{*}={\dot{\epsilon }}_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}/{\dot{\epsilon }}_0$ is the normalized strain rate, and ${\dot{\epsilon }}_0$ is a reference strain-rate.

The total damage $D$ is a combination of the shear damage $D_{\mathrm{S}}$ and the compaction damage $D_{\mathrm{C}}$ , and reads as

$\left(1-D\right)=\sqrt{\left(1-D_{\mathrm{S}}\right)\left(1-D_{\mathrm{C}}\right)}$

(12)

The evolution of the shear damage ${\mathrm{\Delta }D}_{\mathrm{S}}$ and the compaction damage ${\mathrm{\Delta }D}_{\mathrm{C}}$ are calculated as

${\mathrm{\Delta }D}_{\mathrm{S}}=\frac{{\mathrm{\Delta }\epsilon }_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}}{{\epsilon }_{\mathrm{f}}^{\mathrm{p}}},{\mathrm{\Delta }D}_{\mathrm{C}}=\frac{\mathrm{\Delta }{\mu }^{\mathrm{p}}}{{\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}}$

(13)

where ${\mathrm{\Delta }\epsilon }_{\mathrm{e}\mathrm{q}}^{\mathrm{p}}$ is the equivalent plastic strain increment, ${\epsilon }_{\mathrm{f}}^{\mathrm{p}}=\alpha {\left(P^{*}+T^{*}\right)}^{\beta }\ge {\left({\epsilon }_{\mathrm{f}}^{\mathrm{p}}\right)}_{\min }$ is the plastic strain to fracture with $\alpha$ and $\beta$ as material constants and ${\left({\epsilon }_{\mathrm{f}}^{\mathrm{p}}\right)}_{\min }$ is introduced to prevent fracturing from low magnitude tensile stress waves. Equivalently, $\mathrm{\Delta }{\mu }^{\mathrm{p}}$ is the plastic volumetric strain increment.

Karagozian and Case concrete model - Release III

The K&C concrete model - Release III in LS-DYNA, also known as the K&C concrete model, or simply the KCC model, is controlled by three deviatoric stress limit surfaces, that is, the initial yield surface $\mathrm{\Delta }{\sigma }_{\mathrm{y}}$ , the maximum failure surface $\mathrm{\Delta }{\sigma }_{\mathrm{m}}$ and the residual failure surface $\mathrm{\Delta }{\sigma }_{\mathrm{r}}$ , given by

$\begin{array}{c}\mathrm{\Delta }{\sigma }_{\mathrm{y}}={\alpha }_{0\mathrm{y}}+\frac{P}{{\alpha }_{1\mathrm{y}}+{\alpha }_{2\mathrm{y}}P}\\ \mathrm{\Delta }{\sigma }_{\mathrm{m}}={\alpha }_{0\mathrm{m}}+\frac{P}{{\alpha }_{1\mathrm{m}}+{\alpha }_{2\mathrm{m}}P}\\ \mathrm{\Delta }{\sigma }_{\mathrm{r}}=\frac{P}{{\alpha }_{1\mathrm{r}}+{\alpha }_{2\mathrm{r}}P}\end{array}$

(14)

where $P$ again denotes the hydrostatic pressure and the α‐parameters are user inputs. The ${\alpha }_{0m}$ , ${\alpha }_{1m}$ and ${\alpha }_{2m}$ parameters are calibrated from uniaxial unconfined compression tests, while the ${\alpha }_{0y}$ , ${\alpha }_{1y}$ , ${\alpha }_{2y}$ , ${\alpha }_{1r}$ and ${\alpha }_{2r}$ are calibrated from triaxial compression paths. Further, a movable yield surface replaces the initial yield surface once it has been reached, with a linear interpolation from the initial yield surface to the maximum failure surface given as

$\mathrm{\Delta }\sigma =\eta \left(\mathrm{\Delta }{\sigma }_{\mathrm{m}}-\mathrm{\Delta }{\sigma }_{\mathrm{y}}\right)+\mathrm{\Delta }{\sigma }_{\mathrm{y}}$

(15)

where $\eta \in \left[\mathrm{0,1}\right]$ denotes the damage and is a user input series relying on the accumulative effective plastic strain parameter $\lambda$ (also called ‘modified’ effective plastic strain, not to be confused with the actual effective plastic strain). The damage parameter $\eta$ initiates at zero for $\lambda =0$ and increases to unity at a given value $\lambda ={\lambda }_{\mathrm{m}}$ . Afterwards, softening commences and $\eta$ decreases to zero with the increase of $\lambda$ . Therefore, the movable yield surface is interpolated from the maximum failure surface to the residual failure surface for $\lambda >{\lambda }_{\mathrm{m}}$ according to

$\mathrm{\Delta }\sigma =\eta \left(\mathrm{\Delta }{\sigma }_{\mathrm{m}}-\mathrm{\Delta }{\sigma }_{\mathrm{r}}\right)+\mathrm{\Delta }{\sigma }_{\mathrm{r}}$

(16)

The K&C concrete model also includes the Willam and Warnke (1975) reduction factor formulation $\mathcal{R}\left(\theta ,e\right)$ to describe the shape of deviatoric plane with respect to the Lode angle $\theta$ and the normalized shape factor $e$ as defined in equation (6).

The shape factor $e$ , denoted $\psi$ by Malvar et al. (1997), is again the ratio of tensile to compressive meridian radii at a given pressure $P$ , and is expressed as

$\psi \left(P\right)=\{\begin{array}{c}\frac{1}{2},P\le 0\\ \hspace{0.17em}\frac{1}{2}+\frac{3f_{\mathrm{t}}}{2f_{\mathrm{c}}},P=\frac{f_{\mathrm{c}}}{3}\\ \begin{array}{c}\frac{1.15f_{\mathrm{c}}}{{\alpha }_{0\mathrm{m}}+\frac{2\alpha \raisebox{1ex}{$f_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}{{\alpha }_{1\mathrm{m}}+2{\alpha }_{2\mathrm{m}}1.15\raisebox{1ex}{$f_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}},P=\frac{2\alpha f_{\mathrm{c}}}{3}\\ \hspace{0.17em}0.753,P=f_{\mathrm{c}}\\ \hspace{0.17em}1,P\ge 8.45f_{\mathrm{c}}\end{array}\end{array}$

(17)

The volumetric pore collapse evolution is introduced as an equation of state (EoS), which includes a set of tabulated data of the hydrostatic pressure $P$ , the volumetric strain ${\epsilon }_{\mathrm{v}}$ and the unloading bulk modulus $K_{\mathrm{U}}$ . Then, the factor $\phi$ is computed internally as a function of the EoS input data. The bulk modulus $K^{\prime }$ is scaled between the loading bulk modulus $K_{\mathrm{L}}$ and the unloading bulk modulus $K_{\mathrm{U}}$ as a function of the factor $\phi$ and is characterized as

$K^{\text{'}}=\left(K_{\mathrm{L}}-K_{\mathrm{U}}\right)\cdot \exp \left(-5.55\phi \right)+K_{\mathrm{U}}$

(18)

The scaled bulk modulus $K^{\prime }$ and Poisson’s ratio $\nu$ are used to calculate the shear modulus $G$ as

$G=\frac{\left(1.5-3\nu \right)K^{\prime }}{\left(1+\nu \right)}$

(19)

The strain-rate enhancement factor $r_{\mathrm{f}\mathrm{c}}$ for compressive loading and $r_{\mathrm{f}\mathrm{t}}$ for tensile loading are adopted from Malvar and Crawford (1998) the CEB-FIB mode code (1990), that is,

$r_{\mathrm{f}\mathrm{c}}=\{\begin{array}{c}{\left(\frac{\dot{\epsilon }}{3\cdot {10}^{-5}}\right)}^{1.026\alpha }, \left|\dot{\epsilon }\right|\ge 30 {\text{s}}^{-1}\\ \gamma {\left(\frac{\dot{\epsilon }}{3\cdot {10}^{-5}}\right)}^{1/3}, \left|\dot{\epsilon }\right|<30 {\text{s}}^{-1}\end{array}$

(20)

in which $\dot{\epsilon }$ is the strain rate, $\alpha =\frac{1}{5+9\cdot \frac{f_{\mathrm{c}}}{f_{\mathrm{c}\mathrm{o}}^{\prime }}}$ and $\gamma ={10}^{6.15\alpha -2}$ with $f_{\mathrm{c}}$ being the unconfined compressive strength of concrete and $f_{\mathrm{c}\mathrm{o}}^{\prime }=10\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ , and

$r_{\mathrm{f}\mathrm{t}}=\{\begin{array}{c}{\left(\frac{\dot{\epsilon }}{1\cdot {10}^{-6}}\right)}^{\delta }, \dot{\epsilon }\hspace{0.17em}\le 1 {\text{s}}^{-1}\\ \beta {\left(\frac{\dot{\epsilon }}{1\cdot {10}^{-6}}\right)}^{1/3}, \dot{\epsilon }>1 {\text{s}}^{-1}\end{array}$

(21)

where $\delta =\frac{1}{1+\frac{8f_{\mathrm{c}}}{f_{\mathrm{c}\mathrm{o}}^{\text{'}}}}$ and $\beta ={10}^{6\delta -2}$ .

The accumulative effective plastic strain parameter $\lambda$ is determined as a function of the effective plastic strain increment $\mathrm{d}\overline{{\epsilon }^{\mathrm{p}}}=\sqrt{\frac{2}{3}{\mathrm{d}\epsilon }_{ij}^{\mathrm{p}}{\mathrm{d}\epsilon }_{ij}^{\mathrm{p}}}$ and the strain-rate enhancement factor $r_{\mathrm{f}}$ . The two scalar parameters $b_1$ and $b_2$ allow different softening in tension and compression, that is,

$\lambda =\{\begin{array}{c}{\int }_0^{\overline{{\epsilon }^{\mathrm{p}}}}\frac{\mathrm{d}\overline{{\epsilon }^{\mathrm{p}}}}{r_{\mathrm{f}\mathrm{c}}{\left(1+\frac{P}{r_{\mathrm{f}\mathrm{c}}\hspace{0.17em}{f}_{\mathrm{t}}}\right)}^{b_1}},P\ge 0\hspace{0.17em}\left(\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\right)\\ {\int }_0^{\overline{{\epsilon }^{\mathrm{p}}}}\frac{\mathrm{d}\overline{{\epsilon }^{\mathrm{p}}}}{r_{\mathrm{f}\mathrm{t}}{\left(1+\frac{P}{r_{\mathrm{f}\mathrm{t}}\hspace{0.17em}{f}_{\mathrm{t}}}\right)}^{b_2}},P<0\hspace{0.17em}\left(\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\right)\end{array}$

(22)

When the stress path is found along the vicinity of the triaxial tensile meridian, which is measured by the ratio of the deviatoric stress $\sqrt{3J_2}=q$ to the hydrostatic pressure $P$ , a volumetric increment $\mathrm{\Delta }\lambda$ is defined as

$\mathrm{\Delta }\lambda =b_3{f}_{\mathrm{d}}{k}_{\mathrm{d}}\left({\epsilon }_{\mathrm{v}}-{\epsilon }_{\mathrm{v},\mathrm{y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}\right)$

(23)

where $b_3$ describes the softening due to triaxial tensile strain, $k_{\mathrm{d}}$ is an internal scalar multiplier, ${\epsilon }_{\mathrm{v}}$ is the volumetric strain and ${\epsilon }_{\mathrm{v},\mathrm{y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}$ is the volumetric yield strain. The factor $f_{\mathrm{d}}$ is given as

$f_{\mathrm{d}}=\{\begin{array}{c}\frac{1-|\sqrt{3J_2}/P|}{0.1},0\le \left|\frac{\sqrt{3J_2}}{P}\right|<0.1\\ 0,\left|\frac{\sqrt{3J_2}}{P}\right|\ge 0.1\end{array}$

(24)

This means that for stress triaxiality ratios above 0.1, equation (23) does not contribute to increasing $\lambda$ .

Continuous surface cap model

The CSCM provides a plasticity surface $f\left(I_1,J_2,J_3,\kappa \right)$ with a smooth intersection between the shear failure surface $F_f$ and the isotropic hardening surface (cap) $F_c$ , that is,

$f\left(I_1,J_2,J_3,\kappa \right)=J_2-{\mathcal{R}}^2{F_f}^2{F_c}^2$

(25)

Here, $\kappa$ is the hardening parameter that takes the value of the first invariant $I_1$ at the initial intersection between the two surfaces $F_f$ and $F_c$ . The shear failure surface $F_f$ defines the triaxial compressive meridian as a function of the first invariant $I_1$ and the coefficients $\alpha ,\beta ,\lambda ,\theta$ , and is given as

$F_f\left(I_1\right)=\alpha -\lambda e^{-\beta I_1}+\theta I_1$

(26)

The Rubin (1991) scaling function represents the reduction factor $\mathcal{R}$ for the CSCM. The Rubin scaling function determines the strength of the model for any stress state relative to the compressive meridian. A transition from triangular, to irregular hexagonal and ultimately to a circular shape of the deviatoric plane is allowed with the functions $Q_1$ for triaxial torsion and $Q_2$ for triaxial extension, that is,

$\begin{array}{c}{Q}_1\left(I_1\right)={\alpha }_1-{\lambda }_1{e}^{-{\beta }_1{I}_1}+{\theta }_1{I}_1\\ Q_2\left(I_1\right)={\alpha }_2-{\lambda }_2{e}^{-{\beta }_2{I}_1}+{\theta }_2{I}_1\end{array}$

(27)

The hardening cap $F_c$ represents the plastic volumetric deformation due to pore collapse under high confining pressures (compaction) with a non-dimensional function attributed to Pelessone (1989), and is defined as

$F_c\left(I_1,\kappa \right)=\{\begin{array}{c}1,I_1\le L\left(\kappa \right)\\ 1-\frac{{\left[I_1-L\left(\kappa \right)\right]}^2}{2{\left[X\left(\kappa \right)-L\left(\kappa \right)\right]}^2},I_1>L\left(\kappa \right)\end{array}$

(28)

The functions $X\left(\kappa \right)$ and $L\left(\kappa \right)$ control the geometry of the hardening surface. A coefficient ${\kappa }_0$ is defined as the value of $I_1$ at the initial intersection of the hardening and the shear surfaces before hardening, viz.

$L\left(\kappa \right)=\{\begin{array}{c}\kappa ,\kappa >{\kappa }_0\\ {\kappa }_0,\kappa \le {\kappa }_0\end{array}$

(29)

Thus, ${\kappa }_0$ is the value of $I_1$ at the initial intersection location of the shear surface $F_f$ and the cap $F_c$ . The motion of the hardening cap is controlled by the parameter $\kappa$ , which intersects the $I_1$ axis at the value $X\left(\kappa \right)$ and it follows an elliptic path determined by the ratio $R$ . The function of $X\left(\kappa \right)$ is following the relation

$X\left(\kappa \right)=L\left(\kappa \right)+R{F}_f\left[L\left(\kappa \right)\right]$

(30)

The volumetric pore collapse evolution is expressed as a hardening rule of the cap motion by the function

${\epsilon }_{\mathrm{v}}^{\mathrm{p}}=W\left[1-e^{-D_1\left[X\left(\kappa \right)-X\left({\kappa }_0\right)\right]-D_2{\left[X\left(\kappa \right)-X\left({\kappa }_0\right)\right]}^2}\right]$

(31)

where ${\epsilon }_{\mathrm{v}}^{\mathrm{p}}$ is the plastic volumetric strain, $W$ is the maximum plastic volumetric strain of the consolidated material, $D_1$ and $D_2$ are shape factors of the pressure-volumetric strain curve and $X\left({\kappa }_0\right)$ denotes the pre-consolidation pressure. The CSCM does not involve bulk modulus alternation between the elastic and solidification phases.

Further, an initial yield surface $F_y{=N_{\mathrm{H}}F}_f$ with $N_{\mathrm{H}}\le 1$ indicates the initial point of yielding for low confining pressures. The initial yielding surface hardens until it coincides with the failure surface after hardening is engaged. After reaching the failure surface the model follows a softening behaviour with brittle damage in tension and ductile damage in compression, that is,

${\sigma }_{ij}^{\mathrm{d}}=\left(1-d\right){\sigma }_{ij}^{\mathrm{v}\mathrm{p}}$

(32)

where $d$ is the damage parameter, ${\sigma }_{ij}^{\mathrm{v}\mathrm{p}}$ is the viscoplastic stress tensor without damage and ${\sigma }_{ij}^{\mathrm{d}}$ is the stress tensor with damage.

The damage accumulation due to softening is given as

$d=\{\begin{array}{c}\frac{0.999}{D}\left[\frac{1+D}{1+D{e}^{-C\left({\tau }_t-{\tau }_{0t}\right)}}-1\right],\mathrm{B}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{t}\mathrm{l}\mathrm{e}\hspace{0.17em}\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\\ \frac{d_{\max }}{B}\left[\frac{1+B}{1+B{e}^{-A\left({\tau }_c-{\tau }_{0c}\right)}}-1\right],\mathrm{D}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{e}\hspace{0.17em}\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\end{array}$

(33)

where $A,B,C$ and $D$ are softening shape factors, ${\tau }_t$ and ${\tau }_c$ are energy terms computed by the model at the current strain, while ${\tau }_{0t}$ and ${\tau }_{0c}$ are the initial damage energy thresholds calculated from the shear failure surface in tension and compression, respectively.

The fracture energy in uniaxial tension ( $\mathrm{G}\mathrm{F}\mathrm{T}$ ), uniaxial compression ( $\mathrm{G}\mathrm{F}\mathrm{C}$ ) and pure shear ( $\mathrm{G}\mathrm{F}\mathrm{S}$ ) are user input parameters. The model computes the fracture energy of the stress state with interpolation equations and the transition parameters $pwrt$ (shear-to-tension) and $pwrc$ (shear-to-compression) as

$G_f=\{\begin{array}{c}\mathrm{G}\mathrm{F}\mathrm{S}+{\left(-\frac{I_1}{\sqrt{3J_2}}\right)}^{pwrt}\left(\mathrm{G}\mathrm{F}\mathrm{T}-\mathrm{G}\mathrm{F}\mathrm{S}\right),\mathrm{T}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{e}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\\ \mathrm{G}\mathrm{F}\mathrm{S}+{\left(\frac{I_1}{\sqrt{3J_2}}\right)}^{pwrc}\left(\mathrm{G}\mathrm{F}\mathrm{C}-\mathrm{G}\mathrm{F}\mathrm{S}\right),\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{v}\mathrm{e}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\end{array}$

(34)

The strength and the fracture energy of the model increase with increasing strain rate. The CSCM inducts strain-rate sensitivity with a viscoplastic algorithm which interpolates between the trial elastic stress ${\sigma }_{ij}^{\mathrm{T}}$ and the plastic response without rate effects ${\sigma }_{ij}^{\mathrm{p}}$ . The interpolation is controlled by the parameter $\gamma$ , according to

${\sigma }_{ij}^{\mathrm{v}\mathrm{p}}=\left(1-\gamma \right){\sigma }_{ij}^{\mathrm{T}}+{\gamma \sigma }_{ij}^{\mathrm{p}}\hspace{0.17em}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\hspace{0.17em}\gamma =\frac{∆t/\eta }{1+∆t/\eta }$

(35)

The parameter $\gamma$ is dependent on the time step $∆t$ , while the fluidity coefficient $\eta$ is computed internally by a set of coefficients ${\eta }_s,{\eta }_t$ and ${\eta }_c$ in shear, tension and compression as

$\eta =\{\begin{array}{c}{\eta }_s+{\left(-\frac{I_1}{\sqrt{3J_2}}\right)}^{pwrt}\left({\eta }_t-{\eta }_s\right),\mathrm{T}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{e}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\\ {\eta }_s+{\left(\frac{I_1}{\sqrt{3J_2}}\right)}^{pwrc}\left({\eta }_c-{\eta }_s\right),\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{v}\mathrm{e}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\end{array}$

(36)

The three fluidity coefficients are calculated from five user input parameters, where ${\eta }_{0c}$ is the reference fluidity in compression, ${\eta }_{0t}$ is the reference fluidity in tension, $S_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}$ is the ratio of shear-to-tensile fluidity, and the rate effect exponents $NT$ and $NC$ for tension and compression, respectively, given as

${\eta }_t=\frac{{\eta }_{0t}}{{\dot{\epsilon }}^{NT}},{\eta }_c=\frac{{\eta }_{0c}}{{\dot{\epsilon }}^{NC}},{\eta }_s=S_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}\times {\eta }_t\hspace{0.17em}$

(37)

The overstress limits in tension ( ${over}_t$ ) and compression ( ${over}_c$ ) prevent the strain-rate effects at very high strains with

$if\hspace{0.17em}\eta >\frac{{over}_{t/c}}{E\dot{\epsilon }}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}\hspace{0.17em}\eta =\frac{{over}_{t/c}}{E\dot{\epsilon }}$

(38)

Further, the CSCM allows strain-rate sensitivity on the fracture energy as

$G_f^{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}=G_f{\left(1+\frac{E\dot{\epsilon }\eta }{f^{\prime }}\right)}^{repow}$

(39)

where $G_f^{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}$ is the enhanced fracture energy by strain-rate effects, $f^{\prime }$ is the yield strength before applying rate effects and $repow$ is the exponent of the strain-rate effects to the fracture energy.

Riedel–Hiermaier–Thoma model

The RHT model couples a deviatoric yield surface ${\sigma }_{\mathrm{y}}$ with an EoS. Both depend on the hydrostatic pressure $P$ , where the former includes the stress triaxiality factor $\mathcal{R}\left(\theta ,Q\right)$ and the strain-rate increase factor $F_{\mathrm{r}}$ , while the latter considers pore closure. The yield surface ${\sigma }_{\mathrm{y}}$ is given as

${\sigma }_{\mathrm{y}}=f_{\mathrm{c}}{\sigma }_{\mathrm{y}}^{*}\left(P^{*},F_{\mathrm{r}},{\epsilon }_{\mathrm{p}}^{*}\right)\mathcal{R}\left(\theta ,Q\right)$

(40)

It involves the unconfined compressive strength $f_{\mathrm{c}}$ , the normalized compression meridian ${\sigma }_{\mathrm{y}}^{*}={\sigma }_{\mathrm{y}}/f_{\mathrm{c}}$ and the reduction factor $\mathcal{R}\left(\theta ,Q\right)$ . Here, ${\sigma }_{\mathrm{y}}^{*}$ is expressed as a function of the normalized failure surface ${\sigma }_{\mathrm{f}}^{*}={\sigma }_{\mathrm{f}}/f_{\mathrm{c}}$ and the factor $\gamma$ . The factor $\gamma$ is the product of the elastic limit surface $F_{\mathrm{e}}$ and the cap function $F_{\mathrm{c}\mathrm{a}\mathrm{p}}$ and is controlled by the parameter ${\epsilon }_{\mathrm{p}}^{*}$ . Once the initial elastic limit surface $F_{\mathrm{e}}$ is reached the plastic strain ${\epsilon }_{\mathrm{p}}>0$ is initiated. Then, the yield surface ${\sigma }_{\mathrm{y}}$ is using the parameter ${\epsilon }_{\mathrm{p}}^{*}$ to interpolate between the initial yield surface $\left({\epsilon }_{\mathrm{p}}^{*}=0\right)$ and the failure surface $\left({\epsilon }_{\mathrm{p}}^{*}=1\right)$ , such that

$\begin{array}{c}{\sigma }_{\mathrm{y}}^{*}={\sigma }_{\mathrm{f}}^{*}\left(\frac{P^{*}}{\gamma },F_{\mathrm{r}}\right)\gamma \\ \gamma =\left[{\epsilon }_{\mathrm{p}}^{*}+\left(1-{\epsilon }_{\mathrm{p}}^{*}\right)F_{\mathrm{e}}{F}_{\mathrm{c}\mathrm{a}\mathrm{p}}\right]\\ {\epsilon }_{\mathrm{p}}^{*}=\min \left(\frac{{\epsilon }_{\mathrm{p}}}{{\epsilon }_{\mathrm{p}}^{\mathrm{h}}},1\right),{\epsilon }_{\mathrm{p}}^{\mathrm{h}}=\frac{{\sigma }_{\mathrm{y}}\left(1-F_{\mathrm{e}}{F}_{\mathrm{c}\mathrm{a}\mathrm{p}}\right)}{3\gamma G}\end{array}$

(41)

The RHT model also employs the Willam and Warnke (1975) formulation to determine the reduction factor $\mathcal{R}\left(\theta ,Q\right)$ as given in equation (6), where $\theta$ is the Lode angle and $Q$ is the normalized shape factor (denoted $e$ in the MHJC model and $\psi$ in the K&C concrete model) of the deviatoric plane. Like in the MHJC and K&C concrete models, $Q$ depends on the normalized hydrostatic pressure $P^{*}$ , where $Q_0$ and $B$ are shape factor coefficients, that is,

$Q=Q_0+B{P}^{*}$

(42)

The normalized failure surface ${\sigma }_{\mathrm{f}}^{*}$ is determined by equation (43), in which $A$ and $n$ are failure surface parameters, $f_{\mathrm{t}}^{*}$ and $f_{\mathrm{s}}^{*}$ are the tensile and shear strengths relative to the compressive strength $f_{\mathrm{c}}$ , $P_{\mathrm{t}}^{*}=\frac{F_{\mathrm{r}}\hspace{0.17em}Q{f}_{\mathrm{s}}^{*}{f}_{\mathrm{t}}^{*}}{3\left(Q_1{f}_{\mathrm{t}}^{*}-Q{f}_{\mathrm{s}}^{*}\right)}$ is the tensile cut-off pressure, $Q_1=\mathcal{R}\left(\frac{\pi }{6},0\right)$ is the value of the stress reduction factor of the torsion meridian at $P^{*}=0$ .

${\sigma }_{\mathrm{f}}^{*}=\{\begin{array}{c}A{\left[P^{*}-\frac{F_{\mathrm{r}}}{3}+{\left(\frac{A}{F_{\mathrm{r}}}\right)}^{-1/n}\right]}^n,{3P}^{*}\ge F_{\mathrm{r}}\\ \frac{F_{\mathrm{r}}\hspace{0.17em}{f}_{\mathrm{s}}^{*}}{Q_1}+{3P}^{*}\left(1-\frac{f_{\mathrm{s}}^{*}}{Q_1}\right),F_{\mathrm{r}}>{3P}^{*}\ge 0\\ \frac{F_{\mathrm{r}}\hspace{0.17em}{f}_{\mathrm{s}}^{*}}{Q_1}-{3P}^{*}\left(\frac{1}{Q}-\frac{f_{\mathrm{s}}^{*}}{Q_1{f}_{\mathrm{t}}^{*}}\right),0>{3P}^{*}\ge {3P}_{\mathrm{t}}^{*}\\ 0,{3P}_{\mathrm{t}}^{*}>{3P}^{*}\end{array}$

(43)

The elastic limit $F_{\mathrm{e}}$ ranges between the ratios of the elastic compressive stress ${\mathrm{g}}_{\mathrm{c}}^{*}$ and the tensile stress ${\mathrm{g}}_{\mathrm{t}}^{*}$ over the respective ultimate strength, given as

$F_{\mathrm{e}}=\{\begin{array}{c}\hspace{0.17em}{\mathrm{g}}_{\mathrm{c}}^{*},{3P}^{*}\ge F_{\mathrm{r}}^{\mathrm{c}}{\mathrm{g}}_{\mathrm{c}}^{*}\\ {\mathrm{g}}_{\mathrm{c}}^{*}-\frac{{3P}^{*}-F_{\mathrm{r}}^{\mathrm{c}}{\mathrm{g}}_{\mathrm{c}}^{*}}{F_{\mathrm{r}}^{\mathrm{c}}{\mathrm{g}}_{\mathrm{c}}^{*}+F_{\mathrm{r}}^{\mathrm{t}}{\mathrm{g}}_{\mathrm{t}}^{*}{f}_{\mathrm{t}}^{*}}\left({\mathrm{g}}_{\mathrm{t}}^{*}-{\mathrm{g}}_{\mathrm{c}}^{*}\right),{F_{\mathrm{r}}^{\mathrm{c}}{\mathrm{g}}_{\mathrm{c}}^{*}>3P}^{*}\ge -F_{\mathrm{r}}^{\mathrm{t}}{\mathrm{g}}_{\mathrm{t}}^{*}{f}_{\mathrm{t}}^{*}\\ \hspace{0.17em}{\mathrm{g}}_{\mathrm{t}}^{*},-F_{\mathrm{r}}^{\mathrm{t}}\hspace{0.17em}{\mathrm{g}}_{\mathrm{t}}^{*}{f}_{\mathrm{t}}^{*}\ge {3P}^{*}\end{array}$

(44)

The cap function $F_{\mathrm{c}\mathrm{a}\mathrm{p}}$ intersects the elastic limit surface $F_{\mathrm{e}}$ at the lower cap pressure $P_{\mathrm{u}}^{*}$ . The former follows a parabolic path until it meets the pressure axis at the pore crush pressure $P_{\mathrm{c}}^{*}\left(a\right)$ . The cap function is represented by

$F_{\mathrm{c}\mathrm{a}\mathrm{p}}=\{\begin{array}{c}\hspace{0.17em}0,P^{*}\ge P_{\mathrm{c}}^{*}\left(a\right)\\ \sqrt{1-{\left(\frac{P^{*}-P_{\mathrm{u}}^{*}}{P_{\mathrm{c}}^{*}\left(a\right)-P_{\mathrm{u}}^{*}}\right)}^2},{P_{\mathrm{c}}^{*}\left(a\right)>P}^{*}\ge P_{\mathrm{u}}^{*}\\ \hspace{0.17em}1,P_{\mathrm{u}}^{*}\ge P^{*}\end{array}$

(45)

The pore crush pressure $P_{\mathrm{c}}\left(a\right)$ is represented by a compaction model depending on a porosity history variable $a\left(t\right)$ . The compaction model is described by the initial pore crush pressure $P_{\mathrm{e}\mathrm{l}}$ , the compaction pressure $P_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}$ , the hydrostatic pressure $P$ , the porosity exponent $N$ and the initial porosity ${\alpha }_0$ , so that

$\begin{array}{c}{P}_c\left(a\right)=P_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}-\left(P_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}-P_{\mathrm{e}\mathrm{l}}\right){\left(\frac{\alpha \left(t\right)-1}{{\alpha }_0-1}\right)}^{1/N}\\ a\left(t\right)=1+\left({\alpha }_0-1\right){\left(\frac{P_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}-P}{P_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}-P_{\mathrm{e}\mathrm{l}}}\right)}^N\end{array}$

(46)

Pore closure is characterized by a Hugoniot polynomial EoS, where $\eta =\frac{a\rho }{{\alpha }_0{\rho }_0}-1$ is the plastic volumetric strain, ${\rho }_0$ is the initial density and $e$ is the internal energy, ${\mathrm{{\rm B}}}_0$ , ${\mathrm{{\rm B}}}_1$ , ${\mathrm{{\rm A}}}_2$ , ${\mathrm{{\rm A}}}_3$ and $T_2$ are the parameters of the polynomial EoS, while ${\mathrm{{\rm A}}}_1$ and ${\mathrm{{\rm T}}}_1$ are the bulk modulus, in compression and tension, respectively, viz.

$P=\frac{1}{\alpha }\{\begin{array}{c}\hspace{0.17em}\\ \left({\mathrm{{\rm B}}}_0+{\mathrm{{\rm B}}}_1\eta \right)\alpha {\rho }_{0}e+{\mathrm{{\rm A}}}_1\eta +{\mathrm{{\rm A}}}_2{\eta }^2+{\mathrm{{\rm A}}}_3{\eta }^3,\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{v}\mathrm{e}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\\ \hspace{0.17em}{\hspace{0.17em}\mathrm{{\rm B}}}_0\alpha \rho e+{\mathrm{{\rm T}}}_1\eta +{\mathrm{{\rm T}}}_2{\eta }^2,\mathrm{T}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{e}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\end{array}$

(47)

Alternatively, in case of ${\mathrm{{\rm B}}}_0=0$ , the EoS involves a Grüneisen parameter $\mathrm{\Gamma }$ , whereas equations (47) and (48) are identical for ${\mathrm{{\rm B}}}_0={\mathrm{{\rm B}}}_1=\mathrm{\Gamma }=0$ . Thus,

$P=\mathrm{\Gamma }{\rho }_{0}e+\left({\mathrm{{\rm A}}}_1\eta +{\mathrm{{\rm A}}}_2{\eta }^2+{\mathrm{{\rm A}}}_3{\eta }^3\right)\left(1-\frac{1}{2}\mathrm{\Gamma }\eta \right)$

(48)

The strain-rate dependence increase factor is calculated by

$F_{\mathrm{r}}=\{\begin{array}{c}\hspace{0.17em}{F}_{\mathrm{r}}^{\mathrm{c}},{3P}^{*}\ge F_{\mathrm{r}}^{\mathrm{c}}\\ F_{\mathrm{r}}^{\mathrm{c}}-\frac{{3P}^{*}-F_{\mathrm{r}}^{\mathrm{c}}}{F_{\mathrm{r}}^{\mathrm{c}}+F_{\mathrm{r}}^{\mathrm{t}}{f}_{\mathrm{t}}^{*}}\left(F_{\mathrm{r}}^{\mathrm{t}}-F_{\mathrm{r}}^{\mathrm{c}}\right),{F_{\mathrm{r}}^{\mathrm{c}}>3P}^{*}\ge -F_{\mathrm{r}}^{\mathrm{t}}{f}_{\mathrm{t}}^{*}\\ {\hspace{0.17em}F}_{\mathrm{r}}^{\mathrm{t}},-F_{\mathrm{r}}^{\mathrm{t}}{f}_{\mathrm{t}}^{*}\ge {3P}^{*}\end{array}$

(49)

in which $F_{\mathrm{r}}^{\mathrm{c}/\mathrm{t}}$ is computed by a power law, given as

$F_{\mathrm{r}}^{\mathrm{c}/\mathrm{t}}=\{\begin{array}{c}\hspace{0.17em}\\ {\left(\frac{{\dot{\epsilon }}_{\mathrm{p}}}{{{\dot{\epsilon }}_0}^{\mathrm{c}/\mathrm{t}}}\right)}^{{\beta }^{\mathrm{c}/\mathrm{t}}},{\dot{\epsilon }}_{\mathrm{p}}\le {{\dot{\epsilon }}_{\mathrm{p}}}^{\mathrm{c}/\mathrm{t}}\\ {\hspace{0.17em}\gamma }^{\mathrm{c}/\mathrm{t}}\hspace{0.17em}\sqrt[3]{{\dot{\epsilon }}_{\mathrm{p}}},{\dot{\epsilon }}_{\mathrm{p}}>{{\dot{\epsilon }}_{\mathrm{p}}}^{\mathrm{c}/\mathrm{t}}\end{array}$

(50)

This power law is nearly identical to the CEB-FIB mode code (1990). Superscripts ‘c’ and ‘t’ denote compression and tension, ${\dot{\epsilon }}_{\mathrm{p}}$ is the plastic strain rate, ${{\dot{\epsilon }}_0}^{\mathrm{c}/\mathrm{t}}$ is a reference strain rate and ${{\dot{\epsilon }}_{\mathrm{p}}}^{\mathrm{c}/\mathrm{t}}$ is the break strain rate, while ${\beta }^{\mathrm{c}/\mathrm{t}}$ is the strain-rate dependence exponent and ${\gamma }^{\mathrm{c}/\mathrm{t}}$ is the internally calculated coefficient for the continuity of equation (50) at ${{\dot{\epsilon }}_{\mathrm{p}}}^{\mathrm{c}/\mathrm{t}}$ .

Material and ballistic impact tests

To evaluate the chosen concrete models and applied calibration methods, the experimental data presented in Kristoffersen et al. (2021b) will be utilized. In that study, the ballistic perforation resistance of 50 mm thick concrete slabs impacted by 20 mm diameter ogive-nose steel projectiles was investigated both experimentally and numerically. Three types of commercially produced concrete with nominal unconfined compressive strengths of 35, 75 and 110 MPa were used to cast material test specimens and slabs. After curing, ballistic impact tests were carried out in a compressed gas gun facility to determine the ballistic limit curve and velocity for each slab quality. In this study we will only consider data from the C75 concrete tests, and some of the main results are briefly described below for completeness. For a full description of the various material and ballistic impact tests, the reader is referred to Kristoffersen et al. (2021b).

Material tests

The composition of the C75 concrete is given in Table 1. Three common quasi-static concrete material tests – cylinder compression, cube compression and tensile splitting – were carried out. Cylinder compression and tensile splitting are used herein in the calibration of the four concrete material models presented. Ideally, triaxial test and dynamic test data should have been available for a complete calibration of the various models. However, in this study, we investigate whether the accuracy of the models is improved given only the standard material tests and some data from the literature for the calibration of their parameters.

Table 1.

Composition and mechanical properties of the C75 concrete.

Composition		Mechanical properties
Cement (kg/m3)	427	Cylinder compressive strength (MPa)	72.8
Silica fume (kg/m3)	48	Cube compressive strength (MPa)	87.4
Water (kg/m3)	168	Tensile splitting strength (MPa)	5.2
Fine aggregate: 0–8 mm (kg/m3)	917	Density (kg/m3)	2506
Coarse aggregate: 8–16 mm (kg/m3)	881	Young’s modulus $E$ (GPa)	42
W/C ratio	0.39	Poisson’s ratio $\nu$ .	0.2

Simultaneously as the concrete slabs for the ballistic impact tests were cast in custom-made wooden moulds, many cubes and cylinders were cast in reusable steel moulds and kept submerged in water at room temperature until the day of testing. The cylinders were 100 mm in diameter and 200 mm long, while the cubes had side lengths of 100 mm. All tests were carried out in a fully automated Toni Tech 3000 kN load controlled rig at a loading rate of 0.08 MPa/s. During testing, the force was measured and synchronized with pictures from a high-resolution digital camera for two-dimensional digital image correlation (2D-DIC). Table 1 gives average values found in Kristoffersen et al. (2021b), which are based on three tests of the cylinder compressive strength $f_{\mathrm{c}}$ , the cube compressive strength $f_{\mathrm{c}\mathrm{c}}$ , the tensile splitting strength $f_{\mathrm{t}}$ and the mass density ${\rho }_0$ of the C75 concrete 28 days after casting. Note that the elastic constants $E$ and $\nu$ were recalibrated for the present study.

The results from the force and 2D-DIC measurements in terms of average engineering stress-strain curves from cylinder compression, cube compression and tensile splitting tests of the C75 concrete are given in Figure 3. The engineering stress-strain curves from the cylinder compression tests and the tensile splitting tests were used for the calibration of the material constants in the various models based on inverse modelling in numerical simulations.

Figure 3.

Average engineering stress-engineering strain curves from cylinder compression, cube compression, and tensile splitting tests of the C75 concrete (Kristoffersen et al., 2021b).

Ballistic impact tests

The ballistic impact tests in Kristoffersen et al. (2021b) were performed in a compressed gas gun facility described by Børvik et al. (1999). During testing, the ogive-nose projectiles were fired at impact velocities just below and well above the ballistic limit of the concrete slabs by using compressed air. The projectiles were manufactured from tool steel and heat treated to a Rockwell C value of 53 after machining. Nominal diameter ( $d_{\mathrm{p}}$ = 20 mm), length ( $l_{\mathrm{p}}$ = 98 mm), mass ( $m_{\mathrm{p}}$ = 196 g) and critical-radius-head (CRH = 3) of the projectile were kept constant in all tests, and the exact geometry is shown Kristoffersen et al. (2021b). The plain (i.e., no rebar or steel fibres) concrete slabs with dimensions 625 mm × 625 mm × 50 mm were securely fixed to a rigid boundary in the impact camber of the gas gun by massive clamps at each corner.

Accurate optical measurements of the initial velocity $v_{\mathrm{i}}$ and the residual velocity $v_{\mathrm{r}}$ of the projectile were provided by two synchronised Phantom v2511 high-speed cameras operating at 50 000 fps, while lighting was offered by two flashlights in combination with powerful halogen lamps. The high-speed cameras were also used to measure projectile pitch before impact, and to capture high-resolution videos of the penetration and perforation process.

Based on the measured initial and residual projectile velocities, the ballistic limit curve and velocity of the concrete slabs were estimated. This was done by a least squares fit of the model constants in the Recht-Ipson model (Recht and Ipson, 1963) to the measured data. The model is defined as

$v_{\mathrm{r}}={a\left(v_{\mathrm{i}}^p-v_{\mathrm{b}\mathrm{l}}^p\right)}^{1/p}$

(51)

where $v_{\mathrm{i}}$ and $v_{\mathrm{r}}$ are the initial and residual velocities of the projectile, respectively, $v_{\mathrm{b}\mathrm{l}}$ is the ballistic limit velocity, and $a$ and $p$ were taken as empirical constants. Based on the least squares fit it was found that $v_{\mathrm{b}\mathrm{l}}$ = 140.4 m/s, $a$ = 0.98 and $p$ = 1.47 for the C75 concrete slab. Figure 4 shows a plot of the measured initial versus residual velocity data for this concrete slab, together with the ballistic limit curve. Figure 5(a) shows a time lapse of the penetration and perforation process from a typical test with an impact velocity close to the ballistic impact velocity. As seen, a lot of debris was ejected from the concrete slabs during impact, and it was occasionally hard to spot the projectile before it had travelled a certain distance after perforation. Here, the projectile is marked with a white rectangle in the last image. Figure 5(b) shows pictures of the rather irregular spall and scab craters from a typical test, while Figure 5(c) shows pictures of opposing cross-sectional slices of a plate showing the penetration channel and the craters. Note that the spall craters were in general considerably smaller than the scab craters and that a clear tunnelling section connecting the two craters was in general lacking in these tests since the plates were rather thin.

Figure 4.

Ballistic limit curve for the C75 concrete slab.

Figure 5.

(a) Time lapse from high-speed camera images of a typical ballistic impact tests with $v_{\mathrm{i}}$ = 157.1 m/s and $v_{\mathrm{r}}$ = 46.2 m/s, (b) crater patterns (measures in mm) after perforation and (c) opposing cross-sectional slices of a plate showing the penetration channel (outlined in yellow).

Numerical simulations

The accuracy of four common concrete models is evaluated for predicting ballistic impacts on concrete structures. Three standard models (K&C concrete model - RIII, CSCM and RTH) in the explicit FE code LS-DYNA (LST, 2020) include automatic parameter generators based solely on the unconfined compression strength of the concrete. The number of parameters required for these models set significant challenges for their calibration. In addition to the three aforementioned models, the MHJC model is selected because it comprises a comparatively small set of constitutive parameters. Therefore, the calibration of the MHJC is convenient due to its simple constitutive equations.

The constitutive parameters of the MHJC model are identified using material experiments (uniaxial compression, tensile splitting) and data from the literature (triaxial hydrostatic compression and Hopkinson bar tests). The accuracy of the concrete model is validated with the simulations of the perforation experiments. Additionally, the response of the K&C - RIII, CSCM and RHT concrete models with the default generated parameters is evaluated. Hence, the quasi-static and ballistic responses of the three standard LS-DYNA models are compared with the MHJC in numerical simulations. Finally, a study is conducted to readjust their constitutive parameters and re-evaluate their ballistic resistance with perforation simulations. The numerical results are compared with experimental data obtained from ballistic impact test at SIMLab (Kristoffersen et al., 2021b).

A 3D FE model mesh was used for the quasi-static compression (QSC) and tension (QST) tests. A rectangular prism with dimensions 0.1 m × 0.1 m × 0.2 m was created for the concrete specimen. A rectangular prismatic geometry offers a homogeneous mesh with the same size of FEs through the entire volume, which allow us to minimize mesh size dependency of the constitutive parameters. However, the experimental specimens are cylindrical; thus it is necessary to maintain slenderness ratio 2:1 for the specimen in Figure 6 (left). Shape effects have been shown to be negligible for the quasi-static case by Li et al. (2018).

Figure 6.

Finite element setups: 3D model for quasi-static compression test (left) and 2D axisymmetric model for ballistic impact tests.

A rigid material model (*MAT_020 in LS-DYNA) was selected for the steel plates. The following material parameters for the rigid steel plates were used: Young’s modulus $E_{\mathrm{p}\mathrm{l}}=208\hspace{0.17em}\mathrm{G}\mathrm{P}\mathrm{a}$ , Poisson ratio ${\nu }_{\mathrm{p}\mathrm{l}}=0.33$ and density ${\rho }_{\mathrm{p}\mathrm{l}}=7802 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ .

Contacts between the concrete mesh and the rigid plates were imposed by the penalty method with the surface-to-surface option in LS-DYNA. Baltay and Gjelsvik (1990) suggest an average friction coefficient of 0.47 between concrete and steel under quasi-static conditions, while Burak and Tuncan (2014) consider a value of 0.55. Simulations with both friction coefficients produce indistinguishable variations; thus, a value of 0.5 was selected. Hexahedral elements with side lengths of 1 mm were chosen for the concrete mesh, which gives 100 × 100 × 200 elements. The same elements were used for the rigid plates with dimensions 2 mm × 2 mm × 2.5 mm.

The boundary conditions of the QST test were imposed by fixing the bottom surface of the concrete specimen. A velocity curve was imposed on the top surface of the specimen with axial direction along its length, which creates tensile stress. The boundary conditions of the QSC experiment were modelled using two rigid steel plates. The bottom plate is fixed using translational and rotational constrains in all directions, while the top plate is free to translate in the axial direction. A velocity curve describes the motion of the top plate which applies compressive stress to the FE concrete specimen. A time scaling factor of 0.5 × 10⁻⁴ was applied, which means that the simulated time was 10–20 ms depending on the test. The strain-rate sensitivity parameters were set to zero to prevent strain-rate effects for the calibration of the parameters under quasi-static conditions.

Similarly, a 3D FE model was generated for the triaxial hydrostatic (TXH) compression tests. A cube of side 100 mm was used for the concrete specimen. Rigid steel plates were assigned on each surface of the cubic sample to assimilate the hydrostatic pressure. The rigid material model of the QSC model was adopted for the steel plates. The plates were fixed in all degrees of freedom, except the translational direction perpendicular to the contact surface. A velocity profile was applied to the unconstrained degree of freedom of each plate. The velocity direction was towards the specimen. Contact and friction coefficients between the concrete and steel were imposed identically as to the QSC model. Hexahedral elements were used with side lengths of 1 mm, which gives 100³ elements for the concrete mesh. The FE size for the rigid plates was the same as for the QSC model, that is, 2 mm × 2 mm × 2.5 mm.

A 2D axisymmetric FE model was used for the ballistic simulation with element type 15 in LS-DYNA and a size of 1.0 mm, see Figure 6 (right). Hourglass type 6 was chosen with default values to prevent zero-energy deformation modes. The elastic-plastic model with linear isotropic hardening (*MAT_003 in LS-DYNA) was selected for the projectile. The steel material parameters of the projectile were taken as Young’s modulus $E_{\mathrm{p}}=200\hspace{0.17em}\mathrm{G}\mathrm{P}\mathrm{a}$ , Poisson ratio ${\nu }_{\mathrm{p}}=0.3$ , density ${\rho }_{\mathrm{p}}=7802 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ , yield stress ${\sigma }_{\mathrm{y}}=1720 \mathrm{M}\mathrm{P}\mathrm{a}$ and tangent modulus $E_{\mathrm{t}}=15 \mathrm{G}\mathrm{P}\mathrm{a}$ . Contact was created between the concrete slab and the steel projectile by employing the penalty-based 2D automatic single surface contact option with zero friction. The concrete slab was fixed at the boundary 260 mm from the centre. For all perforation simulations, an element erosion criterion based on effective strain was added through *MAT_ADD_EROSION. This was included to remove only the most distorted elements, so for this reason the effective strain value at which the element is eroded was set deliberately high to 1.0. Using a too low value results in element erosion too early and gives inaccurate results. A small sensitivity study was run to verify that 1.0 was sufficiently high.

Modified Holmquist–Johnson–cook parameter calibration

A simplified parameter calibration procedure is proposed for the calibration of the MHJC constitutive parameters. The MHJC parameters were divided into categories and listed in calibration order: (1) the elastic behaviour, (2) the quasi-static strength and yield surface, (3) the pore collapse evolution and (4) the strain-rate sensitivity (see Table 3). The shear modulus $G=E/\left[2\left(1+\nu \right)\right]=17.5$ GPa and the initial bulk modulus $K_0=E/\left[3\left(1-2\nu \right)\right]=23.33$ GPa were calculated directly from the concrete’s mechanical properties in Table 1.

The uniaxial compressive strength parameter $f_{\mathrm{c}}=72.8\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ was introduced from the nominal value of the unconfined cylindrical compression. On the other hand, the uniaxial tensile strength parameter $f_{\mathrm{t}}=3.2\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ was identified by inverse modelling the tensile splitting test to obtain peak stress of 5.2 MPa (Figure 7 right). Quasi-static and triaxial experimental data are required to calibrate the remaining parameters of the yield surface. Still, the coefficients $B_{\mathrm{c}}=1.62$ and $N_{\mathrm{c}}=0.5$ of the yield surface were only calibrated to satisfy the peak stress of the QSC test in Figure 7 (left). The Lode-dependent reduction factor $\mathcal{R}$ is computed internally. Figure 8 (left) presents the yield surface of the compression meridian and Figure 8 (right) the shape of the deviatoric plane from low to high hydrostatic pressure. Figure 9 plots the reduction factor of the tension and shear meridians.

Figure 7.

Engineering stress-strain curves of quasi-static compression (left) and tension (right).

Figure 8.

Normalized compression meridian yield surface-to-compression strength (left) and deviatoric plane (right).

Figure 9.

Reduction factor of tension meridian (left) and shear meridian (right).

A triaxial hydrostatic experiment is essential for the determination of the pore collapse evolution parameters. Triaxial testing of concrete requires apparatus with high-pressure capacity such as a GIGA press described in Malecot et al. (2010). Consequently, the available experimental data are limited. For this reason, we estimated the pore collapse evolution parameters from data found in Malecot et al. (2010) and Magellanes et al. (2010). The pore collapse evolution parameters are the crush pressure $P_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}=\frac{f_{\mathrm{c}}}{3}\approx 24.3\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ , the solidification pressure $P_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}=500\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ , the solidification bulk modulus $K_1=1.5\cdot K_0=35\hspace{0.17em}\mathrm{G}\mathrm{P}\mathrm{a}$ , the crush volumetric strain ${\mu }_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}=P_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}/K_0=1.04\cdot {10}^{-3}$ and the fully compacted volumetric strain ${\mu }_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}=0.057714$ . Further, the solidification coefficients $K_2$ and $K_3$ were set equal to zero for simplicity, thus reducing the cubic polynomial of the solidification phase of equation (8) to a linear function. The hydrostatic test simulation results are given in Figure 10 (left), in which the slope of the first branch is $K_0$ and the slope of the third branch is $K_1$ .

Figure 10.

Triaxial hydrostatic simulation (left) ballistic limit curves (right).

The dynamic mechanical properties of concrete can be obtained by, for example, split-Hopkinson bar tests at different strain rates. The absence of dynamic experimental data for the C75 concrete compels us to use stain-rate sensitivity parameters from the literature. The parameters ${\dot{\epsilon }}_0=1.0\cdot {10}^{-5} {\mathrm{s}}^{-1}$ and $C=0.04$ were adopted from the study by Polanco-Loria et al. (2008). Moreover, we ascertained the accuracy of the strain-rate sensitivity parameters with the experimental data from the CEB-FIP model (1990) and Malvar and Crawford (1998).

Table 2 displays the material card of the MHJC model. Keep in mind that this model was implemented as a user defined material model and is thus not one of the standard models in LS-DYNA. Figure 10 (right) shows the impact velocity versus the predicted residual velocity from numerical simulations of the ballistic impact tests using the numerical mesh shown in Figure 6 (right) and the MHJC model. As seen, the proposed calibration of the MHJC model accurately predicts the residual velocities in the ballistic impact experiments.

Table 2.

Material card of the MHJC model.

*MAT_USER_DEFINED_MATERIAL_MODELS
$#	mid	rho	mt	lmc	nhv	iortho	ibluk	ig
	1	2.47E-9	42	24	25	0	8	1
$#	ivect	ifail	itherm	ihyper	ieos	lmca
	1	1	0	0	0	0
$#	G	A	Bc	C	Nc	fc	ft	K0
	17.5E3	1.0	1.62	0.04	0.5	72.8	3.2	23.33E3
$#	eps0dot	efmin	Smax	Pcrush	mucrush	Plock	mulock
	1.0E-5	0.01	7.0	24.2667	1.04E-3	500.0	0.057714
$#	alpha	beta	K1	K2	K3	fs
	0.41	0.30	35.0E3	0.0	0.0	1.0

Automatic parameter generation option – evaluation

The default parameter generation option determines the constitutive parameters of the concrete models K&C - RIII, CSCM and RHT solely based on the cylinder compressive strength of the concrete. Their elastic response, compressive and tensile strengths are illustrated in Figure 7. The CSCM presents a lower while the RHT and K&C - RIII models give a higher Young’s modulus than the experimental value. All the models capture the peak stress well in compression, although only the K&C - RIII model gives accurate peak stress in tension. The CSCM and RHT model underestimate the tensile strength.

Comparing the yield surface of all the three models to the one from the MHJC (Figure 8 left), we observe that at pressures under the uniaxial compression strength, the yield limit of the three models differs from the MHJC model. The yield surface of the K&C - RIII model remains below at all pressures. However, at higher pressures, the yield surface of the CSCM raises significantly, and the RHT model is slightly increasing. There is no coherence between the deviatoric planes of any of the models (Figure 8 right). Additionally, the reduction factors of the tension and shear meridians appear to be comparatively small for the models RHT and CSCM. The reduction factors of the K&C - RIII and MHJC models are relatively close (Figure 9). It is worth mentioning that the reduction factor of the K&C - RIII model is computed internally. In contrast, constitutive equations control the reduction factors for the RHT and CSCM models.

The hydrostatic response in Figure 10 (left) of the CSCM, K&C - RIII and MHJC models are sufficiently similar until the solidification pressure limit $(P_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}}=500\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ ) of the MHJC model. After that point, the K&C - RIII model retains a smaller bulk modulus until reaching approximately $700\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ , where the solidification phase seems to initiate. At the end, the slope of the K&C - RIII model becomes the same as for the MHJC model. The solidification phase does not arrive for the CSCM, at least until a pressure of $800\hspace{0.17em}\mathrm{M}\mathrm{P}\mathrm{a}$ . However, the slope of the compaction phase seems to be increasing, leading to the third phase at a higher pressure. The RHT model overestimates the compaction phase of the pore evolution under hydrostatic pressure, without any indication of the solidification phase for the pressures investigated herein.

The broad diversity of the performance of the models in the material simulations is reflected in the different ballistic responses as shown in Figure 10 (right). The performance of the RHT model exhibits a reduced ballistic response. Despite the accuracy of the CSCM for impact velocities beyond 250 m/s, it overrates the ballistic limit for the lower impact velocities. The K&C - RIII model reveals a tendency to reliably predict the ballistic response with impact velocities higher than 160 m/s. However, low-velocity simulations undervalue the ballistic capacity.

CSCM and RHT parameter readjustment

The investigation of the response of the LS-DYNA concrete models with the default parameters permits comprehending the necessary parameter readjustments. The K&C concrete model - RIII presents adequate performance throughout all the material simulations (uniaxial compression, uniaxial tension and triaxial hydrostatic compression). Thereby the deviation of the residual velocity at impact velocities below 150 m/s is attributed to the damage constitutive formulae of the K&C concrete model. The assessment of the damage functions is excluded in the current study. Consequently, readjustments of the default parameters for the K&C concrete model - RIII are out of scope.

Readjustments for the parameters of the CSCM and RHT models are necessary because the default parameters: (1) produce weak tensile strength, (2) fail to deliver the solidification phase of the porosity evolution under hydrostatic pressure and (3) provide a considerably low reduction factor for the tension and shear meridians. Therefore, an identical procedure to the MHJC model calibration was followed for the redetermination of the CSCM and RHT parameters.

The parameters of CSCM and RHT models that require readjustment are listed in categories in Table 3 with their calibration order. The recalibration of the following parameters employed the same simplified parameter calibration procedure described for the calibration of the MHJC model.

Table 3.

Categories of MHJC, CSCM and RHT concrete constitutive parameters.

Category	Description	MHJC	CSCM	RHT
Elasticity	Bulk modulus	K0	k	a1(compression), t1(tension)
	Shear modulus	G	g	shear
Compression meridian	Surface parameters	Bc, Nc	alpha, theta, lamda, beta	a, n
	Strength parameters	fc, ft	( - )	fc, fs*, ft*
Shear meridian &Tension meridian	Shape factor coefficients	e, ${\mathrm{P}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{*}$ (Internal constants)	alphai a , thetai a , lamdai a , betai a	qo, b
	Reduction factor	$\mathcal{R}\left(\mathrm{\theta },\mathrm{e}\right)$ (Internally)	$\mathcal{R}$ (Internally)	$\mathcal{R}\left(\mathrm{\theta },\mathrm{e}\right)$ (Internally)
Pore collapse evolution	Pressure limits	Pcrush, Plock	Xd	pel, pco
	Volumetric strain limits	mucrush, mulock	w, d1, d2	np, alpha
	Solidification modulus	K1, K2, K3	( - )	a2, a3
Strain-rate sensitivity	Reference strain rate	eps0dot	eta0c, eta0t	e0c, e0t
	Strain-rate coefficients	C	nc, nt	betac, betat

^aThe index i of the CSCM shape factor coefficients takes the value one for the shear and two for the tension meridians.

The elastic parameters were taken directly from the concrete’s mechanical properties. The compressive strength $f_{\mathrm{c}}$ and surface and parameters of the RHT and MHJC models are identical. The corresponding parameters for CSCM were identified as $\lambda =f_{\mathrm{c}}$ and the other three parameters were computed with a curve adjustment on the MHJC yield surface (Figure 12 left). The RHT strength parameters were found as $f_{\mathrm{s}}^{*}={\sigma }_{\mathrm{e}\mathrm{q}}^{*}\left(P^{*}=0\right)$ from equation (5) and $f_{\mathrm{t}}^{*}$ by inverse modelling the tensile splitting test. The newly defined parameters improve Young’s modulus and the tensile strength of the CSCM and RHT models (see Figure 11). Further, the yield surface of the compression meridian of the RHT model becomes identical to the MHJC model. The constitutive equation of the CSCM yield surface defines a slightly elevated compression meridian, though the cap function of the CSCM creates an elliptic path that compensates for the difference (Figure 12 left).

Figure 11.

Stress-strain curves of quasi-static compression (left) and tension (right).

Figure 12.

Normalized compression meridian yield surface-to-compression strength (left) and deviatoric plane (right).

The CSCM shape factor parameters were defined from equation (6); thus ${\alpha }_{\mathrm{1,2}}-{\lambda }_{\mathrm{1,2}}=\mathcal{R}\left[{\theta }_{\mathrm{1,2}},e\left(P^{*}=0\right)\right]$ and ${\beta }_i$ , ${\theta }_i$ for $P^{*}>0$ , while for the RHT model they were calculated from equation (7) as $Q_0=e\left(P^{*}=0\right)$ and $B=e\left(P^{*}=1\right)-e\left(P^{*}=0\right)$ . Figure 13 confirms that the reduction factors of the tension and shear meridians are equal for all three models CSMC, RHT and MHJC. Accordingly, the deviatoric planes of the three models are presently coherent as plotted in Figure 12 (right).

Figure 13.

Reduction factor of tension meridian (left) and shear meridian (right).

The pressure limits of all three models are equivalent, bear in mind that CSCM constitutive equations are functions of the invariant $I_1$ thus $X_{\mathrm{d}}=3P_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}$ . The MHJC and CSCM crush volumetric strain parameter are the same $W={\mu }_{\mathrm{c}\mathrm{r}\mathrm{u}\mathrm{s}\mathrm{h}}$ . The rest of the pore collapse evolution parameters for the CSCM and RHT models were determined by inverse modelling the compaction and solidification phases of the hydrostatic test. The CSCM and RHT models are now capable of successfully describing all three phases of the pore collapse evolution (Figure 14 left). However, the bulk modulus of the solidification phase of the CSCM is lower than the MHJC model because the CSCM does not include any modification of the bulk modulus. The triaxial hydrostatic response of the RHT model matches perfectly with the MHJC model.

Figure 14.

Triaxial hydrostatic simulation (left) and ballistic limit curves (right).

Finally, the strain-rate sensitivity parameters of each model were identified from the CEB-FIB mode code (1990) and Malvar and Crawford (1998). The parameters involved in the damage formulae are excluded from Table 3, and their default values were adopted. The recalibrated parameters are given in the CSCM material card (see Table 4) and the RHT material card (see Table 5). Figure 14 (right) illustrates that the readjusted parameters of the CSCM and RHT models can effectively provide an accurate prediction of the impact experiments.

Table 4.

Material card of the CSCM.

*MAT_CSCM
$#	mid	ro	nplot	incre	irate	erode	recov	itretrc
	1	2.47E-9	7	1.293E-5	1.0	0.0	0.0	0
$#	pred
	0.0
$#	g	k	alpha	theta	lamda	beta	nh	ch
	17.5E3	23.33E3	83.8	0.0617	72.8	0.0092	1.0	0.0
$#	alpha1	theta1	lamda1	beta1	alpha2	thetal2	lamda2	beta2
	1.0	4.84E-5	0.2675	4.0E-4	1.0	1.58E-4	0.35	−7.6E-6
$#	r	xd	w	d1	d2
	5.0	72.8	0.057714	1E-4	1.0E-6
$#	b	gfc	d	gft	gfs	pwrc	pwrt	pmod
	100	12.09	0.1	0.1229	0.1229	5.0	1.0	0.0
$#	eta0c	nc	etaot	nt	overc	overt	srate	rep0w
	1.48E-4	0.739	1.48E-4	0.739	1.0E21	1.0E21	1.0	1.0

Table 5.

Material card of the RHT model.

*MAT_RHT
$#	mid	ro	shera	onempa	epsf	b0	b1	t1
	1	2.47E-9	17.5E3	−1.0	1.0	0.0	0.0	23.33E3
$#	a	n	fc	fs*	ft*	q0	b	t2
	1.62	0.5	72.8	0.33	0.1	0.65	0.035	0.0
$#	e0c	e0t	ec	et	betac	betat	ptf
	1.0E-5	1.0E-5	1.0E30	1.0E30	0.18	0.18	0.0
$#	gc*	gt*	xi	d1	d2	epm	af	nf
	1.0	1.0	1.0	0.04	1.0	0.01	1.62	0.5
$#	gamma	a1	a2	a3	pel	pco	np	alpha
	0.0	23.33E3	300.0E3	0.0	24.2667	500.0	1.0	1.055

Figure 15 illustrates fringe plots of the maximum volumetric strain ${\mu }_{\max }$ from each element at various penetration depths d, with $v_{\mathrm{i}}=157.1$ m/s and the MHJC concrete model. One can observe the formation of a crater on the impacted face of the target (d = 10 mm), followed by a tunnelling phase (d = 25 mm) through penetration. Next, scabbing appears at the rear face with diagonal cracks (d = 45 mm). Finally, the nose of the projectile perforates the concrete slab (d = 65 mm). Figure 15 (right) displays the concrete slab from the experiment and the maximum volumetric strain profile from the simulation superposed after full perforation. The simulation corresponds well with the experimental results.

Figure 15.

Fringe plot of the maximum volumetric strain ${\mu }_{\max }$ at selected projectile penetration depths d using the MHJC concrete model and $v_{\mathrm{i}}=157.1$ m/s, and a comparison of the simulation with the experiment after perforation by the projectile (right).

Concluding remarks

In this study we have evaluated automatic versus material test-based calibrations of concrete models for ballistic impact simulations. The first objective was to investigate the accuracy of the automatic parameter generators to produce input data for some commonly used material models. The second objective was to establish and evaluate a simplified parameter calibration procedure based on standard material experiments and data from the literature. The results employing automatic parameter generators varied notably between the models while still giving decent ballpark estimates. However, the parameters obtained from inverse modelling of standardized material tests improved the ballistic results significantly. Thus, the findings of this study recommend caution when using automatic parameter generators for ballistic impact simulations. Based on the work herein, the following main conclusions can be drawn:

• The MHJC model calibrated based on experimental data gave good predictions for the quasi-static compressive engineering stress-strain curve and for predicting the residual velocity of an ogive-nose projectile perforating a 50 mm thick C75 concrete slab.

• The CSCM, RHT and K&C models’ built-in parameter generators yielded reasonable results for the quasi-static engineering stress-strain curves from cylinder compression tests. There was, however, notable scatter in the estimated tensile capacity and in the residual velocity from the perforation simulations. Thus, the automatic generators should be used with caution.

• Adjusting the parameters of the CSCM and RHT models based on material tests significantly improved the numerical results for all simulated cases. It is thereby recommended to use at least some experimental data for calibration of these models.

• While dynamic and triaxial material tests are desirable for calibration, it is still possible to improve the accuracy of FE simulations significantly by using data from common material tests like cube compression, cylinder compression and tensile splitting.