# NN-EVP: A physics informed neural network-based elasto-viscoplastic
framework for predictions of grain size-aware flow response under large
deformations

Adnan Eghtesad<EMAIL_ADDRESS>Jan Niklas Fuhg  Nikolaos Bouklas  Sibley
School of Mechanical and Aerospace Engineering, Cornell University, NY 14853,
USA  Center for Applied Mathematics, Cornell University, NY 14853, USA

###### Abstract

We propose a physics informed, neural network-based elasto-viscoplasticity
(NN-EVP) constitutive modeling framework for predicting the flow response in
metals as a function of underlying grain size. The developed NN-EVP algorithm
is based on input convex neural networks as a means to strictly enforce
thermodynamic consistency, while allowing high expressivity towards model
discovery from limited data. It utilizes state-of-the-art machine learning
tools within PyTorch’s high-performance library providing a flexible tool for
data-driven, automated constitutive modeling. To test the performance of the
framework, we generate synthetic stress-strain curves using a power law-based
model with phenomenological hardening at small strains and test the trained
model for strain amplitudes beyond the training data. Next, experimentally
measured flow responses obtained from uniaxial deformations are used to train
the framework under large plastic deformations. Ultimately, the Hall-Petch
relationship corresponding to grain size strengthening is discovered by
training flow response as a function of grain size, also leading to efficient
extrapolation. The present work demonstrates a successful integration of
neural networks into elasto-viscoplastic constitutive laws, providing a robust
automated framework for constitutive model discovery that can efficiently
generalize, while also providing insights into predictions of flow response
and grain size-property relationships in metals and metallic alloys under
large plastic deformations.

###### keywords:

Viscoplasticity , Flow response , Grain size , Hall-Petch , Machine learning ,
Neural networks

##  1 Introduction

Understanding the flow response and viscoplastic behavior of metals
incorporating strain rate sensitivity effects, is essential for the
predictions of the mechanical behavior, improving the material performance,
and designing novel alloys, critical towards revolutionizing the objectives of
the aerospace and automotive manufacturing industries. Owing to their
polycrystalline nature, the anisotropic flow response of metals under large
deformations heavily depends on the average grain size that forms the
underlying microstructure. Grain size strengthening, referred to as the Hall-
Petch effect [1, 2, 3], is associated with the grain boundary resistance to
the crystallographic deformation mechanisms and dislocation motions [4, 5].
Authors in [6] implemented a crystal plasticity model that accounts for grain
size effects and slip system interactions on the deformation of austenitic
stainless steels. The crystal plasticity finite element method (CPFEM) is used
to study the grain size and morphology effects on yield strength [7]. In
addition to crystal plasticity modeling, atomistic simulations are well
utilized to study the Hall-Petch relationships in advanced alloys [8, 9].

High-performance computing (HPC) has improved the efficiency of numerical
techniques over the last few decades. However, despite acceleration of HPC
simulations, predicting multi-scale material deformations is still a time-
consuming task and limited by parallel scalability and hardware performance.
To address this issue, recent research has focused on the integration of
genetic algorithms (GA), machine learning (ML), and deep learning (DL) to
facilitate automated constitutive modeling which has the potential to expedite
multi-scale simulations [10, 11, 12]. In particular, they have been used to
model hyperelasticity [13, 14, 15, 16, 17, 18], viscoelasticity [19, 20, 21,
22] and multiphysics problems [23, 24].

Lately, neural networks (NN) and convolutional neural networks (CNN) have also
been utilized for a variety of applications in modeling viscoplastic
deformations in macro and micro scales [25]. A recurrent neural network (RNN)
model was proposed as a computationally-efficient surrogate for crystal
plasticity simulations [26]. A new NN-based crystal plasticity algorithm was
presented for FCC materials and its application to non-monotonic strain paths
[27] followed by a CNN-based CPFEM model to predict the localized viscoplastic
deformation in aluminum alloys [28]. Some recent work has implemented NN and
CNN to predict yield surfaces from microstructural images and crystal
plasticity simulations [29, 30]. A machine learning-enabled crystal plasticity
model with dislocation density hardening was developed to identify stress and
strain localizations under large viscoplastic deformations [31]. An input
convex NN (ICNN) framework was presented for the prediction of texture-
dependent macroscopic yield functions from crystal plasticity simulations
[32]. On the continuum level, NN algorithms have been implemented in finite
element (FE) solvers to replace classical history-dependent constitutive
models to obtain nonlinear structural responses [33, 34, 35, 36, 37].

While the literature addresses a wide range of machine and deep learning
applications related to elastoplasticity and viscoplasticity in the context of
big data, only a few rely on the more realistic case of limited-data
availability of macroscopic observations from measured tensile tests.
Furthermore, studies that focus on viscoplasticity in the low- and limited-
data regimes, only utilize neural networks for parameter estimation of
established phenomenological constitutive models rather than model discovery.
The latter can focus on establishing automated frameworks that remove the need
for a particular phenomenological parametrization. To address this
shortcoming, and allow robust performance in low- and limited-data regimes,
physics-based ML algorithms have been proposed aiming to directly enforce
physical laws, thermodynamic principles and also established domain knowledge
[38]. Motivated by the earlier work of [39], recent work in [40], enables
modular ML-based elastoplasticity in the context of limited data using
thermodynamics-aware and mechanistically informed neural networks offering a
hybrid framework with each component of the model being selective as classical
phenomenological or a data-driven model depending on the data availability,
allowing for trustworthy prediction and generalization.

The present work builds upon the work in [40] and proposes a novel NN-based
framework for modeling the elasto-viscoplastic (NN-EVP) response of metals as
a function of grain size where strain rate sensitivity effects are taken into
consideration. The proposed algorithm and underlying constitutive model
consider the coupled elasto-viscoplastic response in contrast to the studies
where elastic and plastic regimes are treated separately. In addition, for
consistency with the thermodynamic laws and mechanistic assumptions, neural
network architectures are designed to enforce monotonicity and convexity
requirements. The framework allows for the discovery of laws describing the
rate-sensitive viscoplastic flow response of metals and alloys, enabling a
predictive infrastructure for the flow response as a function of strain
amplitude and grain size. The present paper is structured as follows. Section
2 starts with a brief review of the thermodynamics-based modular formulation,
based on dual potentials, as well as particular functional forms, resulting
from mechanistic assumptions, that describe the rate-sensitive viscoplastic
flow response in metals. Afterward, the key aspects pertaining to the novel
implementations of the NN-EVP framework are discussed. Section 3 describes the
applications of the NN-EVP tool in modeling large viscoplastic deformations
with isotropic hardening for both synthetic and experimental data as a
function of grain size. A summary of the findings and key contributions
concludes the work.

##  2 Methods

###  2.1 Elasto-viscoplasticity constitutive formulation based on dual
potential

Under the assumption of small strains, elasto-viscoplasticity can be modeled
by decomposing the strain into its elastic and viscoplastic parts:

$\bm{\epsilon}=\bm{\epsilon}^{e}+\bm{\epsilon}^{vp},\\\ $ (1)

where $\bm{\epsilon}$, $\bm{\epsilon}^{e}$, and $\bm{\epsilon}^{vp}$ are the
total strain, elastic strain, and viscoplastic strain, respectively.

In a similar fashion, the specific free energy $\Psi$ can be decomposed into
elastic and viscoplastic contributions, denoted by $\Psi^{e}$ and $\Psi^{vp}$
respectively, as follows:

$\Psi=\Psi^{e}(\bm{\epsilon}^{e})+\Psi^{vp}(r,\bm{\alpha},T),$ (2)

where $r$ and $\bm{\alpha}$ are the thermodynamic variables related to
isotropic and kinematic hardening, and $T$ is the temperature. Following the
Coleman-Noll procedure [41] the Cauchy stress $\sigma$ and thermodynamic
forces $R$ and $X$ can then be defined as:

$\bm{\sigma}=\dfrac{\partial\Psi}{\partial{\bm{\epsilon}^{e}}}=-\dfrac{\partial\Psi}{\partial{\bm{\epsilon}^{vp}}},$
(3) $R=\dfrac{\partial\Psi}{\partial{r}}$ (4)
$\mathbf{X}=\dfrac{\partial\Psi}{\partial{\bm{\alpha}}}.$ (5)

Under isothermal conditions, the intrinsic potential is written as:

$\Psi_{I}=\Psi_{I}^{e}(\bm{\epsilon}^{e})+\Psi_{I}^{vp}(\bm{\alpha},r).$ (6)

