Macro-scale modeling of finite strain viscoplasticity in irradiated F/M steels: a continuum thermodynamic framework

Abstract

The European reduced activation ferritic/martensitic steel Eurofer97 in irradiated states displays significant nonlinear material behavior involving irradiation hardening, loss of strain hardening, and uniform elongation, as well as irradiation-induced embrittlement. Nonlinear behavior of irradiated steel modeled on the continuum scale will help to estimate the maximum operating range of the irradiated components beyond the onset of localized plastic flow. In this work, a thermodynamic framework for modeling irradiation-influenced deformation is established based on irradiation defect density and a thermodynamically consistent finite strain formulation of an existing viscoplastic model using the Dual Variables concept is presented. The model is implemented in ABAQUS allowing the simulation of tensile tests conducted on irradiated and unirradiated materials which shows the model’s ability to capture the post-yield and post-necking behavior observed in experiments up to ductile failure.

1 Introduction

The reduced activation ferritic martensitic (RAFM) steel Eurofer97 is used to make the European test blanket modules (TBM) for ITER and the breeding blankets of future fusion reactors like the DEMOnstration fusion reactor DEMO. Eurofer97 (Fe–9Cr–1.1W–0.2V–0.12Ta) has good thermomechanical properties and resists irradiation-induced swelling [1, 2]. After low-temperature neutron irradiation, the material exhibits strong material hardening and a reduction in uniform and total elongation [3, 4]. Irradiated material has higher yield-stress and ultimate-tensile-stress (UTS) than unirradiated material, but softens rapidly after UTS, reducing fracture strain. Irradiation-induced hardening is caused by dislocation pinning by irradiation defects, while softening is caused by dislocation movement and interaction with defects [5]. Dislocations move forward and clear defects in localized areas when a critical stress is reached. This process, known as channeling, creates channels of soft material surrounded by defect-rich hard material [6, 7]. Dislocations at grain boundaries interact, nucleating voids that grow and coalesce into microcracks. Channeling is the primary deformation mechanism, while void coalescence is the primary damage mechanism [6]. Loss of uniform elongation is a component failure criterion in the nuclear industry [8]. Given the high cost of replacing components and the probable reactor downtime, it is important to precisely predict their maximum operating range particularly when uniform elongation/deformation is lost, e.g., due to irradiation. Continuum techniques allow modeling at the component scale, which describes post-yield and flow softening behavior. This would help to estimate component life and build enhanced breeding blanket design criteria. Existing micromechanical continuum-based models use crystal plasticity to describe material properties for single crystals [9, 10] and polycrystals [11], and self-consistent plasticity theory [12]. These approaches combine microstructure observations with computationally expensive simulations at different scales to capture dominant physical mechanisms at a finer length scale and pass this information to a higher length scale model until a coarse FE-mesh can predict the macroscopic mechanical properties with physical fidelity. More complex multiscale techniques may be essential when addressing damage initiation and grain boundary effects. They also require very large-scale simulations to mimic necking due to their mesh dependency to suit the microstructure-related parameters set on the nanometer scale. These methodologies are more appropriate for fundamental investigation, but too expensive and impractical for engineering analysis, which requires economical and quick models for components much larger than a tensile sample.

This work uses constitutive equations based on established material theory for continuum scale modeling and can be calibrated with mechanical tests. This approach is important for design and development because it can be used with common FEM codes to simulate component deformation and damage under operating conditions. A viscoplastic framework is used to define a well-posed initial value problem (IVP) and keep field equations hyperbolic, allowing faster calculations and reduced mesh dependency. The elasto-viscoplastic model from [13] proposes constitutive equations for defect nucleation and annihilation to model irradiation’s small-strain effect (below 5%). In this work, the constitutive system of equations by [13] is extended to the finite strain framework, where large deformations and rotations can be handled. The Dual Variables approach [14] is used to ensure a thermodynamically consistent finite strain model formulation. Currently, there is no thermodynamic framework to account for the changes in the internal energy of irradiated material due to interactions between energetic neutrons and lattice atoms and due to the energy dissipation during defect annihilation by plastic deformation. In this work, a modified thermodynamic framework is proposed to account for this change in internal energy, with defect density introduced as an internal variable in the material’s free energy function.

The paper is organized as follows: Sect. 2 discusses finite strain and Dual Variables. Section 3 presents the thermodynamic framework for modeling irradiated materials and the finite strain model. Also covered is model integration. Section 4 compares tensile test results on unirradiated and irradiated materials with model simulations. Section 5 discusses the results, and Sect. 6 concludes.

2 Finite strain framework

The model developed relies on basic theoretical concepts related to the kinetics of finite strain framework and the Dual Variables approach, which are introduced briefly here. Let us consider a continuum body, which consists of an infinite number of material points, undergoing deformation. In continuum mechanics, the vector ${\varvec{X}}$ is used to define the reference configuration $R_{r}$ and$\, {\varvec{x}}$ is the vector that defines the actual configuration $R_{a}$, with the motion of the body ${\varvec{x}}=\,\mathcal {X}\left( {\varvec{X}},t\right)$ during time t is described using the mapping function $\mathcal {X}\left( {{\varvec{X}}},t \right)$. The derivative of $\mathcal {X}$ with respect to ${\varvec{X}}$ provides a second-order tensor known as the deformation gradient,

$$\begin{aligned} {\varvec{F}}=\frac{\partial \mathcal {X}\left( {\varvec{X}},t \right) {\, }}{\partial {\varvec{X}}}=\, \frac{\partial {\varvec{x}}}{\partial {\varvec{X}}}. \end{aligned}$$

(1)

The polar decomposition of ${\varvec{F}}={\varvec{RU}}={\varvec{VR}}$yields the symmetric right and left stretch tensors ${\varvec{U}}$ and V along with a rotation tensor ${\varvec{R}}$. The spatial velocity gradient tensor ${\varvec{L}}=\dot{{\varvec{F}}}{\varvec{\,}}{\varvec{F}}^{{-}1}$ can be decomposed cumulatively into the symmetric deformation rate tensor ${\varvec{D}}$ and the skew-symmetric spin tensor ${\varvec{W}}$ using

$$\begin{aligned} {\varvec{D}}=\frac{1}{2}\, \left( {\varvec{L}}+{\varvec{L}}^{\, T} \right) \, \, \, ,\, \, \, {\varvec{W}}=\frac{1}{2}\, \left( {\varvec{L}}-{\varvec{L}}^{\, T} \right) . \end{aligned}$$

(2)

The multiplicative decomposition of ${\varvec{F}}={\varvec{F}}_{{e}}\, {\varvec{F}}_{p}$ into the elastic and plastic components allows the definition of the intermediate configuration, $\hat{R}_{a}$, which is a stress-free state resulting from the local unloading process. The invariance postulation of [15] demands the invariance of constitutive relationships against any rotations of the intermediate configuration. The Dual Variables approach [14] proposes that a strain tensor and its best-associated stress tensor, the corresponding stress power, have invariance properties. The stress power per unit volume in $R_{r}$ is $W={\varvec{T}}\cdot \dot{{\varvec{E}}}$, where ${\varvec{E}}=\left( \, {\varvec{F}}^{T}{\varvec{F}}-{\varvec{I}} \right) {\varvec{/}}2\mathrm {\, }$ is the Green Lagrange strain and ${\varvec{I}}$ is the second-order unity tensor. The associated stress tensor is ${\varvec{T}}$, the second Piola-Kirchhoff stress tensor. Stress and strain tensors, which are naturally associated in this way, are called Dual Variables. The transformation of the strain tensor [14] along, $R_{r}\rightarrow \hat{R}_{a}\rightarrow R_{a}\,$ is illustrated in Fig. 1, which is given by

Fig. 1

Deformation decomposition and internal variables of the finite strain model

$$\begin{aligned} {{{\varvec{E}}} \, \, \xrightarrow []{ {{\varvec{F}}}_{p}^{T-1} \, \, {{\varvec{E}}} \, \, {{\varvec{F}}}_{p}^{-1}}\, \, \, \hat{\varvec{\varGamma } } \, \, \, \xrightarrow []{ {{\varvec{F}}}_{e}^{T-1} \, \,{\hat{\varvec{\varGamma } }} \, \, {{\varvec{F}}}_{e}^{-1}}\, \,{{\varvec{A}}}} \end{aligned}$$

(3)

