A Dynamic Finite-Deformation Constitutive Model for Steels Undergoing Slip, Twinning, and Phase Changes

Abstract

Depending on composition, processing, microstructure, and loading mode, steels may demonstrate various inelastic deformation mechanisms, notably dislocation glide, deformation twinning, and solid-solid phase changes. Damage in the form of voids, often leading to failure by macro-cracking, is also of current interest. A finite-deformation constitutive model is constructed in order to address such mechanisms in isotropic polycrystals under static and dynamic loading, where the latter encompasses high pressures and extreme strain rates pertinent to ballistic penetration. Slip and twinning are isochoric, and their relative contributions to kinematics are not explicitly distinguished in the present application. Phase changes among, e.g., face-centered-cubic (FCC), body-centered-cubic (BCC), body-centered-tetragonal (BCT), and hexagonal (HCP) structures are admitted, with corresponding deviatoric and volumetric strains. Porosity from voids contributes to volumetric strain. A new consistent thermodynamic framework incorporating an Eulerian strain tensor and internal state variables is developed, whereby kinetic equations for slip–twinning, phase changes, and damage evolution result in contributions to dissipation. An objective rate form of the model, derived assuming small deviatoric elastic strain, is implemented numerically. The model is applied to three different primarily austenitic, medium-high Mn steels. Specifically, representations of a slip-dominated (SLIP), a TRIP (transformation-induced plasticity) and a TWIP (twinning-induced plasticity) steel are calibrated to quasi-static tension, quasi-static compression, and dynamic compression data, at both room and high temperatures. A novel functional form of material strength distinguishes hardening profiles of the different alloys. Experimental data are reasonably well represented. Model extrapolations for dynamic strength and pressure in regimes pertinent to shock compression are analyzed. Predictions for multi-axial loading of shear with simultaneous expansion or contraction demonstrate competing physical mechanisms among the alloys that could be leveraged for optimal ballistic performance.

Appendices

Appendix 1: Theoretical Derivations

Equation (19) is derived as