The generalization of the concept of equipotential surfaces in stress space,
representing similar dissipation (i.e., strain rate) levels leads to the
definition of a dual potential in the form of [42]:

$\varphi^{*}=\varphi^{*}(\bm{\sigma},R,\mathbf{X};r,\bm{\alpha}).$ (7)

The dual dissipation potential $\varphi^{*}$ satisfies the following
conditions:

  * •

$\varphi^{*}$ is convex w.r.t to all its variables,

  * •

$\varphi^{*}$ is always positive: $\varphi^{*}>0$,

  * •

$\varphi^{*}$ includes the origin: $\varphi^{*}(0,0,0,\bm{\alpha},r)=0$.

Following the normality law of generalized standard materials the plastic
strain rate $\dot{\epsilon^{vp}}$ and the hardening variables $r$ and $\alpha$
can then be obtained as [42]:

${\dot{\epsilon}^{vp}}=\dfrac{\partial\varphi^{*}}{\partial\bm{\sigma}},$ (8)
$\dot{r}=-\dfrac{\partial\varphi^{*}}{\partial{R}},$ (9)
$\dot{\bm{\alpha}}=-\dfrac{\partial\varphi^{*}}{\partial\mathbf{X}},$ (10)

which leads to the second law of thermodynamics as follows:

$\Psi_{I}=\bm{\sigma}:\dfrac{\partial\varphi^{*}}{\partial\bm{\sigma}}+\mathbf{X}:\dfrac{\partial\varphi^{*}}{\partial\mathbf{X}}+R\dfrac{\partial\varphi^{*}}{\partial{R}}\geq\varphi^{*}\geq
0.$ (11)

####  2.1.1 Decomposition of dual potential

Under consideration of the recovery effects caused by dislocation annihilation
we can formulate the dual potential $\varphi^{*}$ to also be decomposed into
viscoplastic hardening and recovery potentials as follows [42]:

$\varphi^{*}=\Omega_{vp}+\Omega_{r},$ (12)

where

$\Omega_{vp}=\Omega_{vp}(\bm{\sigma},R,\mathbf{X};r,\bm{\alpha}),$ (13)
$\Omega_{r}=\Omega_{r}(R,\mathbf{X};r,\bm{\alpha}).$ (14)

In the absence of recovery effects, the dual potential can be considered to be
equal to the hardening potential $\Omega_{p}$ [42]:

$\varphi^{*}=\Omega_{vp}=\varphi^{*}(\bm{\sigma},R,\mathbf{X};r,\bm{\alpha}).$
(15)

####  2.1.2 Particular functional forms and power-law

Assuming linear elastic behavior for the elastic part of the free energy
function we can write

$\Psi_{I}^{e}(\bm{\epsilon}^{e})=\frac{1}{2}\bm{\epsilon}^{e}:\mathbb{C}:\bm{\epsilon}^{e}$
(16)

where $\mathbb{C}$ is the fourth-order elastic stiffness tensor. This yields
the Cauchy stress through Hook’s law as follows

${\bm{\sigma}}=\mathbb{C}:\bm{\epsilon}^{e}.$ (17)

For isotropic materials, the stiffness tensor is uniquely defined by the
Young’s modulus $E$ and Poisson’s ratio $\nu$

$\mathbb{C}_{ijkl}=\frac{E\nu}{(1+\nu)(1-2\nu)}\delta_{ij}\delta_{kl}+\frac{E}{2(1+\nu)}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$
(18)

where $\delta_{ij}$ denotes the Kronecker delta.

Under the assumption of isotropic and pressure-independent material response,
we can reformulate the dual potential as a function of only the second and
third stress invariants

$\varphi^{*}=\varphi^{*}(J_{2}(\bm{\sigma}-\mathbf{X}),J_{3}(\bm{\sigma}-\mathbf{X}),R,\mathbf{X};r,\bm{\alpha}),$
(19)

where

$J_{2}(\bm{\sigma}-\mathbf{X})=\sqrt{\dfrac{3}{2}}{\lVert{\bm{\sigma}^{{}^{\prime}}-\mathbf{X}^{{}^{\prime}}}\rVert},$
(20) $J_{3}(\bm{\sigma}-\mathbf{X})=\dfrac{1}{3}tr(\bm{\sigma}-\mathbf{X}),$
(21)

and

$\bm{\sigma}^{{}^{\prime}}=\bm{\sigma}-\dfrac{1}{3}tr(\bm{\sigma}),$ (22)
$\mathbf{X}^{{}^{\prime}}=\mathbf{X}-\dfrac{1}{3}tr(\mathbf{X}).$ (23)

Here, the apostrophe $\square^{{}^{\prime}}$ indicates the deviatoric
component of the tensor and $tr(\,)$ implies the trace of the argument. After
assuming that due to the incompressible behavior of metals, only the
deviatoric components of the stress contributes to the plastic deformation,
the dual potential reduces to [42]:

$\varphi^{*}=\varphi^{*}(J_{2}(\bm{\sigma}-\mathbf{X}),R,\mathbf{X};r,\bm{\alpha}),$
(24)

The fundamental flow rule in viscoplasticity theory involving the
thermodynamic force $R(r)$ is commonly written in the form of a power-law and
is desirable because it provides uniqueness of solution for the stress value
that accommodates an imposed strain rate [43]:

$\varphi^{*}=\dfrac{\dot{\epsilon}_{0}}{n+1}{\absolutevalue{\dfrac{J_{2}(\bm{\sigma}-\mathbf{X})}{R(r)}}^{n+1}}.$
(25)

Incorporating Eq. (25) and Eq. (20) into Eq. (8) we get that

$\dot{\epsilon^{vp}}=\dfrac{3}{2}\dot{\epsilon}_{0}{\absolutevalue{\dfrac{J_{2}(\bm{\sigma}-\mathbf{X})}{R(r)}}^{n}}\dfrac{\bm{\sigma}^{{}^{\prime}}-\mathbf{X}^{{}^{\prime}}}{\lVert{\bm{\sigma}^{{}^{\prime}}-\mathbf{X}^{{}^{\prime}}}\rVert}.$
(26)

The power-law embeds the introduced strain rate sensitivity and provides
uniqueness of solution for the threshold of equivalent von Mises stress
accommodating an imposed strain rate. In Eqs 24 and 25, $n$ is the rate
sensitivity parameter and $\dot{\epsilon}_{0}$ is a reference strain rate
usually chosen as the norm of applied strain rate tensor
${\lVert{\dot{\bm{\epsilon}}^{app}}\rVert}$. Following Eq. 9, the variable $r$
and its rate $\dot{r}$, implying the effective accumulated visco-plastic
strain and effective visco-plastic strain rate can be written as follows:

$r=\epsilon_{eff}^{vp}=\int_{0}^{t}\sqrt{\dfrac{2}{3}}\lVert{\dot{\bm{\epsilon}^{vp}}(\tau)\rVert}d\tau,$
(27)
$\dot{r}=\dot{\epsilon}_{eff}^{vp}=\sqrt{\dfrac{2}{3}}\lVert{\dot{\bm{\epsilon}^{vp}}\rVert}.$
(28)

####  2.1.3 Perfect viscoplasticity

In case of perfect viscoplasticity and no hardening effects, the yield surface
does not evolve as a function of accumulated viscoplastic strain, and the term
$R(r)$ reduces to an initial yield value $\sigma^{Y}$. The dual potential and
viscoplastic strain rate can then be written as:

$\varphi^{*}=\dfrac{\dot{\epsilon}_{0}}{n+1}{\absolutevalue{\dfrac{\sigma_{eq}}{{\sigma^{Y}}}}^{n+1}},$
(29)
$\dot{\bm{\epsilon}^{vp}}=\dfrac{3}{2}\dot{\epsilon}_{0}{\absolutevalue{\dfrac{\sigma_{eq}}{\sigma^{Y}}}^{n}}\dfrac{\bm{\sigma}^{{}^{\prime}}}{\sigma_{eq}},$
(30)

where
$\sigma_{eq}=\sqrt{\dfrac{3}{2}}{\lVert{\bm{\sigma}^{{}^{\prime}}}\rVert}$.
Depending on the material properties, the rate sensitivity parameter varies in
the range of $10\leq n\leq 400$. Figure 1 shows the effect of rate sensitivity
on perfect viscoplastic flow response of Cu under applied quasi-static strain
rate of 0.001$s^{-1}$ generated by the power law equation mentioned above.
Note that while higher values of $n$ represent the rate sensitivity of the
material more accurately, the computational cost and number of iterations
required to numerically solve Eq. (30) increases.