where the transformation $R_{r}\rightarrow \hat{R}_{a}$ is inelastic and $\hat{R}_{a}\rightarrow R_{a}$ is elastic. Here,$\, \hat{\varvec{\varGamma } }$ is the strain tensor transformed onto $\hat{R}_{a}$ and ${\varvec{A}}$ is the Almansi strain tensor acting on $R_{a}\mathrm {.}$ The corresponding dual-stress tensors with respect to $\hat{R}_{a}$ and $R_{a}$ are the second Piola-Kirchhoff stress tensor $\hat{{\varvec{T}}}$ and the Cauchy stress tensor ${\varvec{S}}={\varvec{F}}_{e}{\varvec{\, }}\hat{{\varvec{T}}}{\varvec{\, }}{\varvec{F}}_{e}^{T}$, respectively. Hence, the invariant stress power in the various configurations can be summarized and written as

$$\begin{aligned} W=\underbrace{{\varvec{T}}\cdot \dot{{\varvec{E}}}}_{Reference}= \underbrace{{\hat{{\varvec{T}}}\cdot } \mathop {{\hat{\varvec{\varGamma }}}}\limits ^{\varDelta }}\limits _{Intermediate}= \mathop {\underbrace{{\varvec{S}} \cdot {\varvec{D}}}}\limits _{Actual}. \end{aligned}$$

(4)

The time derivatives of the strain tensors are given by the lower Oldroyd time derivative [16] of the form

$$ \begin{aligned} \left( \, \, \right) ^{\varDelta }=\dot{\left( \, \, \right) }+\, \hat{{\varvec{L}}}_{p}^{T}\, \left( \, \, \right) +\left( \, \, \right) \hat{{\varvec{L}}}_{p}\, \, \, \, \textrm{for}\, \, \, \hat{R}_{a}\, \, \, \mathrm { \& \, \, }\, \left( \, \, \right) ^{\varDelta }=\dot{\left( \, \, \right) }+{\varvec{\, }}{\varvec{L}}^{{T}}{\varvec{\, }}\left( {\varvec{\, \, }} \right) {\varvec{+}}\left( {\varvec{\, \, }} \right) {\varvec{L\, }}\, \, \, \textrm{for}\, \, \, R_{a}. \end{aligned}$$

(5)

The Oldroyd time derivative of ${\varvec{A}}$ is the deformation rate tensor ${\varvec{D}}$ $(\overset{ {\varDelta }}{{\varvec{A}}} = \varvec{\dot{A}} - {\varvec{L}}{\varvec{A}} - {\varvec{A}}{\varvec{L}}^{T} = {\varvec{D}})$ and that of the plastic strain $\hat{\varvec{\varGamma }}_{p}$ is the plastic part of the deformation rate $\hat{{\varvec{D}}}$

$$\begin{aligned} \left( \hat{\varvec{\varGamma }}_{p} \right) ^{\varDelta }=\dot{\hat{\varvec{\varGamma } }}_{p}+\, \hat{\varvec{L}}_{p}^{\, T}\, \left( \hat{\varvec{\varGamma } }_{p} \right) +\left( \hat{\varvec{\varGamma } }_{p} \right) \, \hat{\varvec{L}}_{p}=\hat{\varvec{D}}_{p}.\, \end{aligned}$$

(6)

The material time derivatives of stress tensors are given by the upper Oldroyd derivatives

$$\begin{aligned} \left( \, \, \right) ^\nabla =\dot{\left( \, \, \right) }-\, \hat{{\varvec{L}}}_{p}^{T}\, \left( \, \, \right) -\left( \, \,\right) \hat{{{\varvec{L}}}}_{p}\, \, \, \, \textrm{for}\, \, \, \hat{R}_{a},\, \, \left( \, \, \right) ^\nabla =\dot{\left( \, \, \right) }-\,{\varvec{L}}^{T}{\varvec{\, }}\left( {\varvec{\, \, }} \right) -\left( {\varvec{\, \, }} \right) {\varvec{L\, }}\, \, \, \textrm{for}\, \, \, R_{a}. \end{aligned}$$

(7)

3 Modelling of irradiated materials

In this section, the thermodynamic framework to model an irradiated metal is constructed. The thermodynamic formulation of the finite strain model and consistency criteria are derived using this framework. A simplified model for small elastic deformations is developed for implementation in FE programs.

3.1 Thermodynamics

The local form of the first law of thermodynamics [16] for an unirradiated material is

$$\begin{aligned} \dot{e}=\frac{1}{\rho }\, {\varvec{T}}\cdot \dot{{\varvec{E}}}-\frac{1}{\rho }\, \nabla \cdot {\varvec{q}}+s \end{aligned}$$

(8)

where e is the specific internal energy, ${\varvec{q}}$ is the heat flux vector, s is the specific heat per unit time, $\rho$ is the mass density. Thermodynamic potentials and the system’s internal energy can be described by a pair of conjugate variables, typically a force and the resulting system displacement. Alterations to a mechanical system’s energy are represented by the product of quantities from these pairs. The stress power, ${\varvec{T}}\cdot \dot{{\varvec{E}}}$, describes the transfer of mechanical or dynamic energy resulting from work done.

A radiation damage event is the transfer of energy from high-velocity particles (neutrons) to a solid, causing atom reorganization [17]. For crystalline materials, radiation damage is measured in displacements per atom (dpa) [18]. Most of the absorbed energy is dissipated as heat, but the rest is consumed to create defects in the lattice [18] such as voids, defect clusters, dislocation loops, etc., at the cost of formation energy unique to each defect type. Tensile tests demonstrate that these defects cause irradiation hardening, loss of ductility, etc. Irradiation-induced hardening, generated by the pinning of dislocations by defects, is indicated by $\sigma _{H}$ [4] and is a function of the irradiation defect density N, given by, $\sigma _{H}=h\, \sqrt{N}$ [19]. h is a material-dependent parameter which varies with temperature. Plastic deformation removes defects and softens irradiated materials [5, 20, 21]. In this work, it is postulated that the system absorbs kinetic energy from incident neutrons to form defects and releases it during plastic deformation. The material’s thermodynamics must account for all these exchanges to guarantee energy conservation. The energy can be calculated by adding the products of defect population and defect formation energy for each type of defect. Due to a lack of reliable methods for precisely estimating the types of defects and their population for a given irradiation dose, this is currently impossible. Alternatively, this work proposes a conjugate pair of variables whose product quantifies the energy transferred to the system. This energy, together with the material’s stress power, is dissipated during inelastic deformation. The extrinsic quantity of irradiation damage dose, measured in dpa, can be employed to measure the thermodynamic change. The thermodynamic force needed to bring about this change is denoted as $\xi$. It is an intrinsic quantity measured in MPa/dpa and a function of $\phi$. The product of this conjugate pair,

$$\begin{aligned} P_{irr}\left( \phi \right) =\xi \, \dot{\phi }, \end{aligned}$$

(9)

is a power term measured in MPa/s or J/s (considering volume) and is therefore referred to as the irradiation-induced power in this work. Introducing $P_{irr}$ into first law (8) satisfies the law of energy conservation

$$\begin{aligned} \dot{e}=\frac{1}{\rho }\, \left( {\varvec{T}}\cdot \dot{{\varvec{E}}}+\xi \, \dot{\phi }\, \right) -\frac{1}{\rho }\, \nabla \cdot {\varvec{q}}+s_{total}.\, \end{aligned}$$

(10)

Here, $s_{total}=s+s_{irr}$ where $s_{irr}$ represents the heat generated due to the neutron–atom and atom–atom interactions as a result of irradiation. In this form, the change in internal energy results not only from the mechanical work and heat exchange in the control volume but also from the work done in the creation of irradiation defects. The corresponding Clausius–Duhem entropy inequality for Eq. (10) for positive entropy production is

$$\begin{aligned} -\, \dot{\Psi }\, -\dot{T}\, \eta \, +\frac{1}{\rho }\, \left( \, {\varvec{T}}\cdot \dot{{\varvec{E}}}+\xi \, \dot{\phi } \right) \, -\, \frac{{\varvec{q}}}{\rho \, T}\cdot \nabla T\, \ge \, 0. \end{aligned}$$

(11)

Here, $\Psi =e-\, T\, \eta$ is the specific free energy, $\eta$ is the specific entropy, and T is the absolute temperature. Some defects, like point defects, do not contribute to material hardening, whereas others do. Although eliminating any defect type dissipates energy, this work will focus on defects that contribute to material hardening. The contribution to free energy from defects that do not affect hardening is removed from $\Psi$ and defined as $\Psi _{irr,\, NH}$. For an isothermal process with uniform heat distribution, the inequality reduces to