$$\begin{aligned} {{{\varvec{P}}}}:{\dot{{{\varvec{F}}}}}&= J \pmb {\sigma }:{{{\varvec{l}}}} = J \pmb {\sigma }:[ {{{\varvec{d}}}}^{E} + {{{\varvec{d}}}}^{P}] \\&= J \pmb {\sigma }:[ {{{\varvec{F}}}}^{E} {\dot{{{\varvec{D}}}}}^{E} {{{\varvec{F}}}}^{E \, {\text {T}}}+ {\bar{{{\varvec{d}}}}}^{P} + {\textstyle {\frac{1}{3}}} {\dot{J}}^{P} J^{P-1} \mathbf{1 }] \\&= J [({{{\varvec{F}}}}^{E \, {\text {T}}} \pmb {\sigma } {{{\varvec{F}}}}^{E}): {\dot{{{\varvec{D}}}}}^{E} + \pmb {\sigma }':{\bar{{{\varvec{d}}}}}^{P} - p \, {\text {tr}} {{{\varvec{d}}}}^{P}]. \end{aligned}$$

(97)

Equation (42) is derived as

$$\begin{aligned} J \sigma _{kk}&= \left[ B {D}^{E}_{\alpha \alpha } - {\textstyle {\frac{1}{2}}} B_{0} (B'_{0} -4) ({D}^{E}_{\alpha \alpha } )^{2} - A_{0} B_{0} (T-T_{0}) \right] \\&\quad {V}^{E \, -1}_{jk} {V}^{E \, -1}_{kj} + 2 G {F}^{E \, -1}_{ \beta k} {\bar{{D}}}^{E}_{\beta \gamma } {F}^{E \, -1}_{\gamma k} \\&= \left[ B {D}^{E}_{\alpha \alpha } - {\textstyle {\frac{1}{2}}} B_{0} (B'_{0} -4) ({D}^{E}_{\alpha \alpha } )^{2} - A_{0} B_{0} (T-T_{0}) \right] \\&\quad J^{E \, -2/3} {\bar{V}}^{E \, -2}_{kk} + 2 G J^{E \, -2/3} {\bar{U}}^{E \, -2}_{\beta \gamma } \bar{{D}}^{E}_{\beta \gamma } \\&\approx 3 J^{E \, -2/3} \\&\quad \left[ B {D}^{E}_{\alpha \alpha } - {\textstyle {\frac{1}{2}}} B_{0} (B'_{0} -4) ({D}^{E}_{\alpha \alpha } )^{2} - A_{0} B_{0} (T-T_{0}) \right] . \end{aligned}$$

(98)

Equation (43) is derived as

$$\begin{aligned} p&\approx -J^{P \, -1} J^{E \, -5/3} \\&\quad \left[ B {D}^{E}_{\alpha \alpha } - {\textstyle {\frac{1}{2}}} B_{0} (B'_{0} -4) ({D}^{E}_{\alpha \alpha } )^{2} - A_{0} B_{0} (T-T_{0}) \right] \\&\approx J^{P \, -1} B_{0} \\&\quad \left[ {\textstyle {\frac{3}{2}}} (J^{E \, -7/3}-J^{E \, -5/3}) \{ \zeta ^{B} - {\textstyle {\frac{3}{4}}} (B'_{0} -4)(1-J^{E \, -2/3} ) \} \right. \\&\qquad \left. + J^{E \, -5/3} A_{0} (T-T_{0}) \right] . \end{aligned}$$

(99)

Equation (44) is derived as

$$\begin{aligned} \pmb {\sigma } '&= \pmb {\sigma } + p {\mathbf{1 }} \approx 2 J^{-1} G \\&\quad \left[ {{{\varvec{F}}}}^{E \, -{\text {T}}} {\bar{{{\varvec{D}}}}}^{E} {{{\varvec{F}}}}^{E \, -1} - {\textstyle {\frac{1}{3}}} {\text {tr}} ({{{\varvec{F}}}}^{E \, -{\text {T}}} {\bar{{{\varvec{D}}}}}^{E} {{{\varvec{F}}}}^{E \, -1}) \right] \\&\approx 2 J^{-1} G \left[ {{{\varvec{F}}}}^{E \, -{\text {T}}} {\bar{{{\varvec{D}}}}}^{E} {{{\varvec{F}}}}^{E \, -1} - {\textstyle {\frac{1}{3}}} {\text {tr}} (J^{E \, -2/3} {\bar{{{\varvec{D}}}}}^{E} ) \right] \\&= 2 J^{-1} G {{{\varvec{F}}}}^{E \, -{\text {T}}} {\bar{{{\varvec{D}}}}}^{E} {{{\varvec{F}}}}^{E \, -1} . \end{aligned}$$

(100)

Equation (45) is derived as

$$\begin{aligned} \frac{\partial }{\partial t} {\bar{{{\varvec{D}}}}}^{E}({{{\varvec{X}}}},t)&= {{{\varvec{F}}}}^{E \,-1} {{{\varvec{d}}}}^{E} {{{\varvec{F}}}}^{E \, -{\text {T}}} - {\textstyle {\frac{1}{3}}} {\text {tr}} ( {{{\varvec{F}}}}^{E \,-1} {{{\varvec{d}}}}^{E} {{{\varvec{F}}}}^{E \, -{\text {T}}}) \mathbf{1 } \\&= {{{\varvec{F}}}}^{E \,-1} \\&\quad \left[ {{{\varvec{d}}}}^{E} - {\textstyle {\frac{1}{3}}} {\text {tr}} ( {{{\varvec{F}}}}^{E \,-1} {{{\varvec{d}}}}^{E} {{\varvec{F}}}^{E \, -{\text {T}}}) {{{\varvec{V}}}}^{E \,2} \right] {{{\varvec{F}}}}^{E \, -{\text {T}}} \\&\approx {{{\varvec{F}}}}^{E \,-1} \\&\quad \left[ {{{\varvec{d}}}}^{E} - {\textstyle {\frac{1}{3}}} ( {\bar{{{\varvec{V}}}}}^{E \,-2}: {{{\varvec{d}}}}^{E}) {\bar{{{\varvec{V}}}}}^{E \,2} \right] {{{\varvec{F}}}}^{E \, -{\text {T}}} \\&\approx {{{\varvec{F}}}}^{E \,-1} {\bar{{{\varvec{d}}}}}^{E} {{{\varvec{F}}}}^{E \, -{\text {T}}}. \end{aligned}$$

(101)

Equation (46) is derived as

$$\begin{aligned} {\dot{\sigma }}'_{ij}&= 2G J^{P \, -1} J^{E \, -7/3} {\bar{V}}^{E \, -1}_{ik} {\bar{V}}^{E \, -1}_{jm} {\bar{d}}^{E}_{km} - \sigma '_{ik} l^{E}_{kj} \\&\quad - l^{E}_{ki} \sigma '_{kj} + \sigma '_{ij} \left( \frac{{\dot{G}}}{G} - \frac{{\dot{J}}}{J} \right) \\&\approx 2G J^{P \, -1} J^{E \, -7/3} {\bar{d}}^E_{ij} - \sigma '_{ik} l^{E}_{kj} - l^{E}_{ki} \sigma '_{kj} \\&\quad + \sigma '_{ij} \left( \frac{1}{\zeta ^{G}} \frac{\text {d}\zeta ^{G}}{\text {d}d} {\dot{d}} - \frac{{\dot{J}}}{J} \right) . \end{aligned}$$

(102)

Equation (72) is derived as

$$\begin{aligned} {\mathfrak {D}}^{tr}&= \pmb {\sigma }':{\bar{{{\varvec{d}}}}}^{tr} - p \frac{{\dot{J}}^{tr}}{J^{tr}} + f^{\upsilon } {\dot{\upsilon }} \\&= {\bar{\sigma }} {\dot{\epsilon }}^{tr} - p {\dot{{\varDelta }}}^{tr} + (f^{{\varLambda }} + f^{{\varUpsilon }}) {\dot{\upsilon }} \\&= {\bar{\sigma }} \left[ \frac{1}{2} A^{tr} \gamma ^{tr} + \sqrt{\frac{3}{2}} \frac{{\varSigma }\delta ^{tr}}{1+\delta ^{tr} \upsilon } \right] {\dot{\upsilon }} \\&\quad + (f^{{\varLambda }} + f^{{\varUpsilon }}) {\dot{\upsilon }}\\&= {\bar{\sigma }} \alpha ^{tr}{\dot{\upsilon }} - \left[ {\varLambda }+ \frac{{\varUpsilon }_{0}}{l_{0}} \right] {\dot{\upsilon }} \\&= [f^{mech} - f^{th} - f^{ath} ]{\dot{\upsilon }}, \end{aligned}$$

(103)

where

$$\alpha ^{tr}({\varSigma },\upsilon ) = \frac{1}{2} A^{tr} \gamma ^{tr} + \sqrt{\frac{3}{2}} \frac{{\varSigma }\delta ^{tr}}{1+\delta ^{tr} \upsilon }.$$

(104)

Equation (79) is derived as

$$\begin{aligned} T (\partial \eta / \partial \pmb {\xi }) \cdot {\dot{\pmb {\xi }}}&= - T (\partial \pmb {\zeta } / \partial T ) \cdot {\dot{\pmb {\xi }}} = -T (\partial f^{\upsilon } / \partial T) {\dot{\upsilon }} \\&= - T (\partial f^{{\varLambda }} / \partial T) {\dot{\upsilon }} = T (\partial {\varLambda }/ \partial T) {\dot{\upsilon }} \\&= \lambda _{T} \frac{T}{T_{T}} {\dot{\upsilon }} = \lambda (T) {\dot{\upsilon }}. \end{aligned}$$

(105)

Equation (90) is derived as

$$\begin{aligned} {\mathfrak {D}}^{d}&= -\frac{p}{J^{d}} {\dot{J}}^{d} + f^{d} {\dot{d}} = -\frac{p}{1-d} {\dot{d}} -\frac{\partial W}{\partial d} {\dot{d}} \\&= \left[ \sqrt{\frac{3}{2}} \frac{ {\varSigma }{\bar{\sigma }}}{1-d} \right. \\&\quad \left. - \left( \frac{\partial \zeta ^{B}}{ \partial d} \cdot \frac{1}{2} B_{0} ({\text {tr}} {{{\varvec{D}}}}^{E})^{2} + \frac{\text {d}\zeta ^{G}}{ \text {d}d} \cdot G_{0} {\bar{{{\varvec{D}}}}}^{E}:{\bar{{{\varvec{D}}}}}^{E} \right) \right] {\dot{d}}\\&\ge 0. \end{aligned}$$

(106)

Appendix 2: Data Analysis

In most cases, 2–3 experiments were performed for each alloy at each initial temperature and loading rate. For every loading protocol, the number of experiments is shown in Table 5 along with experimental variation in effective stress, model prediction of effective stress, and the error in this prediction relative to the experimental mean. Comparisons are made at two discrete strain levels $\epsilon$ in each case, corresponding to availability of test data (e.g., some compression experiments were unloaded prior to attainment of very large deformations). Stress variability is defined as ${\varDelta }\sigma = {\text {max}}_{i} | {\bar{\sigma }} - \sigma _{i} |,$ where ${\bar{\sigma }}$ is the mean among tests and $\sigma _{i}$ is the measured value for experiment i,  where i runs from 1 to the number of tests shown in corresponding rows of the table. Modeling error in the rightmost column is defined, at each $\epsilon ,$ as $|\sigma - {\bar{\sigma }} | / {\bar{\sigma }} ,$ with the model result denoted simply by $\sigma .$ Agreement with tensile data is generally superior to that with compression data since the calibration procedure of “Static and Dynamic Loading” section focuses on the tensile response. Of the 36 cases quantified in Table 5, modeling error exceeds 4% only 11 times. The most inaccurate results, with error exceeding 13%, correspond to static compression of the TWIP steel at room temperature (noting that the model, calibrated to tensile response, fails to capture tension–compression asymmetry in the TWIP alloy), and to dynamic high-temperature compression of the TWIP steel, for which thermal softening is apparently not well represented.

Table 5 Experimental variation and model accuracy at two applied strain levels $\epsilon$ for each loading protocol: experimental range ${\varDelta }\sigma$ (max difference between average and max/min among number of tests), model prediction of $\sigma ,$ and error in model prediction of $\sigma$ at corresponding $\epsilon$

Since experimental data on the three alloys of Table 1 for compression and tension are available, but torsion data are notably absent, model predictions for all three stress states and all three steels are compared versus one another and versus trends for several other steels reported elsewhere. Results for the SLIP and TWIP steels are shown in Fig. 13a. These show negligible effect of triaxiality on stress–strain behavior, with the exception of tensile softening due to porosity that initiates at $\epsilon \gtrsim 0.35$ in the SLIP steel. In each case, effective stress $\sigma$ is von Mises stress, and effective strain is $\epsilon = \sqrt{2/3} \int {{{\varvec{d}}}}: {{{\varvec{d}}}} \text {d}t.$ For tension and compression, these definitions are consistent with true stress and logarithmic strain measures used in “Static and Dynamic Loading” section. For free-end torsion, $\sigma = \sqrt{3} \tau$ and $\epsilon = \gamma / \sqrt{3},$ with $\tau$ and $\gamma$ the usual engineering shear stress and strain measures.

Fig. 13

Model predictions of effect of loading mode on effective stress $\sigma$ under quasi-static, room temperature conditions a SLIP and TWIP alloys and b TRIP alloy

Results in Fig. 13b demonstrate a marked effect of stress state on strength behavior of the TRIP steel. For strain-assisted transformation, stress in tension very slightly exceeds that in torsion, since the contribution of tensile pressure is not enough to drastically accelerate phase changes in the former relative to the latter. In contrast, martensitic transformation is impeded by compressive stress, leading to a softer response. Stress-assisted transformation corresponds to the results shown in Fig. 3a, and is thought more realistic than strain-assisted transformation for static tensile loading. In this case, strain accommodation by phase changes initiated prior to plastic yield induces a load reduction, which in turn retards the transformation rate at low applied strains.

Trends predicted in Fig. 13 are similar, but not identical, to those reported for RHA steel in [59]. Quasi-static room temperature data for RHA show very similar plastic strength behavior for tension, compression, and shear, at least until damage accumulates. At larger strain levels, porosity nucleated at hard particles was more severe in torsion than in tension, leading to softer torsional behavior and earlier failure [59]. This behavior cannot be captured by the traditional Cocks–Ashby type damage model (89), used in the present work, which requires positive mean stress for porosity increases.

Data and/or modeling on 304L stainless steel in [77] show tensile strength exceeding compressive strength, with torsional hardening the lowest among the three stress states. Expansive martensitic transformation from the $\gamma$ to $\alpha '$ phase was reported to be greater in compression loading than torsional loading, despite the negative contribution to mechanical driving force induced by pressure in the former. Less localized slip and different dislocation structures observed in torsion than compression were suggested as responsible for the observed stress-state dependencies [77]. If relative torsional softening occurs in any of the three alloys of present study—which must be confirmed or refuted by future torsion experiments—adjustments to the present model framework are necessary to capture such behavior. For example, a $J_{3}$-based yield function is advocated in [29, 77], which provides satisfactory description of their experimental results.

Table 6 compares experimental data and model calibrations for quasi-static tensile failure behavior of the three alloys in Table 1, at room and elevated temperatures. The logarithmic tensile strain to failure is $\epsilon _{F}.$ This is the average ductility among experimental tests, while it is the exact failure strain imposed in calibration of porosity $\phi _{F}$ for the damage model. The experimental variation of ductility is ${\varDelta }\epsilon _{F} = {\text {max}}_{i} | \epsilon _{F} - \epsilon _{i} |$ with $\epsilon _{i}$ the strain to failure of each experiment i at the given temperature for the alloy under consideration. Similarly, experimental variation of tensile strength at failure is ${\varDelta }\sigma _{F} = {\text {max}}_{i} | {\bar{\sigma }}_{F} - \sigma _{i} |,$ where ${\bar{\sigma }}_{F}$ is the average over experiments i with individual tensile strengths $\sigma _{i}.$ The predicted tensile strength of the model at $\epsilon _{F}$ is $\sigma _{F},$ and the error in this prediction is $|\sigma _{F} - {\bar{\sigma }}_{F} | / {\bar{\sigma }}_{F}.$ In order to achieve perfect agreement in $\epsilon _{F}$ between model and experiment, a composition- and temperature-dependent failure porosity $\phi _{F}$ must be assigned in each case. Determination of whether or not such values are physically descriptive requires microscopy of fractured surfaces as in [59], for example. Agreement between model and experiment for tensile strength at $\epsilon _{F}$ is within 8% error, and within 4.5% when high temperature data for TWIP and TRIP steels are excluded.

Table 6 Tensile failure data: average measured strain to failure $\epsilon _{F}$ with range ${\varDelta }\epsilon _{F}$ (max difference between average and max/min among number of tests), experimental range of tensile strength ${\varDelta }\sigma _{F},$ model prediction of $\sigma _{F}$ at $\epsilon = \epsilon _{F},$ model porosity at failure $\phi _{F}$ to achieve same strain to failure, and error in model prediction of $\sigma _{F}$