Figure 1: Effect of rate sensitivity on perfect viscoplastic deformation of Cu
generated using the power law equation.

####  2.1.4 Isotropic hardening

In the case of isotropic hardening, the yield function $R(r)$ starts with an
initial value $\sigma^{Y}$ ($R(r=0)=\sigma^{Y}$) and evolves as a function of
accumulated viscoplastic strain $r$. The dual potential and viscoplastic
strain rate then yield

$\varphi^{*}=\dfrac{\dot{\epsilon}_{0}}{n+1}{\absolutevalue{\dfrac{\sigma_{eq}}{R(r)}}^{n+1}},$
(31)
$\dot{\bm{\epsilon}^{vp}}=\dfrac{3}{2}\dot{\epsilon}_{0}{\absolutevalue{\dfrac{\sigma_{eq}}{R(r)}}^{n}}\dfrac{\bm{\sigma}^{{}^{\prime}}}{\sigma_{eq}}.$
(32)

Figure 2 shows the viscoplastic flow response in Cu generated by the power law
and as a function of applied strain rate and rate sensitivity parameter
$n=20$, demonstrating significant anisotropy in the flow response as a
function of imposed strain rate.

Figure 2: Effect of applied strain rate in viscoplastic flow response of Cu
generated by the power law with a rate sensitivity of $n$=20.

While the particular functional forms and the power law model described above
enable the constitutive formulation for viscoplastic flow response in metals,
predicting the stress-strain behavior using these models often require manual
calibration of a set of hardening parameters with experimental data obtained
from uniaxial deformation tests. While the automation of the calibration
procedures have been enabled by the existing ML and GA models, training these
models can be still a time-consuming task, especially when discrepancies of
the experimental response and the suggested power law model are present.
Therefore, we propose a data-driven elasto-viscoplasticity framework based on
neural networks, NN-EVP, that replaces the power law with a generic neural
networks based algorithm capable of predicting the flow response in metals and
alloys at large deformations in the context of limited data availability.

###  2.2 NN-EVP

In order to employ a physics-informed data-driven model for the elasto-
viscoplastic constitutive formulation described earlier, the dual potential
$\varphi^{*}$ is represented by the response of a neural
${\mathcal{NN}_{\phi^{*}}}$ whose (scalar) input is chosen based on which
hardening phenomena we aim to be captured. This is discussed in more detail in
the subsequent sections.

To remain consistent with thermodynamics laws, the following conditions are
(implicitly) enforced on the dual potential neural network
${\mathcal{NN}_{\phi^{*}}}$:

  * •

${\mathcal{NN}_{\phi^{*}}}$ is convex and monotonically increasing

  * •

${\mathcal{NN}_{\phi^{*}}}$ is always positive: ${\mathcal{NN}_{\phi^{*}}}$ >
0

  * •

${\mathcal{NN}_{\phi^{*}}}$ includes the origin:
${\mathcal{NN}_{\phi^{*}}(0)}=0$

If ${\mathcal{NN}_{\phi^{*}}}:\mathbb{R}\rightarrow\mathbb{R}$ is a
feedforward neural network with $L$ hidden layers, input ${x_{0}}$ and output
${y_{L}}$, the neural network can be written as follows:

$\displaystyle x_{0}$ $\displaystyle\in\mathbb{R}_{\geq 0},$ (33)
$\displaystyle{x}_{1}=\mathcal{F}_{1}\left(x_{0}{W}_{1}^{T}+{b}_{1}\right)$
$\displaystyle\in\mathbb{R}^{n^{1}},$
$\displaystyle{x}_{l}=\mathcal{F}_{l}\left({x}_{l-1}{W}_{l}^{T}+{b}_{l}\right)$
$\displaystyle\in\mathbb{R}^{n^{l}},\qquad l=1,\ldots,L-1$ $\displaystyle
x_{L}={x}_{L-1}{W}_{L}^{T}+{b}_{L},$ $\displaystyle\in\mathbb{R},$

where the weights ${W}_{l}\in\mathbb{R}^{n^{l}\times n^{l-1}}$ and biases
${b}_{l}\in\mathbb{R}^{n^{l}}$ define the set of trainable parameters and the
activation functions are denoted by $\mathcal{F}_{1}$. The neural network
${\mathcal{NN}_{\phi^{*}}}$ is positive, monotonically increasing, and input
convex when the following conditions are met [40]:

  * •

$x_{0}\geq 0,W_{l}\geq 0,b_{l}\geq 0$,

  * •

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}}_{l}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0},l=1,\ldots,L$,

  * •

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}}_{l}^{\prime}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0},l=1,\ldots,L$,

  * •

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}}_{l}^{\prime\prime}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0},l=1,\ldots,L$.

In order to satisfy the conditions above, we choose a parameterized and
adaptive/scalable Softplus activation function [44, 45] as follows:

$\mathcal{F}_{l}^{\mathcal{NN}_{\phi^{*}}}(x)=\frac{1}{\beta}\log(1+e^{x})-\mathcal{F}_{l}^{\mathcal{NN}_{\phi^{*}}}(0),$
(34)

where $\beta>0$ is a training parameter in addition to the weights and biases
of the neural network, providing a more generic form of the activation
function and thus more flexibility in training. Note that $\beta$ is trainable
for each hidden layer of ${\mathcal{NN}_{\phi^{*}}}$. The deduction of
$\mathcal{F}_{l}^{\mathcal{NN}_{\phi^{*}}}(0)$ from the equation above is to
satisfy the condition ${\mathcal{NN}_{\phi^{*}}(0)}=0$. Once the dual
potential ${\mathcal{NN}_{\phi^{*}}}$ is established, it can be used to
replace the particular power law form introduced earlier in Eqs (29) and (31).
Equivalently, the viscoplastic strain rate ${\dot{\epsilon}^{vp}}$ can be
obtained by taking the derivative of the output of ${\mathcal{NN}_{\phi^{*}}}$
with respect to Cauchy stress:

$\bm{\dot{\epsilon}^{vp}}=\dfrac{\partial{\mathcal{NN}_{\phi^{*}}}}{\partial{\bm{\sigma}}}.$
(35)

####  2.2.1 NN-EVP for perfect viscoplasticity

In the case of perfect viscoplasticity, the yield function $\sigma^{Y}$ is
constant and does not evolve. We can therefore use $\sigma_{eq}$ as the input
to the neural network,
${\mathcal{NN}_{\phi^{*}}}={\mathcal{NN}_{\phi^{*}}(\sigma_{eq})}$. Figure 3
shows a schematic of the particular neural network architecture used for
modeling perfect viscoplasticity in this work. While different configurations
are possible in particular with regards to the selection of the number of
layers and neurons, due to the already low amount of training data (less than
100 data points depending on the resolution of stress-strain data points)
compared to the number of trainable parameters we show the results for a
network with 2 layers and 20 neurons in the following. We remark that slight
changes to this network architecture, i.e. 2 layers with 30 neurons, for
example, have not shown to yield different results.

Figure 3: A schematic of the neural network used for perfect elasto-
viscoplasticity.

####  2.2.2 NN-EVP with isotropic hardening

Once hardening effects are introduced into the constitutive model, the yield
function ${R(r)}$ (which will be not imposed, but discovered) evolves during
the deformation as a function of accumulated plastic strain ${r}$. As a
result, information pertaining to the evolution of yield function is also
required as an additional input to the dual potential neural network
$\mathcal{NN}_{\phi^{*}}$. Thus, we could define the dual potential neural
network as $\mathcal{NN}_{\phi^{*}}$
=${\mathcal{NN}_{\phi^{*}}(\sigma_{eq},{R(r)})}$. To this end, an additional
neural network ${\mathcal{NN}_{R}=\mathcal{NN}_{R}(r)}$ is required to track
the evolution of strain hardening. Notice that the output of hardening neural
network ${\mathcal{NN}_{R}(r)}$ could be used as an input to the dual
potential neural network ${\mathcal{NN}_{\phi^{*}}(\sigma_{eq},{R(r)})}$.
Here, instead of taking $\sigma_{eq}$ and ${R(r)}$ as two independent inputs
to ${\mathcal{NN}_{\phi^{*}}}$, we choose to take the ratio
$\dfrac{\sigma_{eq}}{{R(r)}}$ as a single input. This mechanistic assumption
informs the neural network with the physical constraint that $\sigma_{eq}$
needs to be scaled proportionally with the evolution of $R(r)$, imposing the
strain hardening effects as the ratio $\dfrac{\sigma_{eq}}{{R(r)}}$ controls
the elasto-viscoplastic transition as well as the rate of the evolution of
viscoplastic strain rate. Figure 4 shows a schematic of NN-EVP with isotropic
hardening. To remain consistent with the thermodynamics laws, the following
conditions are enforced on the hardening neural network $\mathcal{NN}_{R}$:

  * •