$$\begin{aligned} -\rho \, \dot{\Psi }+\, {\varvec{T}}\cdot \dot{{\varvec{E}}}+\left( \, \xi \, \dot{\phi }-\, \rho \, \dot{\Psi }_{irr,\, NH} \right) \ge 0. \end{aligned}$$

(12)

Subtracting $\dot{\Psi }_{irr,\, NH}$ from the irradiation-induced power leaves us with the portion of power responsible for hardening that is available for dissipation during plastic deformation. To determine the dissipating power term, we propose the following relation that connects the irradiation-induced power to the irradiation hardening, $\sigma _{H}$ and neutron dose, $\phi$ through the term $\varphi$, which is proposed to be the tangent of the $\sigma _{H}-\phi$ curve

$$\begin{aligned} \xi \, \dot{\phi }-\, \rho \, \dot{\Psi }_{irr,\, NH}=\, \, {\varphi }\, \dot{\phi },\, \, \varphi =\, \frac{d\sigma _{H}}{d\phi }. \end{aligned}$$

(13)

The revised entropy inequality,

$$\begin{aligned} -\rho \, \dot{\Psi }+\, {\varvec{T}}\cdot \dot{{\varvec{E}}}+\varphi \, \dot{\phi }\ge 0, \end{aligned}$$

(14)

is used to determine the finite strain model’s constitutive relations and thermodynamic consistency criteria.

Table 1 Conjugate pairs of variables

3.2 Thermodynamic formulation of finite strain model

The constitutive model is defined on $\hat{R}_{a}$ to be invariant to rotations [22, 23]. First, the strain tensor $\hat{\varvec{\varGamma } }$ is decomposed into its elastic and plastic components

$$\begin{aligned} \mathrm {\, }\hat{\varvec{\varGamma } }\mathrm {=\, } \hat{\varvec{\varGamma } }_{e}+\, \hat{\varvec{\varGamma } }_{p}\mathrm {.} \end{aligned}$$

(15)

The Helmholtz free energy $\Psi$ of the irreversible solid system is defined as a function of a set of internal state variables (ISV) which are chosen based on physically observed behavior as

$$\begin{aligned} \Psi \left( t \right) =\, \hat{\Psi }\left( \hat{\varvec{\varGamma } }_{e},\, \hat{\varvec{\alpha } }_{i}, \, \psi _{i}, N_{i}, \, d \right) . \end{aligned}$$

(16)

Table 1 lists the chosen set of ISVs and their conjugates with respect to each configuration. Among these, the irradiation defect-density, $N_{i},$ and its conjugate, ${\mathcalligra {n}}_{i}$, are proposed to be included to account for the irradiation-induced hardening. The elastic and plastic components of the free energy function are given in Table 2. Using Eqs. (4) and (14), the entropy inequality at $\hat{R}_{a}$ changes to

$$\begin{aligned} \hat{\varvec{T}}\cdot \mathop {\hat{\varvec{\varGamma }}}\limits ^\varDelta +\varphi \, \dot{\phi }-\rho \, \dot{\Psi }\ge 0. \end{aligned}$$

(17)

The partial derivatives of $\Psi \mathrm {\, }$with respect to its ISVs provide the following relations, resulting in the decomposition of the entropy inequality into two parts, describing the elastic and inelastic material behavior

Table 2 Free energy functions

$$\begin{aligned} \hat{{\varvec{T}}}=\rho \frac{\partial \Psi _{e}}{\partial \hat{\varvec{\varGamma } }_{e}\, },\, \, \, \, \, \, \, \hat{{\varvec{Z}}}_{i}\, =\rho \frac{\partial \Psi _{p}^{Kin}}{\partial \hat{\varvec{\alpha } }_{i}\, },\, \, \, \, \, \, r_{i}=\rho \frac{\partial \Psi _{p}^{Iso}}{\partial \psi _{i}\, },\, \, \, \, \, \, n_{i}=\rho \frac{\partial \Psi _{p}^{Irr}}{\partial N_{i}\, },\, \, \, \, \, \, Y=\rho \frac{\partial \Psi _{p}^{Dam}}{\partial d\, }.\, \end{aligned}$$

(18)

The elastic behavior is determined from Eq. (17) using the first term in Eq. (18) and $\Psi _{e}\,$from Table 2. For describing the inelastic deformation, an effective measure of stress, $\hat{{\varvec{T}}}/\psi (1-d),$ is adopted in place of $\hat{{\varvec{T}}}$, considering the relaxation due to the isotropic softening and damage. Simplification of Eq. (17) leads to the definition of the Mandel-type stress tensor,

$$\begin{aligned} \hat{{\varvec{P}}}=\left( {\varvec{I}}+2\hat{\varvec{\varGamma }}_{{e}}\right) \, \hat{{\varvec{T}}}\, \end{aligned}$$

(19)

which reduces to $\hat{{\varvec{T}}}$ under the assumption that elastic deformation is minimal in finite strain plasticity, $\hat{\textbf{V}}_{{e}}\mathrm {\approx }{\varvec{I}}$, $\hat{\varvec{\varGamma } }_{e}\mathrm {\approx }{\varvec{0}}$ [23]. Using Eqs. (17), (18) and the free energy functions from Table 2, the plastic part of entropy inequality (17) is obtained as

$$\begin{aligned} \underbrace{\left( \frac{\hat{{\varvec{T}}}}{\varPsi (1-d)}-{\varvec{\hat{\Omega }}}_{i}\right) \cdot \mathop {{\varvec{\hat{D}}}}\nolimits _p+\varphi \dot{\phi }-r_{i}\dot{\psi _i}-{\mathcalligra {n}}_{i}\, \dot{N}_{i}}_{Isotropic}+ \underbrace{\left( ({{\varvec{I}}}+2{\varvec{\hat{\varPhi }}}){\varvec{\hat{D}}}_p- \dot{\hat{\varvec{\alpha }}}_{i}\right) \cdot {\varvec{\hat{Z}}}_i}_{Kinematic} -\quad \underbrace{Y\,\,\dot{d}}_{Damage}\ge 0 \end{aligned}$$

(20)

where $\hat{\varvec{\Omega } }_{i}=\left( {\varvec{I}}+2\, \hat{\varvec{\varPhi } } \right) \, \hat{{\varvec{Z}}}_{i}$ is a Mandel-type kinematic hardening tensor with $\hat{\varvec{\varPhi } }$ being a strain-type tensor.

The viscoplastic yield potential $\Omega$ of the model [24] is based on the over-stress tensor $\hat{\varvec{\varSigma }}=\hat{{\varvec{T}}}^{D}/\psi (1-d)-\hat{\varvec{\Omega } }$ with $\hat{{\varvec{T}}}^{{{D}}}= \, \hat{{\varvec{T}}}\, - \, tr\left( \hat{{\varvec{T}}} \right) {{\varvec{I}}}$ being the deviatoric part of the stress tensor,

$$\begin{aligned} \Omega =\, \, \mathrm {\varSigma }_{eq}-\sigma _{H}-k \end{aligned}$$

(21)

where $\mathrm {\varSigma }_{eq}$ is the von Mises equivalent of $\hat{{\varvec{\varSigma }} }$, $\sigma _{H}$ is the irradiation hardening variable, and k is the initial yield stress of the unirradiated material. Applying the normality rule, we obtain the flow rule and the inelastic strain rate

$$\begin{aligned} \overset{\varDelta }{\hat{\varvec{\varGamma }}_{p}}= & {} \hat{{\varvec{D}}}_{p}=\, \frac{3}{2}\left\langle \frac{\varSigma _{eq}-\sigma _{H}-k\, }{K} \right\rangle ^{n}\frac{\hat{{\varvec{\varSigma }} }}{\varSigma _{eq}},\, \, \end{aligned}$$

(22)

$$\begin{aligned} \dot{p}= & {} \sqrt{\frac{2}{3}\, \hat{{\varvec{D}}}_{p}\cdot \hat{{\varvec{D}}}_{p}} =\left\langle \frac{\varSigma _{eq}-\sigma _{H}-k\, }{K} \right\rangle ^{n}\, \end{aligned}$$

(23)

where p is the accumulated plastic strain variable and $n,\, K$ are the temperature-dependent viscous parameters. The McAuley brackets, $\langle {}\rangle$ operate as: $\left\langle x \right\rangle =\frac{\, x+\left| x \right| }{2}$. The kinematic hardening tensor should be handled with the finite strain framework. Here, the Chaboche adoption of the nonlinear hardening rule by Armstrong and Frederick [25] is used. Its relationship with its conjugate variable is obtained using Eq. (18) and $\Psi _{p}^{Kin}$ from Table 2 as

$$\begin{aligned} \hat{{\varvec{Z}}}_{i}=A\left( p \right) \, \, \hat{{\varvec{\alpha }} }_{i} \end{aligned}$$

(24)

where $A\left( p \right) =C\, a\, f\left( p \right)$ is a plasticity-dependent function with $C\,$and a as temperature-dependent parameters. The thermodynamic consistency is achieved by deriving the following Oldroyd time derivative of $\hat{{\varvec{Z}}}$ using the kinematic part of the entropy inequality, which is discussed in appendix:

$$\begin{aligned} \mathop {\hat{{\varvec{Z}}}}\limits ^{{\nabla }} =C\, f\left( p \right) \left\{ a~\hat{{{\varvec{D}}}}_{p}+2\, \hat{{\varvec{D}}}_{p}{{\varvec{\, }}}\left( \hat{{\varvec{\varPhi }}} -\, \hat{{\varvec{\alpha }} }_{i} \right) -b\, \dot{p}\, \hat{{\varvec{Z}}}+a\, \left| \hat{{\varvec{Z}}} \right| ^{m-1}\, \hat{{\varvec{Z}}} \right\} \end{aligned}$$

(25)

where b and m are further temperature-dependent material parameters. For small elastic deformations, the assumption of $\hat{\varvec{\varPhi }}=\hat{\varvec{\alpha }}$ is made [23, 26], which leads to a relationship similar to Armstrong-Frederick

$$\begin{aligned} \mathop {\hat{{\varvec{Z}}}}\limits ^{{\nabla }} =\, H_{i}\, \hat{{\varvec{D}}}_{p}-Q_{i}\, \dot{p}\, \hat{{\varvec{Z}}}+R\, \left| {\varvec{Z}} \right| ^{m-1}\, {\varvec{Z}} \end{aligned}$$

(26)

by combining the parameters for simplicity. The isotropic softening variables, being scalar quantities, are not affected by rotation. Therefore, they are described through the evolution of two independent variables that are related to their conjugate variables using Eq. (18) and $\Psi _{p}^{Iso}$ (Table 2):

$$\begin{aligned} r_{1}=\, -\beta \, \psi _{1},\, \, r_{2}=\beta \, \, \psi _{2}, \end{aligned}$$

(27)

with $\beta$ being a positive-valued, temperature-dependent material parameter. With the consideration that the isotropic softening variable is purely plasticity driven, the evolution equations of isotropic softening variables

$$\begin{aligned} \dot{\psi }_{1}= & {} -h\, \dot{p},\nonumber \\ \dot{\psi }_{2}= & {} c\left( \psi _{s}-\, \psi _{2} \right) \dot{p}-r_{\psi }\left| \psi _{2}-\psi _{r} \right| ^{m_{\psi }-1\, }\left( \psi _{2}-\psi _{r} \right) \end{aligned}$$

(28)

are taken from [24]. Here, $\psi _{s}=1-\, \psi _{s,\infty } \mathrm {\, }\left( 1-\exp \left( -c_{s}\, {{\epsilon } }_{eq,\, max}^{in}\, \right) \right)$ represents the increase of the cyclic softening capacity corresponding to increasing inelastic strain, with ${{\epsilon } }_{eq,\, max}^{in}$ being the highest value of the von Mises strain equivalent of $\hat{{\varvec{D}}}_{p}$ during the load history and $c,\, \, {m_{\psi },\, \, \psi }_{r},\, \, \psi _{s,\infty } \mathrm {\, and\, }c_{s}\mathrm {\, are}$ temperature-dependent material parameters. The irradiation behavior being isotropic is described by a scalar variable, defect density $N_{i}$ whose evolution for ‘i’ type of defects is given by the relation from [19]

$$\begin{aligned} \dot{N}_{i}\left( \phi ,\, p \right) =a_{i}\left( N_{s,i}-N_{i} \right) \dot{\phi }-b_{i}\left( N_{i}-N_{l,i} \right) \dot{p}-r_{N,i}\, N^{q_{N,i}}. \end{aligned}$$

(29)

Here, the first term describes the formation of defects, $N_{i}$ due to an irradiation dose, $\phi ,$ reaching a saturation value of $N_{s,i}$. The second term describes the removal of the portion of defects given by the variable $N_{l,i}$. The third term describes the healing of defects at high temperatures. The parameters $a_{i},\, b_{i},\, r_{N,i},\, q_{N,i},N_{s,i}$ are material and temperature dependent. In this work, the following relation is proposed between $N_{l,i}$ and p to better describe the influence of defect removal on stress–strain curve, with $N_{r,i}$ and $N_{e}$ being dimensionless, temperature and material-dependent parameters

$$\begin{aligned} {h_{N,i}^{2}\, N}_{l,i}=\,\left\langle \mathop {\hbox {max}}\limits _{0<\tau <t} \sqrt{h_{N,i}^{2}{\, N}_{i}\left( \tau \right) } \cdot \left( 1-N_{r,i}\left( 2-e^{-N_{e}\, p} \right) \right) \right\rangle ^{2}. \end{aligned}$$

(30)

Finally, the irradiation hardening $\sigma _{H}$ is obtained from $N_{i}$ for $n_{H}$ defect types contributing to the irradiation-induced hardening using the distributed barrier hardening (DBH) model [27]

$$\begin{aligned} \sigma _{H}=\sum \limits _i^{n_{H}} \sigma _{H,i} \, \, \, {\text {with}}\, \, \sigma _{H,i}=\, \sqrt{h_{N,i}^{2}{\, N}_{i}}. \end{aligned}$$

(31)

We assume that irradiation hardening is dominated by irradiation-induced dislocation loops. Therefore, $n_{H}\mathrm {=1}$ is a reasonable assumption [28, 29]. The conjugate ${\mathcalligra {n}}$ of N is got from Eq. (18) and $\Psi _{p}^{irr}$ in Table 2.

$$\begin{aligned} {\mathcalligra {n}}=\gamma \, N\, \end{aligned}$$

(32)

Here, ${\mathcalligra {n}}$ has $\textrm{MPa}$ as its unit, while N has a unit of ${\textrm{mm}}^{\mathrm {-3}}$.

The thermodynamic consistency conditions for the evolution of $\psi _{1},\, \psi _{2}$ and N are determined from the isotropic part of Clausius–Duhem entropy inequality (20), which is simplified using Eqs. (19), (22) and (23) into

$$\begin{aligned} \varSigma _{eq}\, \dot{p}+\varphi \, \dot{\phi }-r_{i}\, \dot{\psi }_{i}-{\mathcalligra {n}}\, \dot{N}\ge 0. \end{aligned}$$

(33)

From Eq. (21), we obtain the relation $\varSigma _{eq}\ge \, \sigma _{H}\, +\, k$. With Eqs. (13), (27), (28), (32) and (33), we deduce the following consistency criteria needed to meet the entropy inequality when $\dot{p}\ge 0\mathrm {\, }$and $\dot{\phi }\ge 0$. A detailed derivation is provided in the Appendix.

• 1.
  $$\begin{aligned} \psi _{1}=0,\, \psi _{2}=1\, {\text {at}}\, t=0 \end{aligned}$$

  (34)

  are the initial boundary conditions for the isotropic softening equations.

• 2.
  $$\begin{aligned} k-\, \, \beta \, \left[ \psi _{1}+\, c\, \psi _{2}\, \left( \psi _{s}-\, \psi _{2} \right) \right] \ge 0\, \rightarrow \, c\ge \frac{k-\, \beta \, \psi _{1}}{\, \beta \, \psi _{2}\left( \psi _{s}-\, \psi _{2} \right) } \end{aligned}$$

  (35)

  here,$\mathrm {\, }\psi _{s}\le \, \psi _{2}$ and this relation determines the condition to be fulfilled by the parameter c must fulfill.

• 3.
  $$\begin{aligned} \gamma _{r}\, \beta \, \psi _{2}\, \vert \psi _{2}-\psi _{r}\left. \right| ^{m_{\psi }}\ge 0\, \Rightarrow \, \gamma _{r}>0,\, \psi _{2}>\psi _{r} \quad \quad \because \, \, \, \, \, \, \, \beta>0,\, \psi _{2}>0 \end{aligned}$$

  (36)

suggests that $\psi _{2}$ lies in the range $\left[ \psi _{r},1 \right]$ ($\psi _{r}\le \psi _{2}\le 1).$