$\mathcal{NN}_{R}$ is monotonically increasing,

  * •

$\mathcal{NN}_{R}$ is always positive: ${\mathcal{NN}_{R}}>0$,

  * •

$\mathcal{NN}_{R}$ does not include the origin: $\mathcal{NN}_{R}(r=0)\neq 0$.

The neural network $\mathcal{NN}_{R}$ is positive and monotonically increasing
when the following conditions are met:

  * •

$x_{0}\geq 0,W_{l}\geq 0,b_{l}\geq 0$,

  * •

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}}_{l}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0},l=1,\ldots,L$,

  * •

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}}_{l}^{\prime}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0},l=1,\ldots,L$.

Notice that in contrast to the dual potential network
$\mathcal{NN}_{\phi^{*}}$, the hardening network $\mathcal{NN}_{R}$ does not
contain the origin and does not require convexity. This is due to the fact
that the yield function ${R(r)}$ has a nonzero value equal to the initial
yield at zero accumulated viscoplastic strain $r$. We also remark that since
$\mathcal{NN}_{R}$ is monotonically increasing, the reciprocal form of
$\dfrac{1}{\mathcal{NN}_{R}}$, used as input to $\mathcal{NN}_{\phi^{*}}$, is
monotonically decreasing [40]. In order to satisfy the conditions above, we
choose combinations of forms of ReLU, adaptive logistic or adaptive tanh
activation functions as follows:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}_{l}}^{\mathcal{NN}_{R}}(x)=\alpha_{1}\max(x,0)+\alpha_{2}\frac{1}{1+e^{-\beta
x}},$ (36)
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}_{l}}^{\mathcal{NN}_{R}}(x)=\alpha_{1}\max(x,0)+\alpha_{2}\frac{e^{\beta
x}-e^{-\beta x}}{e^{\beta x}+e^{-\beta x}},$ (37)
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}_{l}}^{\mathcal{NN}_{R}}(x)=\alpha_{1}\frac{1}{1+e^{-\beta
x}}+\alpha_{2}\frac{e^{\beta x}-e^{-\beta x}}{e^{\beta x}+e^{-\beta x}},$ (38)

where $\alpha_{1}$ and $\alpha_{2}$ denote the weights assigned to the
activation functions and applied to the hidden layers of each neural network.
The reason behind these mixed activation functions is to facilitate learning
the hardening response for larger strain amplitudes. Since the logistic
activation function saturates early at lower strain levels, the addition of
either a ReLu or tanh response compensates for stress-strain curvature at
later stages of hardening. The effect of variation in $\alpha_{1}$ and
$\alpha_{2}$ as well as the particular selection of the activation functions
on the flow response is discussed later in detail.

Figure 4: A schematic of the neural network used for elasto-viscoplasticity
with isotropic hardening.

####  2.2.3 Hall-Petch effects and grain size-aware NN-EVP

Hall-Petch effects and variation in grain size are concomitant with strong
anisotropy in viscoplastic flow responses. In particular, the initial yield
function of the material alters drastically depending on the average grain
size of the underlying microstructure. Particular functional forms for Hall-
Petch relationship vary in the literature. However, the most common form [46,
47, 48, 4, 31] can be written as follows:

$R_{0,HP}=\dfrac{H\mu\sqrt{b}}{\sqrt{d_{grain}}}$ (39)

where $b$ is the Burgers vector, $\mu$ is the shear modulus, $d_{grain}$ is
the average grain size and $H$ stands for the Hall-Petch coefficient usually
obtained from calibrations with experimental observations. Notice that Eq.
(39) provides a nonlinear relationship between the grain size and initial
yield function in non-logarithmic space, however, it can be shown that the
Hall-Petch stress does not always scale as the inverse square root of grain
size [2] and thus Eq. 39 is simply a specific function of many possible forms.
In order to incorporate the Hall-Petch effect into the NN-EVP framework, we
introduce a Hall-Petch neural network
$\mathcal{NN}_{HP}=\mathcal{NN}_{HP}({d_{grain}})$ in addition to the dual
potential $\mathcal{NN}_{\phi^{*}}$ and the hardening neural network
$\mathcal{NN}_{R}$. This allows us to learn and discover the Hall-Petch
relationship. To remain consistent with thermodynamics laws, the following
conditions are enforced on the Hall-Petch neural network $\mathcal{NN}_{HP}$:

  * •

$\mathcal{NN}_{HP}$ is monotonically decreasing,

  * •

$\mathcal{NN}_{HP}$ is always positive: ${\mathcal{NN}_{HP}}$ > 0,

  * •

$\mathcal{NN}_{HP}({d_{grain}\rightarrow 0})\rightarrow\infty$.

The neural network $\mathcal{NN}_{HP}$ is positive when the following
conditions are met:

  * •

$x_{0}\geq 0,W_{l}\geq 0,b_{l}\geq 0$,

  * •

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}}_{l}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0},l=1,\ldots,L$.

Here, we first choose a standard tanh activation function as follows and then
take the reciprocal of the network output to satisfy the conditions above:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{F}_{l}}^{\mathcal{NN}_{{HP}}(x)}=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}},$
(40) $\mathcal{NN}_{HP}\leftarrow\frac{1}{\mathcal{NN}_{HP}}.$ (41)

Figure 5: A schematic of neural network for elasto-viscoplasticity with
isotropic hardening and Hall-Petch effects.

Notice that since the Hall-Petch term only affects the initial yield function
and does not evolve during the deformation, an adaptive activation function
here is not necessary. Figure 5 shows a schematic of grain size-aware NN-EVP
architecture including the Hall-Petch effects. As illustrated in the figure,
the response of the Hall-Petch network $\mathcal{NN}_{HP}$ and the hardening
network $\mathcal{NN}_{R}$ are combined to
$\dfrac{\sigma_{eq}}{{\mathcal{NN}_{R(r)}+\mathcal{NN}_{HP}(d_{grain})}}$ and
then used as an input to the dual potential network $\mathcal{NN}_{\phi^{*}}$.

##  3 Results and Discussion

###  3.1 Implementation highlights

The NN-EVP framework is implemented in PyTorch [49]. PyTorch takes advantage
of automatic differentiation [50] facilitating automatic computation of local
gradients of the output of neural network with respect to its inputs.
Automatic differentiation also allows us to determine the Jacobian required
for solving the update step of implicit time-stepping using the Newton Raphson
(NR) method [51]. Details of the implementation are provided in the form of a
pseudo-code in algorithm 1. Due to its adaptive learning decay algorithm and
robustness, the AdamW optimizer [52] is utilized to optimize the parameters of
the three neural network used in this study. In addition, a cosine annealing
scheduler [53, 54] is utilized to further enhance the learning rate adaptivity
and therefore improve the robustness of the NN-EVP framework. The cosine
annealing scheduler is set with a starting learning rate of $1e^{-2}$ and
saturates to a learning rate of $1e^{-3}$.

The established framework is then trained to model perfect viscoplasticity as
well as an isotropic hardening response under uniaxial tensile loading. We
consider two different scenarios. First, the flow response is synthetically
generated using the power law and phenomenological Johnson-Cook hardening [55]
where in addition to fitting the existing data, the flow response is predicted
to larger strain amplitudes via extrapolations enabled by recovering the
trained neural networks. Next, experimentally measured data pertaining to
large deformations as a function of grain size are used to train the NN-EVP
framework presented herein. In all cases, a constant 11 component of the
strain rate $\dot{\epsilon}^{app}=1e^{-3}s^{-1}$ in X direction is applied. In
order to satisfy the uniaxial loading conditions, the 22 and 33 components of
the stress tensor are enforced to be zero (see algorithm 1). The NR solver
tolerance is set to $\mathcal{TOL_{NR}}=1e^{-6}$ for all simulations which is
reached within 2-4 iterations depending on the deformation history. Finally,
the mean squared error function is used to calculate the loss value in each
training epoch.

Algorithm 1 Grain size aware NN-EVP

Stress-strain response obtained from power law or experimental data:
$\begin{bmatrix}\bm{\epsilon}^{True}&,&\bm{\sigma}^{True}\\\ \end{bmatrix}$

Establish:

  * •