For the irradiated case, in addition to the above-discussed conditions we have the following conditions to fulfill:

• 1.
  $$\begin{aligned}{} & {} \frac{d\, \sigma _{H}}{d\, \phi }-\gamma \, N\cdot a\, \left( N_{s}-N \right) \ge 0\, \, \, \rightarrow \, \, \, a\left( N_{s}-N \right) \left[ \frac{h}{2\sqrt{N} }-\gamma \, N \right] \ge 0\, \end{aligned}$$

  (37)
  $$\begin{aligned}{} & {} \quad a\ge 0,\, \, N_{s}\ge N,\, \, \gamma \le \frac{h}{2\, N\, \sqrt{N} }\, \end{aligned}$$

  (38)

  where $a\ge 0$ is necessary for defect formation. Since $\gamma$ is positive, the condition for the selection of the $\gamma$ is

• 2.
  $$\begin{aligned}{} & {} 0\le \gamma \le \frac{h}{2\, N\, \sqrt{N} }. \end{aligned}$$

  (39)
  $$\begin{aligned}{} & {} \sigma _{H}+\gamma \, N\cdot \, b\, \left( N-N_{l}\, \right) \ge 0\, \, \rightarrow \, \, \, b\ge -\frac{\, \sigma _{H}}{\gamma \, N\, \left( N-N_{l} \right) } \end{aligned}$$

  (40)

  Since only available defects can be cleared, $N_{l}\le N$. Experiments show that $N_{l}$ is always lower than the maximum possible hardening; thus, $N_{l}<N_{s}$. As a result, the right-hand side of the equation will always be less than 0. Therefore, the required condition for the parameter is simply $b\ge 0$.

• 3.
  $$\begin{aligned} \gamma \, r_{N}N^{q_{N}+1}\ge 0\, \, \, \end{aligned}$$

  (41)

Since $\gamma >0,\, N\ge 0$, we obtain the condition that $r_{N}$ must always be positive, i.e., $r_{N}\ge 0$.

For the model to be implemented in FE programs, the evolution equations of tensor variables are transformed to $R_{a}$ using Eq. (3), while the equations of scalar variables remain unaltered. The Hookean hyperelasticity equation, $\mathop {\hat{{\varvec{T}}}}\limits ^\varDelta =\mathbb {C:}\mathop {\hat{\varvec{\varGamma }}_{e}}\limits ^{\varDelta }$, is used to describe stress evolution and, on transformation to $R_{a}$, becomes

$$\begin{aligned} \mathop {{\varvec{S}}}\limits ^\nabla \, =\hat{\varvec{R}}_{e}\mathop {{\hat{{\varvec{T}}}}}\limits ^\nabla \hat{\varvec{R}}_{e}^{{T}}= \frac{\partial \mathbb {C}}{\partial T} \dot{T}\, \text {:}\, {\varvec{A}}_{{\varvec{e}}}+\, \, \mathbb {C}\text {:}\left[ {\varvec{D}}-{\varvec{D}}_{{\varvec{p}}} \right] \end{aligned}$$

(42)

where $\mathbb {C}$ is the elasticity tensor. For the isothermal process, the first term on the right side reduces to zero. The transformation of Eq. (26) provides the evolution equation for the kinematic hardening variable

$$\begin{aligned} \mathop {\textbf{Z}}\limits ^\nabla \, =\hat{\varvec{R}}_{e}\mathop {\hat{{\varvec{Z}}}}\limits ^\nabla \hat{\varvec{R}}_e^{{T}}=\, H_{i}\, {\varvec{D}}_{p}-Q_{i}\, \dot{p}\, {\varvec{Z}}+R\, \left| {\varvec{Z}} \right| ^{m-1}\, {\varvec{Z}}. \end{aligned}$$

(43)

To integrate the model, a suitable integration algorithm is adopted and discussed in the next section.

The damage in irradiated Eurofer97 under monotonic loading is primarily due to accelerated void nucleation and coalescence. To describe this, a ductile damage law is the most suitable. Popular models [30, 31] used to model void coalescence are defined for rate-independent plasticity, which is in conjugation with damage-induced softening and highly mesh dependent. The proposal of a ductile damage law is targeted in future works, as an extension to the viscoplastic framework. Therefore, this paper does not examine thermodynamic relations for ductile damage.

3.3 Corotational integration algorithm

The stress response of deformation models must be incrementally objective, i.e., independent of finite-volume rotations and only contain material response. Therefore, a corotational integration algorithm is required to handle these finite rotations. The model is implemented in ABAQUS using the user subroutine UMAT [32]. Under the NLGEOM option, ABAQUS uses the corotational algorithm of Hughes and Winget [33] to determine the relative rotation, ${\varvec{Q}}_{HW}$, and the strain increment, $\varDelta \varvec{\epsilon }_{HW}= \varDelta t\, {\varvec{D}}_{HW}$, based on the constant deformation rate, ${\varvec{D}}_{HW}$, between the local material configurations at $t_{0}$, the beginning and $t_{1}$, the end of the increment. The inputs provided by ABAQUS to UMAT include the Cauchy stress, ${{\varvec{S}}}_{0}^{\mathbf {*}} ={\varvec{Q}}_{HW}\, {\varvec{S}}_{0}\, {\varvec{Q}}_{HW}^{T}$, and the user state variables (${\varvec{A}}_{p}{\varvec{,\, }}{\varvec{\Omega }}_{i},p,\, \Psi ,\, {N})$ at $t_{0}$. The user is required to provide ${{\varvec{S}}}_{1}^{\mathbf {*}}={{\varvec{S}}}_{0}^{\mathbf {*}}+\mathbb {C}{{\,:}} {\varDelta } {\varvec{\epsilon }}_{HW,\, e}$ at $t_{1}$. Along with ${\varvec{S}}_{1}^{*}$, the UMAT must return a consistent tangent modulus, ${\mathbb {C}}_{CTM}$ for the faster convergence of Newton iterations in the global system. Since the Hughes and Winget algorithm requires very small time increments to produce reliable results, this work employs an alternative integration algorithm developed by [34] that offers a stable and efficient computation of stress and state variables by proposing bar transformation of tensor variables and takes advantage of the available ${\varvec{Q}}_{HW}$ and $\varDelta \varvec{\epsilon }_{HW}$. The Runge–Kutta explicit integration scheme of order 5 is used to obtain results closest to the analytical solution within a specified tolerance.

First, the Oldroyd derivatives from Eqs. (42) and (43) are reformulated using Eqs. (2) and (6) to obtain the objective Jaumann derivatives as used by ABAQUS

$$\begin{aligned}{} & {} \mathop {{\varvec{S}}}\limits ^\circ ={{\varvec{DS}}}+{{\varvec{SD}}}+\mathbb {C}\text {:} \left[ {\varvec{D}}-{\varvec{D}}_{p} \right] \, \end{aligned}$$

(44)

$$\begin{aligned}{} & {} \quad \mathop {{\varvec{Z}}}\limits ^\circ = {{\varvec{DZ}}}+ {{\varvec{ZD}}} +H_{i}\, {\varvec{D}}_{p}-Q_{i}\, \dot{p}\, {\varvec{Z}}+R\, \left| {\varvec{Z}} \right| ^{m-1}\, {\varvec{Z}}. \end{aligned}$$

(45)

Next, the following algorithm is used to integrate the developed model:

• 1. Initial stress$\mathrm {\, }{\varvec{S}}_{0}$ is obtained from ${\varvec{S}}_{0}^{\mathbf {*}}\,$and ${\varvec{Q}}_{1}$, ${\varvec{S}}_{0}={\varvec{Q}}_{1}^{T}\, {\varvec{S}}_{0}^{\mathbf {*}}\, {\varvec{Q}}_{1}$.

• 2. Deformation rate, ${\varvec{D}}$ is calculated using ${\varvec{D}} = {\varDelta \varvec{\epsilon }}_{HW}/\varDelta t$.

• 3. The trial value of stress is calculated, ${\varvec{S}}_{1,\, trial}={\varvec{S}}_{0}+ \mathbb {C}:{\varDelta \varvec{\epsilon }}_{HW}$

• 4. Integrate Eq. (29) to determine irradiation hardening $\sigma _{H}$.

• 5. The yield criterion, $F\, =\, \mathrm {\varSigma }_{eq}-k-\, \sigma _{H}\ge 0$, is used to check for material yield. In the absence of yield, ${\varvec{S}}_{1}={\varvec{S}}_{1,\, trial}$ is set and the state variables remain unchanged. The elastic tensor is passed as the $\mathbb {C}_{CTM}{\varvec{.}}$