Initialize positive, monotonically increasing, input convex neural network
$\mathcal{NN_{\phi^{*}}}$ and its corresponding optimizer
$\mathcal{O_{NN_{\phi^{*}}}}$ and parameters
$\mathcal{\theta_{NN_{\phi^{*}}}}$,

  * •

Initialize positive, monotonically decreasing neural network
$\mathcal{NN_{R}}$ and its corresponding optimizer $\mathcal{O_{NN_{R}}}$ and
parameters $\mathcal{\theta_{NN_{R}}}$,

  * •

Initialize positive, monotonically decreasing neural network
$\mathcal{NN_{HP}}$ and its corresponding optimizer $\mathcal{O_{NN_{HP}}}$
and parameters $\mathcal{\theta_{NN_{HP}}}$ ,

  * •

Define applied strain rate tensor corresponding to the uniaxial deformation :
$\begin{bmatrix}\dot{\epsilon}^{app}&0&0\\\
0&-\dfrac{\dot{\epsilon}^{app}}{2}&0\\\
0&0&-\dfrac{\dot{\epsilon}^{app}}{2}\end{bmatrix}$,

  * •

Define total number of time steps $N_{Steps}$ , training epochs $N_{Epochs}$
and NR tolerance $\mathcal{TOL_{NR}}$,

  * •

Define the material properties, Young’s modulus $E$ and Poisson’s ratio $\nu$
as well as the average grain size $d_{grain}$.

Begin:

Calculate the elastic stiffness $\mathbb{C}$:
$\mathbb{C}_{ijkl}=\frac{E\nu}{(1+\nu)(1-2\nu)}\delta_{ij}\delta_{kl}+\frac{E}{2(1+\nu)}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$.
Initialize strain and stress tensors: $\bm{\epsilon}^{t=0}=\bm{0}$,
$\bm{\epsilon}^{vp,t=0}=\bm{0}$ , $\bm{\sigma}^{t=0}=\bm{0}$

for $epoch=0,N_{epochs}$ do Initialize loss:
$\mathcal{L}(\mathcal{O_{NN_{\phi^{*}}}},\mathcal{O_{NN_{R}}},\mathcal{O_{NN_{HP}}})=0$

for $t=0,N_{steps}$ do Update the total strain at current time step
$\bm{\epsilon}^{t+1}$ using applied strain rate:
$\bm{\epsilon}^{t+1}\leftarrow\bm{\epsilon}^{t}+\bm{\dot{\epsilon}}\Delta t$
Calculate trial stress using elastic strain:
$\bm{\sigma}^{t+1}=\mathbb{C}:\bm{\epsilon}^{t+1}$ Calculate viscoplastic
strain rate using dual potential neural network $\mathcal{NN_{\phi^{*}}}$:
${\dot{\bm{\epsilon}}^{vp,t+\Delta
t}}=\dfrac{\partial{\mathcal{NN}_{\phi^{*}}(\bm{\sigma}^{t+1},R^{t+1},d_{grain})}}{\partial{\bm{\sigma}}}$
$\triangleright$ Using Torch.autograd.grad() Define the residuals:
$\triangleright$ Use Voigt notation
$\sigma_{ij}\rightarrow\sigma_{k},\mathbb{C}_{ijkl}\rightarrow\mathbb{C}_{kl},k,l=1,6$
[]
$\mathcal{X}_{k}=\begin{bmatrix}{\mathcal{R}es_{k}(\sigma^{t+1})}=\sigma_{k}^{t+1}-\mathbb{C}_{kl}({\epsilon_{l}}^{t+\Delta
t}-{\epsilon_{l}}^{vp,t}-{\dot{\epsilon}_{l}}^{vp,t+\Delta t}\Delta t)=0\\\
{\mathcal{R}es_{k}(\sigma^{t+1})}=\sigma^{t+1}_{k}=0\,\qquad\text{for }k\neq
1\\\ \end{bmatrix}$

Deploy NR with Jacobian $\mathcal{J}$ to solve for stress: $\triangleright$ NR
iterations denoted by i

while $\norm{\mathcal{X}}<\mathcal{TOL_{NR}}$ do []
$\displaystyle\sigma_{k}^{t+1,i+1}=\sigma_{k}^{t+1,i}-\mathcal{J}_{kl}^{-1}{\mathcal{X}_{l}(\sigma^{t+1,i})},\mathcal{J}_{kl}=\dfrac{\partial{\mathcal{X}_{k}}}{\partial\sigma_{l}}|_{\sigma^{t+1,i}}$

end while$\triangleright$ End of loop for NR solver Update optimization loss:
$\mathcal{L}\leftarrow\mathcal{L}+\norm{\bm{\sigma}^{t+1}-\bm{\sigma}^{True,t+1}}^{2}$

end for$\triangleright$ End of loop for time increments Calculate gradients:
$\bm{G}_{\mathcal{NN}_{\phi}^{*}}=\nabla_{{\theta_{\phi^{*}}}}\mathcal{L},\bm{G}_{\mathcal{R}}=\nabla_{{\theta_{R}}}\mathcal{L},\bm{G}_{\mathcal{HP}}=\nabla_{{\theta_{HP}}}\mathcal{L}$
Update NN parameters: $\mathcal{O}_{R},\mathcal{O}_{\phi},\mathcal{O}_{M}$
$\triangleright$ Use gradient clipping enforce positive outputs

end for$\triangleright$ End of loop for training epochs

End

###  3.2 Rediscovering the power law for perfect viscoplasticity

In order to predict perfect viscoplasticity in Cu, the power law equation with
a constant yield function $\sigma^{Y}$ and rate sensitivities of $n=$10, 20,
and 100 are used to generate the synthetic stress-strain data. Since strain
hardening is not involved, only the dual potential neural network
${\mathcal{NN}_{\phi^{*}}(\sigma_{eq}=\sigma)}$ is utilized. The left plot of
Fig. 6(a) shows the evolution of normalized loss over 200 training epochs for
a rate sensitivity of $n=$10\. The training error saturates to a value of
$1e^{-5}$ after around 100 iterations. Note that the fluctuations observed in
the loss evolution are innate to the optimizer algorithm regardless of the
learning rate value. The right plot of Fig. 6(a) illustrates the trained model
output using the ground truth provided by the power law. The transparent lines
represent the training history of the model during the optimization process.
Figures 6(b) and 6(c) highlight the loss evolution and training response for
power laws with rate sensitives of $n=$20 and $n=$100 respectively. As the
rate sensitivity increases, a sharper transition from elastic to viscoplastic
deformation is observed, making it more challenging for the neural network to
fit the data at the transition regime. An increased number of training epochs
from 200 to 500 is an indication of such behavior. Note that while the
training loss is still reducing slowly after around 500 epochs, to maintain
consistency and computational efficiency, we hereafter set the maximum number
of epochs $N_{Epochs}$=500 for all the simulations.

Figure 6: Model fit for perfect elasto-viscoplasticity with power-law rate
sensitivities of (a) n=10, (b) n=20 and (c) n=100. Transparent lines indicate
the incremental training history.

###  3.3 Training via phenomenological hardening

Synthetic stress-strain data for isotropic hardening can be generated using
either physics-based or phenomenological hardening models. Physics-based
hardening laws are more suited for mesoscale or micromechanical simulations
because they consider the evolution of dislocation substructures and
dislocation densities [4] and often require an evolved implementation
framework which is computationally intensive. On the other hand,
phenomenological hardening models such as Voce [56], Peirce, Asaro, and
Needleman (PAN) [57] or Johnson-Cook [55] are simpler in implementation and
thus computationally more efficient for modeling macroscopic behavior where
the details pertaining to the underlying microstructure is not considered. In
the present work, we use the Johnson-Cook isotropic hardening model with the
functional form written as follows [55]:

$R(r)=[A+Br^{n}][1+C\log\dfrac{r}{r^{*}}][1-{T^{*}}^{m}],$ (42)
$T^{*}=\dfrac{T-T_{0}}{T_{m}-T_{0}},$ (43)

where $A$, $B$, $C$, $n$ and $m$ are the hardening parameters associated with
the Johnson-Cook model. Note that since we are considering deformations at
room temperature $T_{0}$, the term $T^{*}$ that includes the effects of
melting temperature $T_{m}$ vanishes. The values corresponding to the Johnson-
Cook model parameters for Cu are listed in Table 1 [55]. Here both the dual
potential network ${\mathcal{NN}_{\phi^{*}}(\sigma,R(r))}$ and hardening
neural network ${\mathcal{NN}_{R(r)}}$ are utilized for training. The dual
potential network uses the SoftPlus activation function while the hardening
neural network ${\mathcal{NN}_{R(r)}}$ uses the adaptive logistic activation
function obtained by setting $\alpha_{1}=$0 and $\alpha_{2}=$1 in Eq. 36.
Figure 7 shows the loss evolution and trained model. Young’s modulus and
Poisson’s ratio of 130 GPa and 0.34 are used for the elastic response. The
loss evolution experiences an abrupt drop to a normalized value of around
$5e^{-4}$ and saturates after around 100 epochs, indicating a fast convergence
of the proposed framework.

$A[MPa]$ | $B[MPa]$ | $C$ | $m$ | $n$
---|---|---|---|---
90 | 292 | 0.31 | 0.025 | 1.09
Table 1: Johnson-Cook hardening model parameters for Cu. Figure 7: Model fit
for elasto-viscoplasticity with isotropic Johnson-Cook hardening. Transparent
lines indicate the incremental training history.

####  3.3.1 Extrapolation of flow response to larger strain amplitudes

One of the powerful advantages (and challenges) of algorithms that are
developed via neural network is the ability to extrapolate beyond the data
available to the networks during the training process. Once a model is trained
for a specific set of data, the recovered trained model is able to predict the
behavior for any given input within the trained framework and beyond the data
observed by the model. Here, we restrict the model during training for
predictions up to 0.5% total strain and recover the model and its trained
parameters to examine the ability of the model to extrapolate the flow
response up to 2% total strain. This benchmark provides us with a train-test
configuration consisting of 25% training and 75% testing data. Figure 8 shows
the extrapolations of stress-strain response for different ratios of the mixed
activation functions ReLU, adaptive logistic, and adaptive tanh with variable
combination weights $\alpha_{1}$ and $\alpha_{2}$ in Eqs. (36), (37) and (38).

Figure 8: Extrapolations of flow response with an isotropic hardening response
via phenomenological Johnson-Cook model as a function of selected activation
functions and their combination weights $\alpha_{1}$ and $\alpha_{2}$ with (a)
20% logistic - 80% ReLU, (b) 80% logistic - 20% tanh, (c) 50% logistic - 50%
tanh, (d) 100% tanh and (e) 80% tanh - 20% ReLU. The dashed red line shows the
ground truth. The solid grey line indicates the trained NN model and the solid
blue line implies the extrapolations of the flow response using the NN model.

The logistic activation function saturates rapidly at low strain amplitudes,
hindering the model from properly capturing the hardening response at higher
strain amplitudes. To address this issue, tanh, and ReLU activation functions
are combined with logistic activation functions using different weights to
predict the flow response beyond the training domain. Note that the addition
of the ReLU activation function should be done with caution to avoid excessive
hardening behavior shown in Fig. 8(a). However, it is worthwhile to mention
that sudden increase in strain hardening is well observed in case of
deformation twinning in compression tests [58] and thus, a mixed activation
function with large ReLU weights seems to be a promising approach for
capturing the rapid hardening behavior upon twinning formations. Among various
configurations shown in Fig. 8(a)-(e), a mixed activation function with 80%
adaptive tanh and 20% ReLU best predicts the hardening behavior up to 2%
strain.

###  3.4 Training large plastic deformations via experimental data

So far the capability of the NN-EVP framework in training and extrapolating
the flow response at small deformations is well demonstrated. In this section,
we elaborate on training the flow response on experimentally measured data at
large plastic deformations. One of the challenges involved in training with
experimental data is the limitation of stress-strain data points as well as
their frequency and form of occurrence through the deformation history. Access
to the experimental data reported in the literature usually involves the
extraction of data from stress-strain images using digitizing software.
Typically data extraction is based on manual user input via the coordinates
representing the image. Creating a pair of lists for the model output and
measured data matching at the corresponding strain levels and necessary for a
one-by-one comparison is challenging. Generating data with equal spacing is
also not an option since adaptive time stepping is required for a
computationally efficient framework with fewer time increments that represent
the same curve with large deformations up to 10% strain.

Thus, to address this issue and to enable access to the experimental data
points at arbitrary strain amplitudes with arbitrary and adaptive time
stepping, a nonlinear interpolation preprocessing step via the Scipy [59]
library is used prior to training. Figure 9 shows an example of a continuous
Scipy interpolation for experimentally obtained data in Ni with uneven spacing
and arbitrary strain amplitudes. Once the Scipy interpolation is applied, the
stress values corresponding to the desired strain levels are readily evaluated
using the obtained interpolator. Figures 10(a) and 10(b) present the trained
models using measured data for small and large deformations up to 1% and 10%
strain in Ni. The model shows promising performance in capturing the elasto-
viscoplastic transition in small strains as well as prediction of hardening
curvature under large deformations. Note that since most of the experimental
data existing in the literature ignore the elastic response, an additional
step based on the current value of viscoplastic strain is imposed in algorithm
1 to ignore the elastic response when computing the loss between the training
output and measured response.

Figure 9: Interpolation of measured data in Ni using SciPy.  Figure 10: Model
fit for measured flow response in Ni for (a) small deformation to a
macroscopic strain of 1% and (b) large deformation to a macroscopic strain of
10%.

###  3.5 Grain size-aware flow response and Hall-Petch discovery

After validating the ability of the NN-EVP model in training experimental data
to large strains with arbitrary data frequency, we take a step further to
train the flow response at large deformations as a function of grain size and
aim to discover the Hall-Petch relationship without considering a particular
functional form. To this end, the measured flow response of Cu as a function
of grain size and for average grain sizes of 2.1 $\mu m$, 3.4 $\mu m$, 7.1
$\mu m$ and 15 $\mu m$, as shown in Fig. 11, is utilized for training the
model including the grain size effects. Here, in addition to the dual
potential network ${\mathcal{NN}_{\phi^{*}}(\sigma,R(r))}$ and hardening
network ${\mathcal{NN}_{R(r)}}$, the Hall-Petch neural network
${\mathcal{NN}_{HP(d_{grain})}}$ is also used in the model. In order to
increase the training speed, adaptive time stepping was implemented with
$\Delta t|_{t+\Delta t}=1.15\Delta t|_{t}$ , resulting in 50 increments per
stress-strain curve. Since the time stepping is adaptive, the training losses
corresponding to each data point are scaled proportionally with the time
increment to that data point to balance the loss evaluation. This is required
due to the uneven data spacing resulting from irregular time steps where the
density of training data is much smaller within the large deformation regime
above 2% strain.

Figure 11: Variation in flow response of Cu under large deformation as a
function of grain size.

First, each curve is trained individually and independent of one another to
test the ability of the model in training the flow response with varying
stress levels. Results corresponding to single-curve fits are shown in Fig.
12. Since one curve is trained at a time, it is easier for the model to learn
the flow response and thus the convergence is relatively fast. Next, a
quaternary configuration (with all four stress-strain responses) including 200
data points is utilized for training. Training on binary and ternary
configurations with two and three curves consisting of 100 and 150 data points
are also provided in the supplementary data shown in Figs. 15 and 16. Note
that since multiple curves with contrasting stress levels and hardening
curvatures are trained in parallel, it is more challenging for the model to
converge to an optimal solution as evident in the behavior of loss evolution.
Also, due to more data points the training process is computationally more
demanding.

Ultimately, we aim to discover the Hall-Petch relationship which describes the
dependence of the initial yield on the grain size. To this end, after training
the model for grain sizes of 2.1 $\mu m$, 3.4 $\mu m$, 7.1 $\mu m$ and 15 $\mu
m$, the Hall-Petch neural network ${\mathcal{NN}_{HP(d_{grain})}}$ is
recovered to extrapolate the initial flow response for grain sizes below 2.1
$\mu m$ and beyond 7.1 $\mu m$ with average grain sizes and their
corresponding Hall-Petch stress values listed in Table 2. The discovered Hall-
Petch behavior is shown in Fig. 14. These results indicate the ability of the
model to capture the expected linear Hall-Petch relationship in log-log space,
discovering it without any underlying assumptions on its functional form. Note
that regardless of the particular material used herein, the proposed NN-EVP
model should be able to predict the Hall-Petch relationship for a wide range
of materials with arbitrary grain size and stress-strain response.
Extrapolation of the Hall-Petch relationship is a complicated topic that we
aim to study in the future.