• 6. If material yields, constitutive equations (23), (28), (29), (44) and (45) are integrated to obtain the results ${\varvec{S}}_{1},\, {\varvec{A}}_{p,1},\, {\Omega }_{i,1},p_{1},\, \mathrm {\Psi }_{\textrm{i,1}}\mathrm {,\, }N_{\textrm{1}}$ at $t_{1}$. The consistent tangent modulus, $\mathbb {C}_{CTM}{\varvec{=\, }}\partial \varvec{\sigma /}\partial \varvec{\epsilon }$, is computed using the constitutive equations from a numerical approximation procedure of [35].

• 7. The tensor variables are finally rotated forward to the configuration at $t_{1}$: ${\varvec{S}}_{1}^{\mathbf {*}}={\varvec{Q}}_{HW}\, {\varvec{S}}_{1}\, {\varvec{Q}}_{HW}^{T}$, ${\varvec{Z}}_{i,1}^{\mathbf {*}}={\varvec{Q}}_{HW}\, {\varvec{Z}}_{i,1}\, {\varvec{Q}}_{HW}^{T}, \quad {\varvec{A}}_{p,1}^{\mathbf {*}}={\varvec{Q}}_{HW}\, {\varvec{A}}_{p,1}\, {\varvec{Q}}_{HW}^{T}$. Using a simplified approach by [36], the consistent tangent modulus is also rotated to $t_{1}$: ${\mathbb {C}}_{CTM}^{'}={\varvec{Q}}_{HW}{\varvec{\, }}{\varvec{Q}}_{HW}{\varvec{\, }}\mathbb {C}_{CTM}{\varvec{\, }}{\varvec{Q}}_{HW}^{T}{\varvec{\, }}{\varvec{Q}}_{HW}^{T}$.

• 8. The results of stress ${\varvec{S}}_{1}^{\mathbf {*}}$, internal state variables$\, {\varvec{A}}_{p,1}{\varvec{,\, }}{{\varvec{\Omega }} }_{i,1},p_{1},\, \Psi _{i,1} \mathrm {,\, }{N}_{\textrm{1}}$ and tangent modulus ${\mathbb {C}}_{CTM}^{'}$ are returned to ABAQUS as outputs of UMAT.

Using this implementation, simulations are performed to emulate the results of tensile test experiments.

4 Model application to experiments

Tensile tests were performed on unirradiated and irradiated specimens to calibrate the model. Using the UMAT, the tensile tests are simulated to assess the model’s capability to describe the post-yield and post-necking behavior under large strains while also describing the reduction in the minimum cross-sectional area.

Table 3 List of tests including testing parameters: irradiation dose and temperature, test temperature, and displacement rate

4.1 Tensile test experiment

Camera-monitored tensile tests were performed at SCK CEN as a part of the collaboration under the M4F projects (European Horizon2020) on unirradiated and irradiated samples of Eurofer97 at 300 $^\circ$C with a strain rate of 0.2 mm/min to assess the influence of neutron irradiation on the material’s tensile properties. Table 3 lists the test parameters, which show that the test temperature is 300 $^\circ$C. Figure 2 shows the specimen’s geometry. Figure 3 shows the engineering stress–strain $\left( \sigma _{eng}\mathrm {-\, }{{\epsilon } }_{eng} \right)$ relationships determined from experimental load–displacement data. Compared to the unirradiated specimen, the irradiated specimen has a higher yield stress, a brief strain hardening phase and a steep softening post the ultimate tensile strength (UTS). Fracture strain decreases with irradiation dose, indicating embrittlement. Analyzing optical images recorded at discrete deformation levels yields the mean true stress–strain relationships for unirradiated and 0.65 dpa irradiated samples. The Bridgman correction [37] is applied to determine the true stress–strain. Figure 4 shows these stress–strain relationships. The irradiated specimen’s true fracture strain is only slightly lower than that of the unirradiated specimen, and its maximum Bridgman true stress is very close to that of the unirradiated specimen. The irradiated specimen’s hardening modulus $\left( \partial \sigma _{true}/\partial {{\epsilon } }_{true}\right)$, on the other hand, is significantly lower than that of the unirradiated specimen.

Fig. 2

Geometry of the tensile specimen

Fig. 3

Engineering stress–strain curve of unirradiated and irradiated samples tested at ${T}_{test} = {T}_{{irr}} =$ 300 $^\circ$C. The irradiation conditions are given in the figure legend

Fig. 4

True stress, Bridgman corrected true stress plotted against true strain, along with horizontal and vertical error bars for the a unirradiated sample and b irradiated sample

True stress–strain data is needed to calibrate the finite strain model. This is composed of the uni-axial theoretical true stress $\left( \sigma _{T}=\sigma _{eng}[1+ {{\epsilon } }_{eng}] \right)$ and theoretical true strain $\left( \varepsilon _{T}=\ln \left[ 1+ {{\epsilon } }_{eng} \right] \right)$ up until necking and Bridgman corrected true stress–strain data after necking. The model stress can be decomposed into $\sigma _{T}/\psi =k+\sigma _{H}+\sigma _{vis}+Z\textrm{,}$with $\sigma _{H}=0$ for an unirradiated specimen and $\sigma _{vis}=K\left( \dot{{\epsilon } }_{p} \right) ^{1/n}$, the viscous stress. Since Eurofer97 has low strain rate sensitivity, $\sigma _{vis}=$ 1 MPa is assumed for the unirradiated specimen. From the yield stress $\sigma _{y}$ of the unirradiated specimen, the initial yield stress k is calculated. As the strain hardening of the material is not affected by irradiation hardening [38], the parameters for kinematic hardening and isotropic softening are calibrated using data from the unirradiated specimen. Since the isotropic softening parameters are known [13], $\psi$ can be calculated for the given plastic strain, allowing $\sigma _{T}/\psi$ to be determined. Since the loading is uniaxial and monotonic, parameters $H_{i}$ and $Q_{i}$ are curve-fitted using the estimated Z and the integral form of the Armstrong-Frederick law

$$\begin{aligned} Z=\, \sum \limits _i^3 Z_{i} \, =\sum \limits _i^3 {\, \frac{H_{i}}{Q_{i}}\, \left( 1-e^{-Q_{i}\, p} \right) }. \end{aligned}$$

Since the model’s viscous effects regulate the development of the neck, the viscous parameters K and n must be computed iteratively for the model’s load–displacement result to match the experimental data. Due to the low strain-rate dependence of the material, the parameters must be kept to a low value [39].

Table 4 Irradiation hardening ($T_{irr}=$300 $^\circ$C) in irradiated specimen, measured at the end of the linear elastic regime (yield), and at 0.2% strain

Table 5 Material parameters of the finite strain model determined for Eurofer97 at 300 $^\circ$C

Next, the parameters related to the irradiation-induced material changes are determined. The irradiation hardening, $\sigma _{H}$ (shown in Table 4), was determined by comparing the $R_{0.2}$ of both unirradiated and irradiated samples. Based on the $\sigma _{H}$ values obtained in this study and the saturation value (554 MPa) at higher doses (>10 dpa) from [19, 40], the parameters a and $h\sqrt{N}_{s}$ are fitted using the relation $\sigma _{H}=\, h\sqrt{N}_{s} \, \left( \, 1-\exp \left( -a\, \phi \right) \right) ^{0.5}$ [19]. As the parameter b is known [19], $N_{r}$ is obtained iteratively until the UTS predicted by the model matches the experimental results. Since irradiated materials have lower strain rate sensitivity [41], K and n are determined again iteratively. Table 5 shows the list of model parameters at 300 $^\circ$C.

4.2 Numerical simulation of tensile test

To simulate the tensile tests, an axisymmetric model is chosen to take advantage of the specimen symmetry, as illustrated in Fig. 5. The solid axisymmetric element type CAX4 with 4 nodes is chosen. To limit element distortion, a fine mesh with an aspect ratio of 1:5 is utilized in locations where large deformation is expected. The material parameters listed in Table 5 are utilized to simulate the model. A velocity displacement of $\mathrm {3.33\times }{\textrm{10}}^{\mathrm {-3}}$ mm/s is applied on the top of the shoulder along the y-axis, as seen in Fig. 5. The $\sigma _{eng}$ is calculated from the load distribution on the loading edge, $L\mathrm {=}\sigma _{loading\mathrm {\, }edge}\mathrm {\times }A_{shoulder}$ and the cross-sectional area of the gauge $A_{\textrm{0}}$. The engineering strain, ${{\epsilon } }_{eng}$ is obtained from the change in gauge length $\mathrm {\varDelta }l_{\textrm{0}}$ and the initial gauge length $l_{\textrm{0}}$

$$\begin{aligned}{} & {} \sigma _{eng}=\frac{L}{A_{0}}\, ,\, \, {{\epsilon } }_{eng}=\frac{\varDelta l_{0}}{l_{0}}. \end{aligned}$$

(46)

The mean true stress–strain is calculated using the minimum cross-sectional area at the neck, $A_{min}$

$$\begin{aligned}{} & {} \quad \sigma _{true}=\, \frac{L}{A_{min}},\, \, {{\epsilon } }_{true}=\ln {\left( \frac{A_{0}}{{\, A}_{min}}\, \right) \, }\, \end{aligned}$$

(47)

Since the model cannot predict ductile damage, the results of the simulation up to the experimental fracture strain are extracted and presented in the following discussions.

Fig. 5

Axisymmetric model of one-fourth tensile specimen cross section with boundary conditions

Fig. 6

Simulation results of the unirradiated sample with tension applied along the Y-axis: a deformed and undeformed body, b axial displacement (mm), c radial displacement (mm) and d accumulated plastic strain

4.2.1 Unirradiated specimen

The simulated deformation of the unirradiated sample tested at 300 $^\circ$C is shown in Fig. 6, including axial displacement, radius reduction, and accumulated plastic strain. Figure 7 compares the simulation’s stress–strain relationship with experimental results. The post-yield strain hardening and flow softening after necking are well represented. The predicted minimum neck diameter agrees with experimental results, allowing an accurate prediction of mean true stress–strain. Figure 8a shows the evolution of kinematic hardening, ${\varvec{Z}}$, using three variables. While two variables are sufficient to describe deformation up to the UTS, a third variable is required to describe material hardening until failure, as shown in the mean true stress–strain plot. Figure 8b shows the evolution of isotropic softening. Since material softening is not detected during monotonous loading, the suitably calibrated kinematic hardening compensates for the isotropic softening.

Fig. 7

Comparison of simulated tensile test results with experiment results of unirradiated sample at 300 $^\circ$C: a engineering stress–strain, b mean true stress–strain

Fig. 8

Evolution of a effective kinematic hardening with an inset plot of the constituting variables, and b isotropic variables and effective isotropic softening

4.2.2 Irradiated specimen

The specimen irradiation to 0.65 dpa and 1.18 dpa at 300 $^\circ$C followed by their tensile test at 300 $^\circ$C is simulated. Figure 9 shows the deformation of the specimen (0.65 dpa), illustrating the neck development at the middle of the gauge section. The same is observed from the simulation for 1.18 dpa. Figure 9c shows the distribution of defect density in the specimen after deformation, with no change in regions with no inelastic deformation and a minimum in regions with severe inelastic deformation, as seen in Fig. 9b. The simulated engineering stress–strain curves of both irradiated specimen are plotted in Figs. 10a and 11a, which compare extremely well to the experimental curves. The predicted mean true stress–strain is compared with experimental results for an irradiation of 0.65 dpa in Fig. 10b, confirming that the material hardening due to irradiation and the material softening due to defect annihilation post-yield are well represented by the model. The kinematic hardening and isotropic softening are the same as they are in the unirradiated case. The evolution of irradiation hardening concerning defect nucleation during irradiation and defect annihilation during plastic deformation is shown in Fig. 12. For $\phi =0.65$ dpa, a maximum hardening of 173 MPa is reached, and for $\phi =1.18$ dpa, a maximum hardening of 229 MPa is attained. During the initial part of the load history, inelastic deformation removes a portion of this hardening at an accelerated rate. When the annihilation reaches a threshold controlled by $N_{r}$, it slows down to a steady rate until the end of the load history. In other words, at the experimentally established fracture strain, residual irradiation-induced defects would still be present.

Fig. 9

Tensile test simulation of the irradiated sample (0.65 dpa, ${T}_{irr}={T}_{test}=$ 300 $^\circ$C) loaded along y-axis: a gauge radius reduction (mm), b accumulated plastic strain, c defect density distribution ($\mathrm {{{MP}}}\mathrm {{a}}^{\mathrm {{2}}}{)}$ and d) irradiation hardening distribution ($\mathrm {{MPa}}{\varvec{)}}$

Fig. 10

Tensile test simulation results of specimen irradiated to 0.65 dpa at 300 $^\circ$C, a engineering stress–strain, b mean true stress–strain

Fig. 11

Tensile test simulation results of specimen irradiated to 1.18 dpa at 300 $^\circ$C: a engineering stress–strain, b mean true stress–strain

Fig. 12

Evolution of irradiation hardening with a defect nucleation during irradiation compared with experimental results and b defect annihilation during plastic deformation

Fig. 13

Comparison of stress with respect to plastic strain at each integration point in unirradiated and irradiated materials

Fig. 14

Development of mean true strain with respect to engineering strain

5 Discussion

Experiment and simulation results are used to draw conclusions about Eurofer97’s behavior in unirradiated and irradiated states, with a focus on its deformation mechanism.

True stress–strain: Fig. 13 summarizes the true stress–strain predicted by the model at integration points from the center of the neck where triaxiality and plastic deformation are highest, demonstrating the influence of irradiation defect removal on hardening modulus under inelastic deformation. The small hill-like response near the UTS is due to the combination of strain hardening and irradiation hardening reduction. Experiments [38, 42] and simulations of micromechanical models [9, 10, 43] confirm these phenomena in irradiated specimens.

Minimum diameter prediction: The predicted mean true stress–strain plots compared against experimental results in Figs. 7, 10 and 11 indicate that reduction in gauge radius is predicted accurately. Figure 14 compares true strain to engineering strain to examine neck development in unirradiated and irradiated specimens. The absence of runaway behavior indicates that neck development is diffused in nature. In irradiated samples, the neck develops significantly earlier, and the rate of gauge diameter reduction increases proportionally with irradiation dose, resulting in higher local strain rates at smaller total elongations.

Triaxiality: Due to a concentration of inelastic deformation and defect removal, triaxiality is greatest near the neck. Figure 15 demonstrates that the presence of irradiation defects has a substantial impact on the triaxiality evolution in the neck. Future research will investigate the role of irradiation-induced material changes in reducing material ductility, given that ductile damage is among others triaxiality-driven.

Fig. 15

Triaxiality development at the necking region with defect removal under inelastic deformation with respect to a accumulated plastic strain and b irradiation hardening

Irradiation hardening: The projected irradiation hardening prior to plastic deformation exhibits some scattering when compared to the experimental results (Fig. 12). Due to the complexity of the defect type and population formed at low doses, irradiation hardening is extremely nonlinear at low irradiation doses, resulting in significant experimental scatter [44, 45]. As irradiation hardening saturates (~554 MPa) at higher doses (>10 dpa) [19, 40], the model predictions will be more consistent with experiments performed on specimens exposed to higher irradiation doses.

Strain rate sensitivity: Viscosity controls neck development and stress–strain curves. When the Considère instability criterion is met, local strain rate and viscous stress increase, affecting the progression of the instability and the formation of the neck. Overestimating viscous stress delays the neck development and overestimates the total elongation for the same minimum gauge diameter. Because Eurofer97 has a low strain rate sensitivity, [39] suggests maintaining low values for the viscous parameters n and K. Optimized parameters allow for the emergence of instabilities and the formation of a diffused and localized neck.

Corotational system: By using corotational formulation and integration, the engineering stress–strain curve relaxes significantly in the later load history when the neck elements undergo finite rotations.

Performance: The model was integrated using an explicit algorithm that limited time increments. Smaller time steps restrict plastic strain increments, ensuring system stability and convergence while achieving a solution close to the analytical solution.

Damage: The presented model is lacking the component that describes the ductile damage mechanism, which would aid in predicting the specimen’s fracture. In addition, the damage model will help us to determine the effect of irradiation on triaxiality and, consequently, the fracture prediction.

6 Conclusion

A finite strain framework is established in order to develop a model for describing the time-dependent plasticity and damage in the irradiated material.

The proposed inclusion of irradiation-related power term in the law of energy conservation, which is based on the introduction of instantaneous defect density as an internal state variable of the material’s free energy, accounts for the energy exchange arising from the nucleation and annihilation of irradiation defects. The resulting equation of the modified-entropy inequality is used to derive the finite strain formulation of an existing deformation model. The conditions necessary to fulfill the thermodynamic consistency criteria are presented. The model is implemented in the ABAQUS finite element code to simulate tensile tests on unirradiated and irradiated specimens at elevated temperatures. The post-yield and post-necking behavior is well captured by the engineering stress–strain and true stress–strain curves produced from the simulation results. The role of viscosity effects on diffuse and localized neck formation is found to be a crucial component. Since the model does not yet include the part describing ductile damage, it shall be considered in subsequent works by coupling the model with a suitable damage model for void nucleation and coalescence caused by viscoplastic deformation.

Appendix

1.1 Thermodynamics

The balance of entropy in its local form is:

$$\begin{aligned} \dot{\eta }=\frac{s_{total}}{\theta }-\frac{1}{\rho }\, \nabla \cdot \frac{{\varvec{q}}}{T}+\gamma .\, \, \end{aligned}$$

(48)

Using the specific entropy, $\eta$, the absolute temperature $T\textrm{,}$ and the specific free energy relation $\Psi =e-T\, \eta$, modified first law Eq. (10) can be rewritten for the specific rate of entropy production, $\gamma$ as

$$\begin{aligned} \theta \, \gamma =\, -\dot{\Psi }-\dot{T}\eta +\frac{1}{\rho }\, \left( {\varvec{T}}\cdot \dot{{\varvec{E}}}+\xi \, \dot{\phi }\, \right) -\, \frac{{\varvec{q}}}{\rho \, T}\cdot \, \nabla T. \end{aligned}$$

(49)

With the positive entropy production according to the second law of thermodynamics, i.e., $\gamma \ge 0$, we have the following modified form of Clausius–Duhem entropy inequality

$$\begin{aligned} -\, \dot{\Psi }\, -\dot{T}\, \eta \, +\frac{1}{\rho }\, \left( \, {\varvec{T}}\cdot \dot{{\varvec{E}}}+\xi \, \dot{\phi } \right) \, -\, \frac{{\varvec{q}}}{\rho \, T}\cdot \nabla T\, \ge \, 0. \end{aligned}$$

(50)

1.2 Kinematic hardening

$$\begin{aligned}{} & {} \Psi _{p}^{Kin}\left( t \right) =\frac{A\left( p \right) }{2\rho _{R}}\, \hat{\varvec{\alpha }}_{i}\cdot \hat{\varvec{\alpha }}_{i} \end{aligned}$$

(51)

$$\begin{aligned}{} & {} \quad \hat{{\varvec{Z}}}_{i}=\, \rho _{R}\frac{\partial \Psi }{\partial \hat{\alpha }_{i}\, }=A\left( p \right) \, \hat{\varvec{\alpha }}_{i} \end{aligned}$$

(52)

Kinematic part of entropy inequality (20) using the Oldroyd derivative from Eq. (7):

$$\begin{aligned}{} & {} \left( \hat{{\varvec{D}}}_{p}-\frac{\left( \hat{{\varvec{Z}}} \right) ^ {\mathrm {\nabla }}}{A\left( p \right) }+2\, \hat{{\varvec{D}}}_{p}\, \left( \hat{{\varvec{\varPhi }} }\mathbf {-\, }\hat{{\varvec{\alpha }} }_{i} \right) \right) \cdot \hat{{\varvec{Z}}}_{i}\ge 0 \end{aligned}$$

(53)

$$\begin{aligned}{} & {} \quad \hat{{\varvec{D}}}_{p}\mathbf {-}\frac{\left( \hat{{\varvec{Z}}} \right) ^{\mathrm {\mathbf {\nabla }}}}{A\left( p \right) }\mathbf {+}2{\, }\hat{{\varvec{D}}}_{p}\,\left( \hat{{\varvec{\varPhi }}} -\,\hat{{\varvec{\alpha }}}_{i} \right) =\, \kappa _{p}\, \hat{{\varvec{Z}}}_{i}\, \end{aligned}$$

(54)

Choosing $A\left( p \right) =C\, a\, f\left( p \right)$ and $\kappa _{p}=\frac{b}{a}\dot{p}+\left| \hat{{\varvec{Z}}} \right| ^{m-1}\ge 0$.

$$\begin{aligned} \left( \hat{{\varvec{Z}}} \right) ^{\mathrm {\nabla }}=C\, f\left( p \right) \left\{ a\hat{{\varvec{D}}}_{p}+2\, \hat{{\varvec{D}}}_{p}\, \left( \hat{{\varvec{\varPhi }} }-\, \hat{{\varvec{\alpha }} }_{i} \right) -b\, \dot{p}\, \hat{{\varvec{Z}}}+a\, \left| \hat{{\varvec{Z}}} \right| ^{m-1}\, \hat{{\varvec{Z}}} \right\} \end{aligned}$$

(55)

1.3 Isotropic softening and irradiation defect density

$$\begin{aligned} \Psi ^{iso}\left( t \right) \, \, =\, \, \frac{1}{2\rho _{R}\, }\, \beta \, \left( \psi _{2}^{2}-\, \psi _{1}^{2} \right) \, \, +\, \, \frac{1}{2\rho }\cdot \gamma \, N^{2} \end{aligned}$$

(56)

Here $\beta$ is a positive definite material parameter.

$$\begin{aligned}&\displaystyle r_{1}=\, \rho _{R}\frac{\partial \Psi ^{iso}}{\partial \psi _{1}\, }=-\, \beta \, \psi _{1},\, \, r_{2}=\, \rho _{R}\frac{\partial \Psi ^{iso}}{\partial \psi _{2}\, }=\, \beta \, \, \psi _{2} \end{aligned}$$

(57)

$$\begin{aligned}&\displaystyle {\mathcalligra {n}} =\, \rho _{R}\frac{\partial \Psi ^{iso}}{\partial N}\, =\gamma \, N \end{aligned}$$

(58)

Substituting the above in the isotropic part of the entropy inequality from Eq. (20)

$$\begin{aligned}{} & {} \dot{p}\, \left( k+\, \sigma _{H} \right) +\frac{d\sigma _{H}}{d\phi }\, \dot{\phi }-\beta \, \psi _{1}\, h\, \dot{p}-\, \beta \, \psi _{2}\left[ c\, \left( \psi _{s}-\psi _{2} \right) \dot{p}-r_{\psi }\, \right. \nonumber \\{} & {} \quad \left. \left( \psi _{2}-\psi _{r} \right) \right] -\gamma \, N\, \left[ a\, \left( N_{s}-N\, \right) \dot{\phi }-b\, \left( N-N_{l} \right) \, \dot{p}-r_{N}N^{q_{N}}\, \right] \ge 0\, \end{aligned}$$

(59)

$$\begin{aligned}{} & {} \dot{p}\, \left[ k-\, \psi _{1}\beta \, h-c\, \psi _{2}\, \beta \, \left( \psi _{s}-\psi _{2} \right) \right] +\, \, \beta \, \psi _{2}\, r_{\psi }\left( \psi _{2}-\psi _{r} \right) +\dot{p}\left[ \sigma _{H}+\gamma \, N\cdot \, b\, \left( N-N_{l}\, \right) \right] \nonumber \\{} & {} \quad +\dot{\phi }\, \left( \frac{d\sigma _{H}}{d\phi }-\gamma \, N\cdot a\, \left( N_{s}-N\, \right) \right) +\gamma \, r_{N}N^{q_{N}+1}\ge 0\, \, \, \end{aligned}$$

(60)

The term $d\sigma _{H}/d\phi$ is determined as follows

$$\begin{aligned} \sigma _{H}=h\sqrt{N} \rightarrow \frac{d\sigma _{H}}{d\phi }=\frac{h}{2\sqrt{N} }\, \, \frac{dN}{d\phi } \end{aligned}$$

(61)

$$\begin{aligned} N\left( \phi ,\, p \right) \rightarrow \, \dot{N}=\, \frac{dN}{d\phi }\dot{\phi }\, +\, \, \frac{dN}{dp}\dot{p} \end{aligned}$$

(62)

Comparing with $\dot{N}$ evolution equation (29), we have

$$\begin{aligned}{} & {} \frac{dN}{d\phi }=a\left( N_{s}-N \right) ,\, \, \frac{dN}{dp} =-b\left( N-N_{l} \right) \end{aligned}$$

(63)

$$\begin{aligned}{} & {} \quad \frac{d\sigma _{H}}{d\phi }=\frac{h}{2\sqrt{N} }\, a\left( N_{s}-N \right) \end{aligned}$$

(64)