Figure 12: Single-fit training of flow response in Cu for grain sizes of (a) 2.1 $\mu m$, (b) 3.4 $\mu m$, (c) 7.1 $\mu m$ and (d) 15 $\mu m$. Figure 13: Multi-fit training of flow response in Cu for quaternary configuration of 2.1-3.4-7.1-15 $\mu m$. Grain size $[\mu m]$ | Hall-Petch stress $[MPa]$
---|---
0.5 | 182.3671
1.0 | 91.1886
2.1 | 43.4293
3.4 | 26.8331
5.0 | 18.2582
7.1 | 12.8734
10.0 | 9.1615
15.0 | 6.1435
20.0 | 4.645
50.0 | 2.0286
100 | 1.2712
250 | 0.9902
500 | 0.8907
Table 2: Hall-Petch contributions in Cu via extrapolations of initial flow
response as function of grain size. Figure 14: Discovery of Hall-Petch
relationship in Copper using NN-EVP in non-logarithmic (on the left) and log-
log (on the right) scales.

##  4 Summary and conclusions

The flow response in metals and metallic alloys differs drastically depending
on the average grain size that forms the underlying microstructure. Modeling
the viscoplastic flow response of metals is usually concomitant with the
assumption of particular constitutive models in the form of a power law and
phenomenological hardening formulations. These hardening laws often include
parameters that involve manual calibrations with experimental data obtained
from uniaxial tests. Machine learning models that can automate such
calibrations generally rely on large amounts of data and often require a large
number of iterations in order for them to be sufficiently accurate and
predictive. For instance, the fitting process via genetic algorithms involve a
large number of functional calls within four orders of magnitude to the black-
box in order to obtain the flow response [60, 61]. To remove the need for
these big datasets and large number of iterations in computationally expensive
training, we propose a data-driven elasto-viscoplasticity framework based on
neural networks called NN-EVP that leverages PyTorch’s high-performance ML
library. The proposed approach is tested and trained using both synthetic and
experimentally measured uniaxial data in the context of limited data
availability.

The developed NN-EVP framework was adopted to train elasto-viscoplastic flow
response of metals at both small and large deformations and as a function of
grain size. First, rate sensitivity-dependent perfect viscoplasticity was
discovered using the power law with no hardening effects. Next, synthetically
generated deformation responses via a Johnson-Cook phenomenological hardening
law were used to train and test the approach on its ability to extrapolate the
flow response at strain amplitudes beyond the observed training data. Next,
the framework was applied to train large deformations obtained from
experimental tests with limited stress-strain data availability. Finally,
simultaneous training of multiple grain size-dependent stress-strain curves
enabled us to obtain a grain size-aware flow behavior and discover the Hall-
Petch relationship pertaining to the grain size strengthening effects.

The proposed NN-EVP model presented herein takes a further step in the
prediction of flow responses in metals and improves the computational
efficiency of structure-property-relationship simulations of metallic
materials. The findings of the present work provide insights into the
versatility and flexibility of data-driven constitutive modeling which can be
adapted for a wide range of materials exhibiting complex large deformation
behavior. This motivates future work in incorporating the current model into
finite elements for efficient full-field modeling of large deformations under
arbitrary loading and boundary conditions.

##  5 Data availability

Data and Python script supporting the findings of this study are available
from the corresponding author upon request.

##  6 Acknowledgments

JF, AE and NB gratefully acknowledge support by the Air Force Office of
Scientific Research under award number FA9550-22-1-0075.

##  Appendix A Supplementary data

Figure 15: Multi-fit training of flow response in Cu for binary configurations
of (a) 2.1-3.4 $\mu m$, (b) 3.4-15 $\mu m$, (c) 2.1-7.1 $\mu m$ and (d) 2.1-15
$\mu m$. Figure 16: Multi-fit training of flow response in Cu for Ternary
configurations of (a) 2.1-3.4-7.1 $\mu m$ and (b) 3.4-7.1-15 $\mu m$.

## References

  * [1] N. Hansen, Hall–petch relation and boundary strengthening, Scripta Materialia 51 (8) (2004) 801–806, viewpoint set no. 35. Metals and alloys with a structural scale from the micrometer to the atomic dimensions.  doi:https://doi.org/10.1016/j.scriptamat.2004.06.002.
URL https://www.sciencedirect.com/science/article/pii/S1359646204003434

  * [2] D. Dunstan, A. Bushby, Grain size dependence of the strength of metals: The hall–petch effect does not scale as the inverse square root of grain size, International Journal of Plasticity 53 (2014) 56–65.
  * [3] R. B. Figueiredo, M. Kawasaki, T. G. Langdon, Seventy years of hall-petch, ninety years of superplasticity and a generalized approach to the effect of grain size on flow stress, Progress in Materials Science (2023) 101131.
  * [4] A. Eghtesad, J. D. Shimanek, S.-L. Shang, R. Lebensohn, M. Knezevic, Z.-K. Liu, A. M. Beese, Density functional theory-informed dislocation density hardening within crystal plasticity: Application to modeling deformation of ni polycrystals, Computational Materials Science 215 (2022) 111803.
  * [5] R. B. Figueiredo, T. G. Langdon, Deformation mechanisms in ultrafine-grained metals with an emphasis on the hall–petch relationship and strain rate sensitivity, Journal of Materials Research and Technology 14 (2021) 137–159.
  * [6] D. Agius, A. Kareer, A. Al Mamun, C. Truman, D. M. Collins, M. Mostafavi, D. Knowles, A crystal plasticity model that accounts for grain size effects and slip system interactions on the deformation of austenitic stainless steels, International Journal of Plasticity 152 (2022) 103249.
  * [7] A. Lakshmanan, M. Yaghoobi, K. S. Stopka, V. Sundararaghavan, Crystal plasticity finite element modeling of grain size and morphology effects on yield strength and extreme value fatigue response, journal of materials research and technology 19 (2022) 3337–3354.
  * [8] L. Zhang, Y. Shibuta, Inverse hall-petch relationship of high-entropy alloy by atomistic simulation, Materials Letters 274 (2020) 128024.
  * [9] M. Zhang, Z. Rao, T. Xu, L. Fang, Quantifying the influence of grain boundary activities on hall-petch relation in nanocrystalline cu by using phase field and atomistic simulations, International Journal of Plasticity 135 (2020) 102846\.
  * [10] B. Le, J. Yvonnet, Q.-C. He, Computational homogenization of nonlinear elastic materials using neural networks, International Journal for Numerical Methods in Engineering 104 (12) (2015) 1061–1084.
  * [11] J. N. Fuhg, C. Böhm, N. Bouklas, A. Fau, P. Wriggers, M. Marino, Model-data-driven constitutive responses: application to a multiscale computational framework, International Journal of Engineering Science 167 (2021) 103522.
  * [12] K. A. Kalina, L. Linden, J. Brummund, M. Kästner, Fe ann: an efficient data-driven multiscale approach based on physics-constrained neural networks and automated data mining, Computational Mechanics 71 (5) (2023) 827–851.
  * [13] N. N. Vlassis, R. Ma, W. Sun, Geometric deep learning for computational mechanics part i: anisotropic hyperelasticity, Computer Methods in Applied Mechanics and Engineering 371 (2020) 113299.
  * [14] J. N. Fuhg, N. Bouklas, On physics-informed data-driven isotropic and anisotropic constitutive models through probabilistic machine learning and space-filling sampling, Computer Methods in Applied Mechanics and Engineering 394 (2022) 114915.
  * [15] D. K. Klein, M. Fernández, R. J. Martin, P. Neff, O. Weeger, Polyconvex anisotropic hyperelasticity with neural networks, Journal of the Mechanics and Physics of Solids 159 (2022) 104703.
  * [16] K. A. Kalina, L. Linden, J. Brummund, P. Metsch, M. Kästner, Automated constitutive modeling of isotropic hyperelasticity based on artificial neural networks, Computational Mechanics (2022) 1–20.
  * [17] J. N. Fuhg, N. Bouklas, R. E. Jones, Learning hyperelastic anisotropy from data via a tensor basis neural network, Journal of the Mechanics and Physics of Solids 168 (2022) 105022.
  * [18] V. Taç, K. Linka, F. Sahli-Costabal, E. Kuhl, A. B. Tepole, Benchmarking physics-informed frameworks for data-driven hyperelasticity, Computational Mechanics (2023) 1–17.
  * [19] V. Taç, M. K. Rausch, F. S. Costabal, A. B. Tepole, Data-driven anisotropic finite viscoelasticity using neural ordinary differential equations, Computer Methods in Applied Mechanics and Engineering 411 (2023) 116046\.
  * [20] F. As’ ad, C. Farhat, A mechanics-informed neural network framework for data-driven nonlinear viscoelasticity, in: AIAA SCITECH 2023 Forum, 2023, p. 0949\.
  * [21] M. Flaschel, S. Kumar, L. De Lorenzis, Automated discovery of generalized standard material models with euclid, Computer Methods in Applied Mechanics and Engineering 405 (2023) 115867.
  * [22] K. Upadhyay, J. N. Fuhg, N. Bouklas, K. Ramesh, Physics-informed data-driven discovery of constitutive models with application to strain-rate-sensitive soft materials, arXiv preprint arXiv:2304.13897 (2023).
  * [23] T. Kadeethum, D. O’Malley, J. N. Fuhg, Y. Choi, J. Lee, H. S. Viswanathan, N. Bouklas, A framework for data-driven solution and parameter estimation of pdes using conditional generative adversarial networks, Nature Computational Science 1 (12) (2021) 819–829.
  * [24] D. K. Klein, R. Ortigosa, J. Martínez-Frutos, O. Weeger, Finite electro-elasticity with physics-augmented neural networks, Computer Methods in Applied Mechanics and Engineering 400 (2022) 115501.
  * [25] D. Yao, Y.-c. Duan, M.-y. Li, Y.-p. Guan, Hybrid identification method of coupled viscoplastic-damage constitutive parameters based on bp neural network and genetic algorithm, Engineering Fracture Mechanics 257 (2021) 108027\.
  * [26] C. Bonatti, B. Berisha, D. Mohr, From cp-fft to cp-rnn: Recurrent neural network surrogate model of crystal plasticity, International Journal of Plasticity 158 (2022) 103430.
  * [27] O. Ibragimova, A. Brahme, W. Muhammad, J. Levesque, K. Inal, A new ann based crystal plasticity model for fcc materials and its application to non-monotonic strain paths, International Journal of Plasticity 144 (2021) 103059\.
  * [28] O. Ibragimova, A. Brahme, W. Muhammad, D. Connolly, J. Lévesque, K. Inal, A convolutional neural network based crystal plasticity finite element framework to predict localised deformation in metals, International Journal of Plasticity 157 (2022) 103374.
  * [29] J. N. Heidenreich, M. B. Gorji, D. Mohr, Modeling structure-property relationships with convolutional neural networks: Yield surface prediction based on microstructure images, International Journal of Plasticity 163 (2023) 103506.
  * [30] A. Nascimento, S. Roongta, M. Diehl, I. J. Beyerlein, A machine learning model to predict yield surfaces from crystal plasticity simulations, International Journal of Plasticity 161 (2023) 103507.
  * [31] A. Eghtesad, Q. Luo, S.-L. Shang, R. A. Lebensohn, M. Knezevic, Z.-K. Liu, A. M. Beese, Machine learning-enabled identification of micromechanical stress and strain hotspots predicted via dislocation density-based crystal plasticity simulations, International Journal of Plasticity (2023) 103646.
  * [32] J. N. Fuhg, L. van Wees, M. Obstalecki, P. Shade, N. Bouklas, M. Kasemer, Machine-learning convex and texture-dependent macroscopic yield from crystal plasticity simulations, Materialia 23 (2022) 101446.
  * [33] M. Stoffel, F. Bamer, B. Markert, Neural network based constitutive modeling of nonlinear viscoplastic structural response, Mechanics Research Communications 95 (2019) 85–88.
  * [34] D. Huang, J. N. Fuhg, C. Weißenfels, P. Wriggers, A machine learning based plasticity model using proper orthogonal decomposition, Computer Methods in Applied Mechanics and Engineering 365 (2020) 113008.
  * [35] M. Stoffel, F. Bamer, B. Markert, Deep convolutional neural networks in structural dynamics under consideration of viscoplastic material behaviour, Mechanics Research Communications 108 (2020) 103565.
  * [36] S. B. Tandale, F. Bamer, B. Markert, M. Stoffel, Physics-based self-learning recurrent neural network enhanced time integration scheme for computing viscoplastic structural finite element response, Computer Methods in Applied Mechanics and Engineering 401 (2022) 115668.
  * [37] L. Benabou, Implementation of a neural network into a user-material subroutine for finite element simulation of material viscoplasticity, Journal of Engineering Materials and Technology 143 (4) (2021).
  * [38] N. Brodnik, C. Muir, N. Tulshibagwale, J. Rossin, M. Echlin, C. Hamel, S. Kramer, T. Pollock, J. Kiser, C. Smith, et al., Perspective: Machine learning in experimental solid mechanics, Journal of the Mechanics and Physics of Solids 173 (2023) 105231.
  * [39] N. N. Vlassis, W. Sun, Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening, Computer Methods in Applied Mechanics and Engineering 377 (2021) 113695.
  * [40] J. N. Fuhg, C. M. Hamel, K. Johnson, R. Jones, N. Bouklas, Modular machine learning-based elastoplasticity: Generalization in the context of limited data, Computer Methods in Applied Mechanics and Engineering 407 (2023) 115930\.
  * [41] W. Noll, B. D. Coleman, W. Noll, The thermodynamics of elastic materials with heat conduction and viscosity, The Foundations of Mechanics and Thermodynamics: Selected Papers (1974) 145–156.
  * [42] J. Lemaitre, J.-L. Chaboche, Mechanics of solid materials, Cambridge university press, 1994.
  * [43] M. Knezevic, M. Zecevic, I. J. Beyerlein, R. A. Lebensohn, A numerical procedure enabling accurate descriptions of strain rate-sensitive flow of polycrystals within crystal visco-plasticity theory, Computer Methods in Applied Mechanics and Engineering 308 (2016) 468–482.
  * [44] R. Gnanasambandam, B. Shen, J. Chung, X. Yue, et al., Self-scalable tanh (stan): Faster convergence and better generalization in physics-informed neural networks, arXiv preprint arXiv:2204.12589 (2022).
  * [45] A. D. Jagtap, G. E. Karniadakis, How important are activation functions in regression and classification? a survey, performance comparison, and future directions, arXiv preprint arXiv:2209.02681 (2022).
  * [46] N. Hansen, Hall–petch relation and boundary strengthening, Scripta materialia 51 (8) (2004) 801–806.
  * [47] G. Saada, Hall–petch revisited, Materials Science and Engineering: A 400 (2005) 146–149.
  * [48] Z. C. Cordero, B. E. Knight, C. A. Schuh, Six decades of the hall–petch effect–a survey of grain-size strengthening studies on pure metals, International Materials Reviews 61 (8) (2016) 495–512.
  * [49] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems 32 (2019).
  * [50] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, A. Lerer, Automatic differentiation in pytorch (2017).
  * [51] T. J. Ypma, Historical development of the newton–raphson method, SIAM review 37 (4) (1995) 531–551.
  * [52] I. Loshchilov, F. Hutter, Decoupled weight decay regularization, arXiv preprint arXiv:1711.05101 (2017).
  * [53] I. Loshchilov, F. Hutter, Sgdr: Stochastic gradient descent with warm restarts, arXiv preprint arXiv:1608.03983 (2016).
  * [54] Z. Lu, J. Wang, J. Song, Multi-resolution csi feedback with deep learning in massive mimo system, in: ICC 2020-2020 IEEE International Conference on Communications (ICC), IEEE, 2020, pp. 1–6.
  * [55] G. R. Johnson, A constitutive model and data for materials subjected to large strains, high strain rates, and high temperatures, Proc. 7th Inf. Sympo. Ballistics (1983) 541–547.
  * [56] E. Voce, The relationship between stress and strain for homogeneous deformation, Journal of the Institute of Metals 74 (1948) 537–562.
  * [57] D. Peirce, R. Asaro, A. Needleman, An analysis of nonuniform and localized deformation in ductile single crystals, Acta metallurgica 30 (6) (1982) 1087–1119.
  * [58] A. A. Salem, S. R. Kalidindi, R. D. Doherty, Strain hardening of titanium: role of deformation twinning, Acta Materialia 51 (14) (2003) 4225–4237.
  * [59] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17 (3) (2020) 261–272.
  * [60] T. Furukawa, G. Yagawa, Inelastic constitutive parameter identification using an evolutionary algorithm with continuous individuals, International Journal for Numerical Methods in Engineering 40 (6) (1997) 1071–1090.
  * [61] A. Jenab, I. S. Sarraf, D. E. Green, T. Rahmaan, M. J. Worswick, The use of genetic algorithm and neural network to predict rate-dependent tensile flow behaviour of aa5182-o sheets, Materials & Design 94 (2016) 262–273.