# Deterministic calculation of elasto-plastic stress-strain behavior from
arbitrary deformation fields

Benjamin C. Cameron
Department of Civil and Environmental Engineering
Massachusetts Institute of Technology
<EMAIL_ADDRESS>
C. Cem. Tasan
Department of Materials Science and Engineering
Massachusetts Institute of Technology
<EMAIL_ADDRESS>

###### Abstract

A system of partial differential equations is derived to compute the full-
field stress from an observed strain field when the plastic component of the
material constitutive equation is unknown. These equations generalize
previously proposed equations that are valid in the case where the elastic
strain is negligible. The system of equations assume that the constitutive
relations are isotropic, but otherwise make few assumptions and can be
directly applied (without modification) to cases of finite deformation, non-
linear elasticity and plasticity, compressible materials, rate dependent
materials, and a variety of different yield surface shapes and hardening laws.
Unlike the prior linear hyperbolic partial differential equations, this system
of equations is non-linear and time dependent. The newly proposed equations
can be used to solve for the stress field (and hence constitutive equation),
for arbitrary geometries and loading conditions when the full-field strain is
known. This generalization resolves a significant challenge in computing the
stress during elasto-plastic deformations: some regions of the material will
be elastically unloading while others are undergoing significant plastic
deformation. This problem significantly limited applications of the previously
proposed equations, but is fully resolved here.

A two-dimensional case study of necking in a uniaxial tensile test specimen is
investigated to illustrate the method. The governing equations are numerically
solved using an algorithm based on on the finite volume method developed in an
accompanying work. This is validated against the solution to the forward
problem and shown to give accurate results (within numerical error of the true
solution). This case study thus demonstrates that the developed approach has
unique capabilities, including the first theoretically exact solution to this
well studied problem.

_K_ eywords B elastic-plastic material $\cdot$ B constitutive behavior $\cdot$
C numerical algorithms $\cdot$ inverse problem

##  1 Introduction

The rapid increase in the availability and quality of full-field deformation
measurements has created a major imbalance in the available strain and stress
data. Full two- or three-dimensional deformation fields can be obtained by
digital image or volume correlation analyses [1, 2, 3], benefiting from
tomographic imaging techniques where needed [4, 5, 6]. However, force or
stress measurements are largely restricted to measurements using load cells,
in-situ stress transducers [7], photoelasticity approaches [8] or destructive
measurements [9]. This poses major limitations in the analysis of
heterogeneous deformations, as in the case of necking, shear band formation,
or high strain-rate deformation. The deformation itself provides only limited
insights regarding material behavior or properties, whereas knowing both the
strain and the corresponding stress fields would determine the constitutive
relations that fully define the mechanical response.

This motivates the study of the inverse problem: computing the stress from the
strain when the constitutive equations remain unknown. There has been some
success in solving this problem in cases with homogeneous constitutive
equations where the functional form is known in advance and is parameterized
by a only a few variables. As the number of parameters required to specify the
full strain field is much larger, the problem is generally over determined and
no solution is available that is consistent with force balance, the
constitutive equation form, and the observed strain field. To address this,
optimization approaches have been applied to minimize error in one or more of
these constraints. There have been a wide range of methods that take this
approach, including the virtual field method [10, 11, 12, 13], the equilibrium
gap method [14, 15], the constitutive equation gap method [11], finite element
model updating [16, 17, 18], and others (e.g. [19, 20, 21, 22]). These
approaches appear to be well suited to noisy or incomplete data sets, however,
depending on the application the errors can be significant (typically ranging
between 1 and 50 % depending on the application) and the computational cost
can be high. In particular these challenges are exacerbated when the material
properties are heterogeneous, as is often the case when considering material
microstructures, graded materials, composite materials, or materials with
heterogeneities arising from various manufacturing processes.

In a previous work [23], the authors demonstrated an alternate approach:
allowing the constitutive equation to be heterogeneous and imposing
significantly weaker requirements on its functional form, such that a unique
deterministic solution can be derived. Unlike other approaches, this approach
is theoretically exact. Furthermore, computational cost is not a major issue,
unlike many alternatives, as the governing equations are linear (it takes on
the order of one second to solve a two-dimensional problem using a laptop
computer). The key assumption made is that the principal directions of the
Cauchy stress $\mathbf{T}$ are the same the observed strain (in the case of
elastic deformation), or strain rate (in the case of plastic deformation).
This leads to the following system of partial differential equations:

$\displaystyle\mathrm{div}{\mathbf{T}}+\mathbf{b}=\rho\ddot{\mathbf{x}},$ (1a)
$\displaystyle[\mathbf{Q^{T}TQ}]_{ij}=0\qquad\mathrm{for}\quad i\neq j.$ (1b)

where $\mathbf{b}$ is the body force, $\ddot{\mathbf{x}}$ is the acceleration,
and $\rho$ is the density. Here, $\mathbf{Q}$ is the matrix of stress
eigenvectors, i.e. the spectral decomposition of $\mathbf{T}$ is
$\mathbf{T}=\mathbf{Q\Lambda_{T}\mathbf{Q}^{T}}$ 111Notation: Bold upper-case
letters to denote second order tensors ($\mathbf{A},\mathbf{B}$,…), bold
lowercase letters to denote vectors ($\mathbf{a},\mathbf{b}$,…), unbolded
letters correspond to scalars ($a$, $b$, $...$), and the superscript $T$
denote the transpose ($\mathbf{A}^{T}$, $\mathbf{a}^{T}$, $...$). We use
standard matrix notation, for example, we have
$\mathbf{a}\cdot\mathbf{b}=\mathbf{a}^{T}\mathbf{b}$ and
$\mathbf{a}\otimes\mathbf{b}=\mathbf{a}\mathbf{b}^{T}$. $\mathbf{T}$ denotes
the Cauchy stress, $\mathbf{b}$ denotes the body force, $\mathbf{x}$ donates
the position of a material element in the deformed frame, and $\rho$ denotes
the density (we introduce other variables as they arise in the paper). We
denote the first and second material time derivative of a second order tensor
$\mathbf{\dot{A}}$ and $\mathbf{\ddot{A}}$ respectively.
$\mathrm{Div}\mathbf{A}$ and $\mathrm{div}\mathbf{A}$ correspond to the
divergence $\mathbf{A}$ with respect to the initial and deformed reference
frame respectively. $|\mathbf{A}|=\sqrt{A_{ij}A_{ij}}$ gives the magnitude of
$\mathbf{A}$ (Einstein summation notation is used). We denote the deviatoric
part of a second order tensor $\mathbf{A}$ as $\mathbf{A_{0}}$, and the trace
as $\mathrm{tr}(\mathbf{A})$. The spectral decomposition of a symmetric tensor
$\mathbf{A}$ is written as
$\mathbf{A}=\mathbf{Q_{\mathbf{A}}\Lambda_{\mathbf{A}}Q_{\mathbf{A}}^{T}}$
where $\mathbf{Q_{\mathbf{A}}}$ has columns with the eigenvectors of
$\mathbf{A}$ denoted $\\{\mathbf{q}^{(i)}|i=1,2,3\\}$ (principal directions)
and $\mathbf{\Lambda_{A}}$ is a diagonal matrix containing the eigenvalues
$\\{\lambda^{i}|i=1,2,3\\}$ (principal values). This decomposition is also
expressed as
$\mathbf{A}=\sum_{i}\lambda_{\mathbf{A}}^{i}\mathbf{q}_{\mathbf{A}}^{(i)}\mathbf{q}_{\mathbf{A}}^{(i)T}$.
As the spectral decomposition is used extensively for the stress, we use the
special notation:
$\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T}$ where
$\\{\sigma^{(i)}|i=1,2,3\\}$ are the principal stresses and
$\\{\mathbf{q}^{(i)}|i=1,2,3\\}$.. Eq. 1a is simply force balance and Eq. 1b
constrains $\mathbf{T}$ to have the principal directions given by
$\mathbf{Q}$. This can be intuitively understood as expressing $\mathbf{T}$ in
the coordinate system aligned with its principal directions and then
specifying the shear components in this coordinate system to be zero. Assuming
$\mathbf{Q}$, $\ddot{\mathbf{x}}$, $\rho$, and $\mathbf{b}$ are known, this
gives rise to a deterministic, linear, variable coefficient second-order,
system of hyperbolic partial differential equations. The core idea is that if
$\mathbf{Q}$ can be determined at every point in space from the measured
strain (or strain rate), these equations can be directly solved for
$\mathbf{T}$. It is then not necessary to apply standard simulation approaches
where material properties must be specified.

It is straight forward to determine $\mathbf{Q}$ in the case of isotropic
(potentially non-linear) elasticity where the stress and strain have the same
principal directions. That is, one can (i) use an approach such as DIC to
measure the full-field strain, (ii) compute the principal directions of the
strain at each point, and (iii) compute the stress by solving Eq. 1. The same
approach can be taken for plastic deformation in the case where the elastic
strain is negligible. Here, one determines $\mathbf{Q}$ from the strain rate
(which can also be computed from the deformation).

One of the limitations of this approach is that it does not directly apply to
materials where the principal directions of the stress do not align with
either the strain or strain rate. This is particularly problematic in many
cases of isotropic elasto-plastic deformation, even when the elastic component
of the deformation is small. Frequently, regions of the sample will
elastically unload while others plastically deform, and it is not be possible
to determine the principal directions of the stress without additional
information. This challenge is illustrated another work currently in
submission [24]. The authors applied Eq. 1 to an experimental data tensile
test on a metallic "dog bone" sample, however, the combination of elastic
unloading and plastic necking prevented direct application of Eq. 1 as one
cannot determine $\mathbf{Q}$. Instead approximate assumptions regarding the
principal directions were required. The topic of this paper will be an
extension of this approach to the case of elasto-plastic deformation where
$\mathbf{Q}$ cannot be directly obtained from the strain or strain rate. We
make the additional assumption that the elastic component of the constitutive
equation is known a priori. Interestingly, the derived system of equations is
non-linear and time dependent.

After deriving a new system of equations, we consider the case study of
necking in a two-dimensional uniaxial tensile test specimen. We use the
derived equations to extract the material stress-strain relationship in the
necking region (at strains on the order of 100 %). We compare the solution
with the solution of the corresponding forward problem (i.e. solving for the
stress and strain when the constitutive equations are known) obtained with a
commercial finite element solver. There is already existing literature on
estimating the stress in a necking tensile specimen. Two geometries have been
of particular interest: cylindrical specimens (for the analyses of thick
plates or other bulk material), or a flat specimens (for the analyses of sheet
materials). The most well known approximation used for cylindrical specimens
is the Bridgman solution which is based on the assumption of constant
equivalent strain in the cross section [25]. However, no such well-established
solution has been developed for flat rectangular specimens [26]. A common
assumption used in this latter case, which we adopt later in the present
paper, is that the material is in a state of plane stress. This assumption
will be valid during the formation of a diffuse neck but will break down when
localized necking occurs [27]. In this regime, the full-field strain can be
measured for flat rectangular specimens using digital image correlation [28]
and various statistical inverse problem solutions have been applied. Iterative
finite element methods minimize the difference in assumed strain and measured
strain (for example, [29]). Alternatively, the virtual field method minimizes
the error in the force balance equation [30]. Other approaches minimize error
in work applied to the neck (e.g. [31]), minimize the error in the force-
displacement curve measured for the sample [32], or estimate and apply
correction factors [26]. However, all of these approaches involve assumptions
about the hardening behavior and do not give exact solutions. Errors vary
depending on the material under investigation, the data available and the
strain level, but are typically lower than more general applications of the
inverse problem as the material properties can be assumed homogeneous
(typically 1-10%). For a more comprehensive review the reader should refer to
[33].

##  2 Theory

###  2.1 Preliminaries

We use standard notation to specify the kinematics of large deformation. We
have a smooth one to one mapping from the reference body $\mathbf{X}$ to the
deformed body is given by $\mathbf{x}=\bar{\mathbf{x}}(\mathbf{X},t)$. This
can be used to specify the deformation gradient

$\mathbf{F}=\mathrm{Grad}(\mathbf{x})$ (2)

and the velocity gradient

$\mathbf{L}=\dot{\mathbf{F}}\mathbf{F^{-1}}.$ (3)

This can be decomposed into symmetric and asymmetric parts

$\mathbf{D}=(\mathbf{L}+\mathbf{L}^{T})/2,\qquad\mathbf{W}=(\mathbf{L}-\mathbf{L}^{T})/2.$
(4)

The left Cauchy-Green tensor is given by

$\mathbf{B}=\mathbf{F}\mathbf{F}^{T}.$ (5)

For elasto-plastic deformation we take the standard Kröner-Lee multiplicative
decomposition

$\mathbf{F}=\mathbf{F}^{e}\mathbf{F}^{p}.$ (6)

Similarly we have

$\mathbf{L}^{e}=\dot{\mathbf{F}}^{e}\mathbf{F}^{e-1},\qquad\mathbf{L}^{p}=\dot{\mathbf{F}}^{p}\mathbf{F}^{p-1},$
(7)
$\mathbf{D}^{e}=(\mathbf{L}^{e}+\mathbf{L}^{eT})/2,\qquad\mathbf{W}^{e}=(\mathbf{L}^{e}-\mathbf{L}^{eT})/2,$
(8)
$\mathbf{D}^{p}=(\mathbf{L}^{p}+\mathbf{L}^{pT})/2,\qquad\mathbf{W}^{p}=(\mathbf{L}^{p}-\mathbf{L}^{pT})/2.$
(9)

In standard theories of isotropic elasto-plasticity we have that $\mathbf{T}$
has the same eigenvectors as the elastic left Cauchy deformation tensor
$\mathbf{B}^{e}=\mathbf{F}^{e}\mathbf{F}^{eT}$. We also have that the stress
(Cauchy stress) has the same principal directions as
$\mathbf{F}^{e}\mathbf{D}^{p}\mathbf{F}^{eT}$.

The key challenge in determining the principal directions of $\mathbf{T}$ in
the case of elasto-plastic deformation is that one can only directly obtain
$\mathbf{B}$ and $\mathbf{D}$ from the observed deformation. $\mathbf{B}^{e}$,
$\mathbf{F}^{e}$ and $\mathbf{D}^{p}$ are unknowns. This is not a problem in
the simplified case of elastic deformation as one can simply take
$\mathbf{B}^{e}=\mathbf{B}$. Similarly in the case of plastic deformation
where the elastic strain is negligible one can take
$\mathbf{F}^{e}\mathbf{D}^{p}\mathbf{F}^{eT}=\mathbf{D}$. Note that this
translates to infinitesimal deformation if one replaces
($\mathbf{B},\mathbf{D},\mathbf{B}^{e},\mathbf{F}^{e}\mathbf{D}^{p}\mathbf{F}^{eT}$)
with the infinitesimal deformation measures
($\mathbf{E},\dot{\mathbf{E}},\mathbf{E}^{e},\dot{\mathbf{E}}^{p}$).

###  2.2 Additional assumption of known elasticity relationship

Allowing the constitutive equation to be elasto-plastic introduces an
additional degree of freedom into the problem. Hence, additional assumptions
are required in order to obtain a deterministic system which can be solved.
This is a general problem in the application of this inverse problem approach:
as one introduces additional degrees of freedom or complexity into the
constitutive equations, one must also introduce additional assumptions to
retain a deterministic system (see [23] for further discussion on this). The
challenge is making assumptions or using models that enable any free
parameters to be obtained from empirical observation.

The primary additional assumption made is that the elastic component of the
constitutive equation is known and invertible, i.e. the function
$\mathbf{B}^{e}=\bar{\mathbf{B}^{e}}(\mathbf{T})$ is known. This will often
the case, for example, many materials are governed by linear elasticity and
more complex plastic deformation relationships. If the material is
homogeneous, the elastic constants can be simply determined using techniques
such as ultrasonic testing. Alternatively, if it is heterogeneous, the elastic
constants can be computed from the initial stage of deformation by solving Eq.
1 (assuming the yield stress is finite) [23]. Once these are determined, one
could calculate the stress throughout the much larger elasto-plastic
deformation by solving the equations developed in this paper. This is true
even in cases where the material is heterogeneous, since the elastic constants
can be obtained by solving the inverse problem in the simpler case of pure
elasticity.

One limitation of this assumption may arise when one is considering a
macroscopic problem and microscopic damage occurs altering the macroscopic
elastic constants [34]. It is likely possible to introduce a model to account
for this, although this would require additional assumptions and material
models which are beyond the present scope. In addition to damage, other
plasticity processes may alter the elastic constants [35], though these are
often not significant.

We note however, that this assumption does not require linear elasticity, and
the approach may find applications in other cases of elasto-plasticity where
the elastic behavior is non-linear and known, while the plastic behavior is
unknown. For example, there may be applications for elastomers [36], or shape
memory alloys [37]. However, it is unclear if there are scenarios where one
can determine the non-linear elastic component of the constitutive equation at
high strains, but where it would still be useful to determine the plastic
component using the proposed approach.

###  2.3 Derivation of governing equations

Here, we derive the elasto-plastic governing equations for the general case of
finite deformation, however, we also present the simpler derivation for the
infinitesimal case in the Appendix. The reader may wish to consult this first
as the main lines of argument are simpler. In order to derive the elasto-
plastic governing equations, we consider the deformation of a material element
from time $t=\tau$ to $t=\tau+\Delta t$. It is assumed that the Cauchy stress
$\mathbf{T}$ is known at time $\tau$ and we wish to determine the stress at
time $\tau+\Delta t$ based on the observed deformation (Fig. 1). We introduce
subscripts to denote quantities evaluated at these times, e.g.
$\mathbf{T}_{\tau}=\bar{\mathbf{T}}(\tau)$ and $\mathbf{T}_{\tau+\Delta
t}=\bar{\mathbf{T}}(\tau+\Delta t)$. This increment of deformation is
characterized by the deformation gradient $\mathbf{F}^{*}$:

Figure 1: The references frames and deformation gradients used to map between
them. The sample is initially deformed with
$\mathbf{F}_{\tau}=\mathbf{F}_{\tau}^{e}\mathbf{F}_{\tau}^{p}$, then an
additional deformation $\mathbf{F}^{*}$ is observed which is related to the
unobserved $\mathbf{F}^{p*}$ and $\mathbf{F}^{e}_{\tau+\Delta t}$. The triple
red lines indicate that the principal directions can be directly related to
one another. This is analogous to being parallel as the sum of two tensors
with the same principal directions also have the same principal directions.

$\mathbf{F}_{\tau+\Delta t}=\mathbf{F}^{*}\mathbf{F}_{\tau}.$ (10)

We define $\mathbf{F}^{p*}$ in the similar way:

$\mathbf{F}^{p}_{\tau+\Delta t}=\mathbf{F}^{p*}\mathbf{F}^{p}_{\tau}$ (11)

hence, we have:

$\mathbf{F}^{*}\mathbf{F}^{e}_{\tau}\mathbf{F}^{p}_{\tau}=\mathbf{F}^{e}_{\tau+\Delta
t}\mathbf{F}^{p*}\mathbf{F}^{p}_{\tau},$ (12)
$\mathbf{F}^{*}\mathbf{F}^{e}_{t}=\mathbf{F}^{e}_{\tau+\Delta
t}\mathbf{F}^{p*}.$ (13)

This is the relationship depicted in the triangle in Fig. 1. We express the
spectral decomposition of $\mathbf{T}$ as

$\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},$ (14)

where $\\{\mathbf{q}^{(i)}|i\in 1,2,3\\}$ are the eigenvectors and
$\\{\sigma^{(i)}|i\in 1,2,3\\}$ are the eigenvalues such that:

$\mathbf{T}\mathbf{q}^{(i)}=\sigma^{(i)}\mathbf{q}^{(i)}.$ (15)

We also assume the eigenvalues do not repeat (see Section 2.5 for discussion).
Due to the assumption that the eigenvectors of the elastic left Cauchy
deformation tensor are the same as $\mathbf{T}$ (see Section 1), we have:

$\mathbf{B}^{e}\mathbf{q}^{(i)}=\mathbf{F}^{e}\mathbf{F}^{eT}\mathbf{q}^{(i)}=\lambda^{(i)}\mathbf{q}^{(i)}.$
(16)

Here, we use the convention that $\\{\mathbf{q}^{(i)}\\}$ always corresponds
to the eigenvector of the stress, however, $\\{\lambda^{(i)}\\}$ are just
generic eigenvalues. If we wish to refer to the eigenvectors and eigenvalues
of some other specific tensor $\mathbf{A}$ we will add a subscript so we have
$\\{\mathbf{q}_{\mathbf{A}}^{(i)}\\}$ and $\\{\lambda_{\mathbf{A}}^{(i)}\\}$.
We also take the eigenvectors to be normalized. Hence, we have
$\mathbf{q}^{(i)T}\mathbf{q}^{(j)}=\delta_{ij}$ where $\delta_{ij}$ is the
Kronecker delta that takes the value 1 if $i=j$ and 0 otherwise. As the
principal directions of stress are aligned with
$\mathbf{F}^{e}\mathbf{D}^{p}\mathbf{F}^{eT}$ we have

$\mathbf{F}^{e}\mathbf{D}^{p}\mathbf{F}^{eT}\mathbf{q}^{(i)}=\lambda^{(i)}\mathbf{q}^{(i)}.$
(17)

where the $\\{\lambda^{(i)}\\}$ are different from Eq. 16 but
$\\{\mathbf{q}^{(i)}\\}$ are still the eigenvectors of the stress. Consider
the following expression:

$\mathbf{F}^{e}_{\tau+\Delta t}\mathbf{F}^{eT}_{\tau+\Delta t}+2\Delta
t\mathbf{F}^{e}_{\tau+\Delta t}\mathbf{D}^{p}_{\tau+\Delta
t}\mathbf{F}^{eT}_{\tau+\Delta t}.$ (18)

This will have the same eigenvectors as $\mathbf{T}_{\tau+\Delta t}$ because
the sum of two tensors with the same eigenvectors also have the same
eigenvectors and multiplication by a scalar preserves the eigenvectors. Hence
we have:

$\mathopen{}\mathclose{{}\left(\mathbf{F}^{e}_{\tau+\Delta
t}\mathbf{F}^{eT}_{\tau+\Delta t}+2\Delta t\mathbf{F}^{e}_{\tau+\Delta
t}\mathbf{D}^{p}_{\tau+\Delta t}\mathbf{F}^{eT}_{\tau+\Delta
t}}\right)\mathbf{q}^{(i)}_{\tau+\Delta
t}=\lambda^{(i)}\mathbf{q}^{(i)}_{\tau+\Delta t}.$ (19)

We manipulate this expression in order to find an expression for the
eigenvectors in terms of observed quantities. We have:

$\mathopen{}\mathclose{{}\left(\mathbf{F}^{e}_{\tau+\Delta
t}(\mathbf{I}+2\Delta t\mathbf{D}^{p}_{\tau+\Delta
t})\mathbf{F}^{eT}_{\tau+\Delta t}}\right)\mathbf{q}^{(i)}_{\tau+\Delta
t}=\lambda^{(i)}\mathbf{q}^{(i)}_{\tau+\Delta t}.$ (20)
$\mathopen{}\mathclose{{}\left(\mathbf{F}^{e}_{\tau+\Delta
t}(\mathbf{I}+\Delta t\mathbf{L}^{p}_{\tau+\Delta t}+\Delta
t\mathbf{L}^{pT}_{\tau+\Delta t})\mathbf{F}^{eT}_{\tau+\Delta
t}}\right)\mathbf{q}^{(i)}_{\tau+\Delta
t}=\lambda^{(i)}\mathbf{q}^{(i)}_{\tau+\Delta t}.$ (21)

where,

$\displaystyle\mathbf{I}+\Delta t\mathbf{L}^{p}_{\tau+\Delta t}+\Delta
t\mathbf{L}^{pT}_{\tau+\Delta t}$ $\displaystyle=(\mathbf{I}+{\Delta
t}\mathbf{L}^{p}_{\tau+\Delta t})(\mathbf{I}+{\Delta
t}\mathbf{L}^{p}_{\tau+\Delta t})^{T}-{\Delta
t}^{2}\mathbf{L}^{p}_{\tau+\Delta t}\mathbf{L}^{pT}_{\tau+\Delta t}$ (22)
$\displaystyle=(\mathbf{I}+{\Delta t}\mathbf{L}^{p}_{\tau+\Delta
t})(\mathbf{I}+{\Delta t}\mathbf{L}^{p}_{\tau+\Delta t})^{T}+O(\Delta t^{2})$
$\displaystyle=\mathbf{F}^{p*}\mathbf{F}^{p*T}+O(\Delta t^{2}).$

This last expression follows from

$\displaystyle\mathbf{F}^{p*}$ $\displaystyle=\mathbf{F}^{p}_{\tau+\Delta
t}\mathbf{F}^{p-1}_{\tau},$ (23)
$\displaystyle=\mathopen{}\mathclose{{}\left(\mathbf{F}^{p}_{\tau}+\dot{\mathbf{F}}_{\tau}^{p}\Delta
t+O({\Delta t}^{2})}\right)\mathbf{F}^{p-1}_{\tau},$
$\displaystyle=\mathbf{I}+\dot{\mathbf{F}}_{\tau}^{p}\mathbf{F}^{p-1}_{\tau}\Delta
t+O({\Delta t}^{2}),$ $\displaystyle=\mathbf{I}+\mathbf{L}^{p}_{\tau}\Delta
t+O({\Delta t}^{2}),$ $\displaystyle=\mathbf{I}+\mathbf{L}^{p}_{\tau+\Delta
t}\Delta t+O({\Delta t}^{2}).$

Substituting Eq. 22 into Eq. 21 gives:

$\mathopen{}\mathclose{{}\left(\mathbf{F}^{e}_{\tau+\Delta
t}\mathbf{F}^{p*}\mathbf{F}^{p*T}\mathbf{F}^{eT}_{\tau+\Delta
t}}\right)\mathbf{q}^{(i)}_{\tau+\Delta t}+O(\Delta
t^{2})=\lambda^{(i)}\mathbf{q}^{(i)}_{\tau+\Delta t}.$ (24)

From Eq. 12 we have $\mathbf{F}^{e}_{\tau+\Delta
t}\mathbf{F}^{p*}=\mathbf{F}^{*}\mathbf{F}^{e}_{\tau}$, hence:

$\mathopen{}\mathclose{{}\left(\mathbf{F}^{*}\mathbf{F}^{e}_{\tau}\mathbf{F}^{eT}_{\tau}\mathbf{F}^{*T}}\right)\mathbf{q}^{(i)}_{\tau+\Delta
t}+O(\Delta t^{2})=\lambda^{(i)}\mathbf{q}^{(i)}_{\tau+\Delta t}.$ (25)

This shows that when ${\Delta t}$ is small we can determine the eigenvectors
for $\mathbf{T}_{\tau+\Delta t}$ from the observed $\mathbf{F}^{*}$ and the
quantity $\mathbf{F}^{e}_{\tau}\mathbf{F}^{eT}_{\tau}$ which known from
$\mathbf{T}_{\tau}$ and the function $\bar{\mathbf{B}^{e}}(\mathbf{T})$. We
can use this expression to derive a differential equation for the evolution of
$\\{\mathbf{q}^{(i)}\\}$.

For some arbitrary tensor $\mathbf{A}$ with eigenvectors
$\\{\mathbf{q}_{\mathbf{A}}^{(i)}\\}$ and eigenvalues
$\\{\lambda_{\mathbf{A}}^{(i)}\\}$ that varies with time,
$\mathbf{A}=\bar{\mathbf{A}}(t)$, we have the following general expression for
the derivative of its eigenvectors:

$\dot{\mathbf{q}}^{(i)}_{\mathbf{A}}=\sum_{j,j\neq
i}\frac{\mathbf{q}_{\mathbf{A}}^{(j)}\mathbf{q}_{\mathbf{A}}^{(j)T}\dot{\mathbf{A}}\mathbf{q}_{\mathbf{A}}^{(i)}}{\lambda^{(i)}_{\mathbf{A}}-\lambda^{(j)}_{\mathbf{A}}},$
(26)

where this assumes there are no repeated eigenvalues and ensures that
$\\{\mathbf{q}_{\mathbf{A}}^{(i)}\\}$ remain normalized.

Let

$\mathbf{A}=\mathbf{F}^{*}\mathbf{F}^{e}_{\tau}\mathbf{F}^{eT}_{\tau}\mathbf{F}^{*T}.$
(27)

This will have the same eigenvectors as the stress up to error of $O(\Delta
t^{2})$ (Eq. 25). Hence, in the limit of $\Delta t\rightarrow 0$, we have
$\mathbf{q}^{(i)}_{\mathbf{A}}=\mathbf{q}^{(i)}$ and
$\dot{\mathbf{q}}^{(i)}_{\mathbf{A}}=\dot{\mathbf{q}}^{(i)}$. We will
determine an expression for $\dot{\mathbf{A}}$ and
$\\{\lambda_{\mathbf{A}}^{(i)}\\}$ to derive an expression for
$\dot{\mathbf{q}}^{(i)}$. Technically, $\mathbf{A}$ will be a function $\tau$
and $\tau+{\Delta t}$ as $\mathbf{F}^{e}_{\tau}$ is evaluated at $\tau$ and
$\mathbf{F}^{*}$ depends on both $\tau$ and $\tau+{\Delta t}$. However,
because $\mathbf{A}$ is used to determine $\\{\mathbf{q}^{(i)}\\}$ at times
just after $\tau$, we consider $\tau$ fixed and allow $t=\tau+\Delta t$ to
vary. Hence we have $\mathbf{A}=\bar{\mathbf{A}}(t)$ and
$\dot{\mathbf{A}}=d\mathbf{A}/dt$. We then have

$\dot{\mathbf{A}}=\lim_{\Delta t\rightarrow
0}\frac{\bar{\mathbf{A}}(t=\tau+\Delta t)-\bar{\mathbf{A}}(t=\tau)}{\Delta t}$
(28)

$\dot{\mathbf{A}}=\lim_{\Delta t\rightarrow
0}\frac{(\mathbf{I}+\mathbf{L}_{\tau}\Delta
t)\mathbf{F}^{e}_{\tau}\mathbf{F}^{eT}_{\tau}(\mathbf{I}+\mathbf{L}_{\tau}^{T}\Delta
t)-\mathbf{F}^{e}_{\tau}\mathbf{F}^{eT}_{\tau}+O(\Delta t^{2})}{\Delta t}$
(29)

where this substitution arises from a similar argument to that made in Eq. 23.
Substituting for $\mathbf{B}^{e}$ and dropping the subscripts $\tau$ we have:

$\dot{\mathbf{A}}=\lim_{\Delta t\rightarrow
0}\frac{(\mathbf{L}\mathbf{B}^{e}+\mathbf{B}^{e}\mathbf{L}^{T})\Delta
t+O(\Delta t^{2})}{\Delta t}$ (30)

$\dot{\mathbf{A}}=\mathbf{L}\mathbf{B}^{e}+\mathbf{B}^{e}\mathbf{L}^{T}.$ (31)

Evaluating $\\{\lambda_{\mathbf{A}}^{(i)}\\}$ at time $t=\tau$ gives
$\lambda_{\mathbf{A}}^{(i)}=\lambda^{e(i)}_{\mathbf{B}}$ where
$\\{\lambda_{\mathbf{B}}^{e(i)}\\}$ are the eigenvalues of $\mathbf{B}^{e}$.
Substituting these quantities into Eq. 26 gives:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\mathopen{}\mathclose{{}\left(\mathbf{L}\mathbf{B}^{e}+\mathbf{B}^{e}\mathbf{L}^{T}}\right)\mathbf{q}^{(i)}}{\lambda^{e(i)}_{\mathbf{B}}-\lambda^{e(j)}_{\mathbf{B}}},$
(32)

When combined with force balance, the assumption that $\mathbf{B}^{e}$ can be
determined by $\mathbf{T}$, and that $\mathbf{L}$ can be directly obtained
from the observed deformation, this is sufficient to give a complete system of
equations and unknowns.

We manipulate this expression further to give an alternate equation in terms
of the velocity gradient $\mathbf{D}$ and spin $\mathbf{W}$. Using Eq. 16 we
have:

$\displaystyle\mathbf{q}^{(j)T}(\mathbf{L}\mathbf{B}^{e}_{t}+\mathbf{B}^{e}_{t}\mathbf{L}^{T})\mathbf{q^{(i)}}$
$\displaystyle=\lambda^{e(i)}_{\mathbf{B}}\mathbf{q}^{(j)T}\mathbf{L}\mathbf{q}^{(i)}+\lambda^{e(j)}_{\mathbf{B}}\mathbf{q}^{(j)T}\mathbf{L}^{T}\mathbf{q}^{(i)}$
(33)
$\displaystyle=\mathbf{q}^{(j)T}(\lambda^{e(i)}_{\mathbf{B}}\mathbf{L}+\lambda^{e(j)}_{\mathbf{B}}\mathbf{L}^{T})\mathbf{q}^{(i)}$
$\displaystyle=\frac{(\lambda^{e(i)}_{\mathbf{B}}+\lambda^{e(j)}_{\mathbf{B}})}{2}\mathbf{q}^{(j)T}(\mathbf{L}+\mathbf{L}^{T})\mathbf{q}^{(i)}+\frac{(\lambda^{e(i)}_{\mathbf{B}}-\lambda^{e(j)}_{\mathbf{B}})}{2}\mathbf{q}^{(j)T}(\mathbf{L}-\mathbf{L}^{T})\mathbf{q}^{(i)}$
$\displaystyle=(\lambda^{e(i)}_{\mathbf{B}}+\lambda^{e(j)}_{\mathbf{B}})\mathbf{q}^{(j)T}\mathbf{D}\mathbf{q}^{(i)}+(\lambda^{e(j)}_{\mathbf{B}}-\lambda^{e(i)}_{\mathbf{B}})\mathbf{q}^{(j)T}\mathbf{W}\mathbf{q}^{(i)}$

Substituting back into Eq. 32 gives:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\mathbf{D}\mathbf{q}^{(i)}(\lambda^{e(j)}_{\mathbf{B}}+\lambda^{e(i)}_{\mathbf{B}})}{\lambda^{e(j)}_{\mathbf{B}}-\lambda^{e(i)}_{\mathbf{B}}}+\sum_{j,j\neq
i}\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\mathbf{W}\mathbf{q}^{(i)}$ (34)

Note that for $i=j$, $\mathbf{q}^{(j)T}\mathbf{W}\mathbf{q}^{(i)}=0$ as
$\mathbf{W}$ is antisymmetric. Hence, we can replace the second sum in Eq. 34
with a sum over all $j$. Then the term
$\sum_{j}\mathbf{q}^{(j)}\mathbf{q}^{(j)T}$ becomes the identity $\mathbf{I}$.
This gives:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\mathbf{D}\mathbf{q}^{(i)}(\lambda^{e(j)}_{\mathbf{B}}+\lambda^{e(i)}_{\mathbf{B}})}{\lambda^{e(j)}_{\mathbf{B}}-\lambda^{e(i)}_{\mathbf{B}}}+\mathbf{W}\mathbf{q}^{(i)}$
(35)

To summarize these results, we now have a complete set of partial differential
equations

$\displaystyle\mathrm{div}{\mathbf{T}}+\mathbf{b}=\rho\ddot{\mathbf{x}},$
(36a) $\displaystyle\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\mathbf{D}\mathbf{q}^{(i)}(\lambda^{e(j)}_{\mathbf{B}}+\lambda^{e(i)}_{\mathbf{B}})}{\lambda^{e(j)}_{\mathbf{B}}-\lambda^{e(i)}_{\mathbf{B}}}+\mathbf{W}\mathbf{q}^{(i)},$
(36b)
$\displaystyle\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},\qquad\mathbf{B}^{e}=\sum_{i}\lambda^{e(i)}_{\mathbf{B}}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},$
(36c) $\displaystyle\mathbf{B}^{e}=\bar{\mathbf{B}^{e}}(\mathbf{T})$ (36d)

There are two major difference between this system of equations and Eq. 1.
First, this system of equations is non-linear (due to Eq. 35). Second, this
system of equations is time dependent.

In the case of linear elasticity where the elastic constants are known, Eq.
36d can be replaced with

$\mathbf{B}^{e}=\frac{1}{2\mu}\mathopen{}\mathclose{{}\left(\mathbf{T}-\frac{\lambda}{2\mu+3\lambda}\mathrm{tr}(\mathbf{T})\mathbf{I}}\right).$
(37)

###  2.4 Limiting cases

We consider several special cases and show how we arrive at the equations for
small deformation, purely elastic, and purely visco-plastic deformation.

Consider first the case of infinitesimal deformation. Here,
$\mathbf{D}\approx\dot{\mathbf{E}}$. We also have that
$\lambda_{\mathbf{B}}^{(i)e}$ will be close to 1, we express this as
$\lambda_{\mathbf{B}}^{e(i)}=1+O(\epsilon)$ where $\epsilon<<1$. Eq. 35
becomes:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\dot{\mathbf{E}}\mathbf{q}^{(i)}(2+O(\epsilon))}{\lambda^{e(i)}_{\mathbf{B}}-\lambda^{e(j)}_{\mathbf{B}}}+\mathbf{W}\mathbf{q}^{(i)}$
(38)

We neglect the $O(\epsilon)$ term in the numerator as it is small compared to
2. We note that
$\lambda^{e(i)}_{\mathbf{B}}=1+2\lambda^{e(i)}_{\mathbf{E}}+O(\epsilon^{2})$
where $\lambda^{e(i)}_{\mathbf{E}}$ is $O(\epsilon)$. As the denominator of
this sum expression will be $O(\epsilon)$, the second
$\mathbf{W}\mathbf{q}^{(i)}$ term will be negligible in comparison. This gives
the following small strain formulation:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\dot{\mathbf{E}}\mathbf{q}^{(i)}}{\lambda^{e(i)}_{\mathbf{E}}-\lambda^{e(j)}_{\mathbf{E}}}.$
(39)

Note this is directly derived in the Appendix using an additive strain
decomposition.

We next consider a case for a purely elastic deformation. Obviously here the
system could be trivially solved using
$\mathbf{T}=\bar{\mathbf{T}}(\mathbf{B}^{e})$, however, we consider how the
other equations in the system become equivalent to Eq. 1 in the elastic limit,
so the equations can be solved without a known elasticity equation. In the
case $\mathbf{B}^{e}=\mathbf{B}$, Eq. 32 becomes:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\mathopen{}\mathclose{{}\left(\mathbf{L}\mathbf{B}+\mathbf{B}\mathbf{L}^{T}}\right)\mathbf{q}^{(i)}}{\lambda^{e(i)}_{\mathbf{B}}-\lambda^{e(j)}_{\mathbf{B}}}.$
(40)

Note the following relationship

$\dot{\mathbf{B}}=\mathbf{L}\mathbf{B}+\mathbf{B}\mathbf{L}^{T}.$ (41)

We now have

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\dot{\mathbf{B}}\mathbf{q}^{(i)}}{\lambda^{(i)}_{\mathbf{B}}-\lambda^{(j)}_{\mathbf{B}}}$
(42)

This is simply the same as the expression for the rate of change of the
eigenvector of $\mathbf{B}$ (see Eq. 26). Hence, choosing $\mathbf{q}$ to be
an eigenvector of $\mathbf{B}$ will solve the equation and we can replace Eq.
35 with

$\mathbf{B}\mathbf{q}^{(i)}=\lambda^{(i)}{\mathbf{B}}\mathbf{q}^{(i)}.$ (43)

This directly reproduces the result in the elastic case.

Last, we consider the the visco-plastic case. We take the limit of Eq. 35 as
$\mathbf{F}^{e}\rightarrow\mathbf{I}$ indicating no elastic strain. Here all
$\lambda^{e}_{\mathbf{B}i}\rightarrow 1$ and the denominator of the first term
on the right hand side approaches zero. Hence, in order to satisfy the
equation in the limiting case, the numerator must approach 0. Hence, we have

$\mathbf{q}^{(j)T}\mathbf{D}\mathbf{q}^{(i)}=0\qquad\textrm{for all}\quad
i\neq j$ (44)

Since the $\mathbf{q}^{(i)}$ vectors are orthogonal, i.e.
$\mathbf{q}^{(i)T}\mathbf{q}^{(j)}=0$ for all $i\neq j$, one can see that
$\mathbf{q}^{(i)}$ are the eigenvectors of $\mathbf{D}$. This directly
reproduces the result for the visco-plastic case. Hence, we can see that in
Eq. 36 is the most general form of governing equations as the other forms,
presented in Section 1, can be derived as limiting cases of it.

###  2.5 Repeated principal values, initial conditions, and boundary
conditions

In various deformations, multiple principal stresses may take on the same
value so Eq. 36 cannot be directly applied. Let $\\{\sigma^{(k)}\\}$ denote
the set of principal stresses that are repeated over some region of
$(\mathbf{x},t)$. The corresponding principal directions,
$\\{\mathbf{q}^{(k)}\\}$, are not unique and any linear combination of
$\\{\mathbf{q}^{(k)}\\}$ is also a principal direction. Note that this also
corresponds to the repeated principal values of $\mathbf{B}^{e}$,
$\\{\lambda_{\mathbf{B}}^{e(k)}\\}$, as the function
$\bar{\mathbf{B}^{e}}(\mathbf{T})$ is isotropic.

The general case is non-trivial as closed characteristic loops may arise
surrounding the region of repeated eigenvalues. Hence, we do not fully address
this here and it may be the subject of future work. Nevertheless, we consider
the special case of stress free initial conditions, which is relevant to a
large number of deformations.

At time $t=0$, for stress free initial conditions, the choice of
$\\{\mathbf{q}^{(i)}\\}$ will be non-unique. To resolve this ambiguity we
additionally require that each $\mathbf{q}^{(i)}$ is continuously varying and
differentiable with respect to time. This restriction gives rise to unique
$\\{\mathbf{q}_{t0}^{(i)}\\}$ [38, 39]. Here the subscript $t0$ is used to
denote variables at time $t=0$. In the case of stress free initial conditions,
one can simply specify that $\\{\mathbf{q}_{t0}^{(i)}\\}$ is equal to the
principal directions at some infinitesimal later time $t=\Delta t$,
$\\{\mathbf{q}_{\Delta t}^{(i)}\\}$. We can compute this by requiring that
$\\{\mathbf{q}_{\Delta t}^{(i)}\\}$ is consistent with $\mathbf{B}^{e}_{\Delta
t}$ and $\mathbf{F}_{t0}^{e}\mathbf{D}_{t0}^{p}\mathbf{F}_{t0}^{eT}$.

First, we have that the principal directions $\\{\mathbf{q}^{(i)}_{t0}\\}$
will be consistent with $\mathbf{B}^{e}_{\Delta t}$. The Taylor expansion for
$\mathbf{B}^{e}_{\Delta t}$ is given by:

$\mathbf{B}^{e}_{\Delta
t}=\mathbf{B}^{e}_{t0}+\dot{\mathbf{B}}^{e}_{t0}{\Delta t}+O({\Delta t}^{2}).$
(45)

However, as the material is unstressed at time $t=0$, we have that
$\mathbf{B}_{t0}=\mathbf{I}$. Therefore $\\{\mathbf{q}^{(i)}_{t0}\\}$ will be
eigenvectors of $\dot{\mathbf{B}}^{e}_{t0}$. We then have:

$\displaystyle\dot{\mathbf{B}}_{t0}^{e}$
$\displaystyle=\mathbf{L}_{t0}^{e}\mathbf{B}_{t0}^{e}+\mathbf{B}_{t0}^{e}\mathbf{L}_{t0}^{eT},$
(46) $\displaystyle=\mathbf{L}_{t0}^{e}+\mathbf{L}_{t0}^{eT}.$

So $\\{\mathbf{q}^{(i)}_{t0}\\}$ will be consistent with
$\mathbf{L}_{t0}^{e}+\mathbf{L}_{t0}^{eT}$.

Second, we have that $\\{\mathbf{q}^{(i)}_{t0}\\}$ will be consistent with
$\mathbf{F}_{t0}^{e}\mathbf{D}_{t0}^{p}\mathbf{F}_{t0}^{eT}$. However, as
$\mathbf{F}_{t0}^{e}=\mathbf{I}$, the stress will be aligned with
$\mathbf{D}^{p}$. Due to the definition of $\mathbf{D}^{p}$,
$\\{\mathbf{q}^{(i)}_{t0}\\}$ will also be consistent with
$\mathbf{L}_{t0}^{p}+\mathbf{L}_{t0}^{pT}$.

As $\\{\mathbf{q}^{(i)}_{t0}\\}$ is consistent with both
$\mathbf{L}_{t0}^{e}+\mathbf{L}_{t0}^{eT}$ and
$\mathbf{L}_{t0}^{p}+\mathbf{L}_{t0}^{pT}$, we have that
$\\{\mathbf{q}^{(i)}_{t0}\\}$ must be consistent with
$\mathbf{D}_{t0}=\mathbf{L}_{t0}+\mathbf{L}_{t0}^{T}$ given:

$\displaystyle\mathbf{L}_{t0}$
$\displaystyle=\mathbf{L}_{t0}^{e}+\mathbf{F}_{t0}^{e}\mathbf{L}_{t0}^{p}\mathbf{F}_{t0}^{e-1},$
(47) $\displaystyle=\mathbf{L}_{t0}^{e}+\mathbf{L}_{t0}^{p}.$

Hence we have

$\mathbf{D}_{t0}\mathbf{q}_{t0}^{(i)}=\lambda^{(i)}\mathbf{q}_{t0}^{(i)}.$
(48)

Therefore, the initial conditions for $\\{\mathbf{q}^{(i)}\\}$ can be obtained
from the observed quantity $\mathbf{D}_{t0}$.

To summarize the initial conditions, if the material is initially in a stress-
free state, one can specify the initial stress to be zero, and the principal
directions to be aligned with $\mathbf{D}$. If the material is not initially
stress free, the stress must be known. In this case one can directly obtain
the principal directions from the known stress.

For a detailed discussion on the boundary conditions, the reader should refer
to [24]. The number of constraints on the stress imposed at the boundary,
$n^{e}$, directly corresponds to the number of characteristic lines entering
the domain. If one has knowledge of the traction vector at the boundary
$\mathbf{t}$, one can use the equations explained in [24] to specify these
components.

$\sum_{i}\sigma^{(i)}(\mathbf{d}^{(j)*T}\mathbf{q}^{(i)})(\mathbf{q}^{(i)T}\mathbf{n})=\mathbf{d}^{(j)*T}\mathbf{t}.$
(49)

where $\mathbf{n}$ is the domain boundary normal, $\mathbf{t}$ is the known
traction, and $\\{\mathbf{d}^{(j)*}\\}$ is a subset of size $n^{e}$ of an
arbitrarily chosen orthonormal basis $\\{\mathbf{d}^{(i)}\\}$. These $n^{e}$
equations can be used to solve for the unknown principal stresses at the
boundary.

##  3 Linearization

As the governing system of equations is non-linear, it is useful to derive
their linearized form about some equilibrium solution which may be useful in a
number of analytical or computational settings (e.g. [40]). We consider the
linearization of Eq. 36, about an equilibrium solution at time $t_{[0]}$
during an arbitrary deformation. For simplicity we assume quasi-static
deformation with no body force, i.e. in the force balance equation we take
$\ddot{\mathbf{x}}=0$ and $\mathbf{b}=0$. Where needed, incorporation of these
terms into the linearization can be carried out without significant
challenges.

We introduce a small quantity $\epsilon$, where $\epsilon<<1$, to capture
small changes from this equilibrium solution. We introduce subscripts $[0]$ to
denote equilibrium quantities and subscripts $[1]$ to denote deviations away
from equilibrium. I.e. we have

$\displaystyle t$ $\displaystyle=t_{[0]}+\epsilon t_{[1]},$ (50a)
$\displaystyle\mathbf{q}^{(i)}$
$\displaystyle=\mathbf{q}^{(i)}_{[0]}+\epsilon\mathbf{q}^{(i)}_{[1]}+O(\epsilon^{2})=\mathbf{q}^{(i)}_{[0]}+\epsilon\dot{\mathbf{q}}^{(i)}t_{[1]}+O(\epsilon^{2}),$
(50b) $\displaystyle\sigma^{(i)}$
$\displaystyle=\sigma^{(i)}_{[0]}+\epsilon\sigma^{(i)}_{[1]}+O(\epsilon^{2}).$
(50c)

In general we treat $\\{\sigma^{(i)}\\}$ and $\\{\mathbf{q}^{(i)}\\}$ as our
independent variables, though we will often refer to $\mathbf{T}$ and
$\mathbf{T}_{[0]}$ where convenient with

$\mathbf{T}_{[0]}=\sum_{i}\sigma_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)T}.$
(51)

Furthermore, we simply refer to $\dot{\mathbf{F}}$,
$\\{\dot{\mathbf{q}}^{(i)}\\}$ and $\mathbf{L}$ with no subscripts as second
derivatives will be small in comparison to other terms. We have

$\mathbf{F}=\mathbf{F}_{[0]}+\epsilon\dot{\mathbf{F}}t_{[1]}+O(\epsilon^{2}).$
(52)

We also use $\mathbf{x}_{[0]}$, $\Omega_{[0]}$ and $\partial\Omega_{[0]}$ to
denote the positions of material elements, the domain, and domain boundary at
time $t_{[0]}$.

First, we consider the equilibrium solution at some arbitrary stage in a
deformation, where, by definition, the quantities are not changing with time.
We assume that $\\{\mathbf{q}^{(i)}_{[0]}\\}$ is known. The quantities at time
$t_{[0]}$ will satisfy

$\mathrm{div}_{[0]}(\mathbf{T}_{[0]})=0,\qquad\mathrm{in}\qquad\Omega_{[0]},$
(53)
$\mathbf{T}\mathbf{n}=\mathbf{t}_{[0]},\qquad\mathrm{on}\qquad\partial\Omega_{[0]},$
(54)

where $\mathrm{div}_{[0]}$ is the divergence with respect $\mathbf{x}_{[0]}$,
i.e. $\mathrm{div}_{[0]}(\mathbf{A})=\partial A_{ij}/\partial x_{[0]j}$.

Next we consider the linearized form of Eq. 36b. Substituting Eq. 32 into Eq.
50b gives:

$\mathbf{q}^{(i)}=\mathbf{q}^{(i)}_{[0]}+\epsilon\mathbf{t}_{[1]}\sum_{j,j\neq
i}\frac{\mathbf{q}_{[0]}^{(j)}\mathbf{q}^{(j)T}_{[0]}\mathopen{}\mathclose{{}\left(\mathbf{L}\mathbf{B}^{e}_{[0]}+\mathbf{B}^{e}_{[0]}\mathbf{L}^{T}}\right)\mathbf{q}_{[0]}^{(i)}}{\lambda^{e(i)}_{\mathbf{B}[0]}-\lambda^{e(j)}_{\mathbf{B}[0]}}+O(\epsilon^{2}),$
(55)

where $\mathbf{B}^{e}_{[0]}=\mathbf{F}_{[0]}^{e}\mathbf{F}_{[0]}^{eT}$ and
$\\{\lambda^{e(i)}_{\mathbf{B}[0]}\\}$ are the corresponding eigenvalues. We
also have:

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}_{[0]}^{(j)}\mathbf{q}^{(j)T}_{[0]}\mathopen{}\mathclose{{}\left(\mathbf{L}\mathbf{B}^{e}_{[0]}+\mathbf{B}^{e}_{[0]}\mathbf{L}^{T}}\right)\mathbf{q}_{[0]}^{(i)}}{\lambda^{e(i)}_{\mathbf{B}[0]}-\lambda^{e(j)}_{\mathbf{B}[0]}}+O(\epsilon).$
(56)

which gives the linearized form of Eq. 36b.

Prior to deriving the linearized form of the force balance equation it is
necessary to consider several relationships regarding the deformation
gradient. We have:

$\mathbf{F}=\mathbf{F}_{[0]}+\epsilon
t_{[1]}\mathbf{L}\mathbf{F}_{[0]}+O(\epsilon^{2}).$ (57)

We also define $\mathbf{F}^{*}=\partial\mathbf{x}/\partial\mathbf{x}_{[0]}$.
We have:

$\displaystyle\mathbf{F}^{*}=\mathbf{F}\mathbf{F}^{-1}_{[0]},$ (58a)
$\displaystyle\mathbf{F}^{*}=\mathbf{I}+\epsilon
t_{[1]}\mathbf{L}+O(\epsilon^{2})$ (58b)
$\displaystyle\mathbf{F}^{*-1}=\mathbf{I}-\epsilon
t_{[1]}\mathbf{L}+O(\epsilon^{2}),$ (58c)
$\displaystyle\mathrm{det}(\mathbf{F}^{*})=1+\epsilon
t_{[1]}\mathrm{tr}(\mathbf{L})+O(\epsilon^{2})$ (58d)

The referential form of the force balance equation gives

$0=\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\mathrm{det}(\mathbf{F}^{*})\mathbf{T}\mathbf{F}^{*-T}}\right),$
(59)
$0=\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\mathrm{det}(\mathbf{F}^{*})\mathopen{}\mathclose{{}\left(\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T}}\right)\mathbf{F}^{*-T}}\right),$
(60)
$0=\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(1+\epsilon
t_{[1]}\mathrm{tr}(\mathbf{L})}\right)\mathopen{}\mathclose{{}\left(\sum_{i}\mathopen{}\mathclose{{}\left(\sigma_{[0]}^{(i)}+\epsilon\sigma_{[1]}^{(i)}}\right)\mathopen{}\mathclose{{}\left(\mathbf{q}^{(i)}_{[0]}+\epsilon\mathbf{q}^{(i)}_{[1]}}\right)\mathopen{}\mathclose{{}\left(\mathbf{q}^{(i)T}_{[0]}+\epsilon\mathbf{q}^{(i)T}_{[1]}}\right)}\right)\mathopen{}\mathclose{{}\left(1-\epsilon
t_{[1]}\mathbf{L}^{T}}\right)}\right)+O(\epsilon^{2}).$ (61)

Consider the summation component of this equation:

$\sum_{i}\mathopen{}\mathclose{{}\left(\sigma_{[0]}^{(i)}+\epsilon\sigma_{[1]}^{(i)}}\right)\mathopen{}\mathclose{{}\left(\mathbf{q}^{(i)}_{[0]}+\epsilon\mathbf{q}^{(i)}_{[1]}}\right)\mathopen{}\mathclose{{}\left(\mathbf{q}^{(i)T}_{[0]}+\epsilon\mathbf{q}^{(i)T}_{[1]}}\right)\\\
=\sum_{i}\sigma^{(i)}_{[0]}\mathbf{q}^{(i)}_{[0]}\mathbf{q}^{(i)T}_{[0]}+\epsilon\sum_{i}\sigma^{(i)}_{[1]}\mathbf{q}^{(i)}_{[0]}\mathbf{q}^{(i)T}_{[0]}+\epsilon\sum_{i}\sigma^{(i)}_{[0]}\mathbf{q}^{(i)}_{[1]}\mathbf{q}^{(i)T}_{[0]}+\epsilon\sum_{i}\sigma^{(i)}_{[0]}\mathbf{q}^{(i)}_{[0]}\mathbf{q}^{(i)T}_{[1]}+O(\epsilon^{2}).$
(62)

We note that $\mathbf{q}^{(i)}_{[1]}\mathbf{q}^{(i)T}_{[0]}$ and
$\mathbf{q}^{(i)}_{[0]}\mathbf{q}^{(i)T}_{[1]}$ will be $O(\epsilon)$ since
Eq. 56 requires $\dot{\mathbf{q}}^{(i)}$ to be perpendicular to
$\mathbf{q}^{(i)}$ (with $O(\epsilon)$ error). Hence, this simplifies to

$\sum_{i}\mathopen{}\mathclose{{}\left(\sigma_{[0]}^{(i)}+\epsilon\sigma_{[1]}^{(i)}}\right)\mathopen{}\mathclose{{}\left(\mathbf{q}^{(i)}_{[0]}+\epsilon\mathbf{q}^{(i)}_{[1]}}\right)\mathopen{}\mathclose{{}\left(\mathbf{q}^{(i)T}_{[0]}+\epsilon\mathbf{q}^{(i)T}_{[1]}}\right)\\\
=\sum_{i}\sigma^{(i)}_{[0]}\mathbf{q}^{(i)}_{[0]}\mathbf{q}^{(i)T}_{[0]}+\epsilon\sum_{i}\sigma^{(i)}_{[1]}\mathbf{q}^{(i)}_{[0]}\mathbf{q}^{(i)T}_{[0]}+O(\epsilon^{2}).$
(63)

Hence, Eq. 61 becomes:

$\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\sum_{i}\sigma^{(i)}\mathbf{q}_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)T}+\epsilon
t_{[1]}\mathrm{tr}(\mathbf{L})\mathbf{T}_{[0]}-\epsilon
t_{[1]}\mathbf{T}_{[0]}\mathbf{L}^{T})}\right)+O(\epsilon^{2})=0.$ (64)

Then substituting for $\epsilon t_{[1]}\mathbf{L}$ using Eq. 58b we have:

$\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\sum_{i}\sigma^{(i)}\mathbf{q}_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)T}+\mathrm{tr}(\mathbf{F}^{*}-\mathbf{I})\mathbf{T}_{[0]}-\mathbf{T}_{[0]}\mathbf{F}^{*T}+\mathbf{T}_{[0]}}\right)+O(\epsilon^{2})=0.$
(65)

Using Eq. 53 gives:

$\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\sum_{i}\sigma^{(i)}\mathbf{q}_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)T}+\mathrm{tr}(\mathbf{F}^{*})\mathbf{T}_{[0]}-\mathbf{T}_{[0]}\mathbf{F}^{*T}}\right)+O(\epsilon^{2})=0.$
(66)

This gives the linearized form of force balance.

We denote $\mathbf{t}^{*}$ as the first Piola traction vector with respect to
the $\mathbf{x}_{[0]}$ reference frame such that:

$\mathrm{det}(\mathbf{F}^{*})\mathbf{T}\mathbf{F}^{*-T}\mathbf{n}_{[0]}=\mathbf{t}^{*}.$
(67)

We now consider the limit as $\epsilon\rightarrow 0$ to obtain the linearized
system:

$\displaystyle\mathrm{div}_{[0]}\mathopen{}\mathclose{{}\left(\sum_{i}\sigma^{(i)}\mathbf{q}_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)T}+\mathrm{tr}(\mathbf{F}^{*})\mathbf{T}_{[0]}-\mathbf{T}_{[0]}\mathbf{F}^{*T}}\right)=0,$
(68a) $\displaystyle\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}_{[0]}^{(j)}\mathbf{q}^{(j)T}_{[0]}\mathopen{}\mathclose{{}\left(\mathbf{L}\mathbf{B}^{e}_{[0]}+\mathbf{B}^{e}_{[0]}\mathbf{L}^{T}}\right)\mathbf{q}_{[0]}^{(i)}}{\lambda^{e(i)}_{\mathbf{B}[0]}-\lambda^{e(j)}_{\mathbf{B}[0]}},$
(68b)
$\displaystyle\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},\qquad\mathbf{B}^{e}=\sum_{i}\lambda^{e(i)}_{\mathbf{B}}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},$
(68c)
$\displaystyle\mathbf{B}^{e}=\bar{\mathbf{B}^{e}}(\mathbf{T})\qquad\textrm{in}\qquad\Omega_{[0]},$
(68d)
$\displaystyle\mathrm{det}(\mathbf{F}^{*})\mathbf{T}\mathbf{F}^{*-T}\mathbf{n}_{[0]}=\mathbf{t}^{*}\qquad\mathrm{on}\qquad\partial\Omega_{[0]}.$
(68e)

The linearized form corresponding to infinitesimal deformation is also
important because in general, deformations where linearization is valid will
be small. They can be obtained via the same procedure as above or by
simplifying the above equations. This gives the following linearized system
for infinitesimal deformation

$\displaystyle\mathrm{div}\mathopen{}\mathclose{{}\left(\sum_{i}\sigma^{(i)}\mathbf{q}_{[0]}^{(i)}\mathbf{q}_{[0]}^{(i)T}}\right)=0,$
(69a) $\displaystyle\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}_{[0]}^{(j)}\mathbf{q}^{(j)T}_{[0]}\dot{\mathbf{E}}\mathbf{q}_{[0]}^{(i)}}{\lambda^{e(i)}_{\mathbf{E}[0]}-\lambda^{e(j)}_{\mathbf{E}[0]}},$
(69b)
$\displaystyle\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},\qquad\mathbf{E}^{e}=\sum_{i}\lambda^{e(i)}_{\mathbf{E}}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},$
(69c)
$\displaystyle\mathbf{E}^{e}=\bar{\mathbf{E}^{e}}(\mathbf{T})\qquad\textrm{in}\qquad\Omega,$
(69d)
$\displaystyle\mathbf{T}\mathbf{n}=\mathbf{t}\qquad\mathrm{on}\qquad\partial\Omega.$
(69e)

##  4 Discussion

The system of governing equations derived has utility characterizing material
properties in the case of heterogeneous deformations. These may occur due to
instabilities and strain softening (e.g. see necking Section 5), heterogeneous
material properties (e.g. metal microstructures [41]), or during complex
deformations (e.g. crack growth). The necessary full field deformation
measurements required to apply the approach can often be obtained in two
dimensions using digital image correlation [1], the grid method [42], or
particle image velocimetry [43]. Furthermore, these full field measurement
approaches may be extended to three dimensions using tomographic techniques
[4, 5, 6]. We refer the reader to [23] for more detailed discussion of the
potential applications of these equations, as there is a significant overlap.
However, the equations presented here are more general and overcome a key
challenge: elastic unloading and non-negligible elastic strain. This is
important for a wide range of elasto-plastic deformations, even when the total
elastic strain is significantly less than the plastic strain. Consider that,
during a deformation, one region of the domain may be elastically unloading
while another is undergoing significant plastic deformation. It is not
possible to directly apply the equations developed in [23] here as one must
assume either elastic deformation or plastic deformation with negligible
elastic strain. This limitation of the prior method is illustrated for a
necking tensile specimen in [24]. However, the equations developed here can be
directly applied, as shown in the next section.

The modified formulation can be applied to elasto-visco-plastic deformation
when the elastic component of the constitutive equation is known. To apply the
approach, there are several requirements that must be satisfied:

  * •

The full-field deformation must be known
($\mathbf{x}=\bar{\mathbf{x}}(\mathbf{X},t)$ is sufficient as one can compute
$\mathbf{F}$ and other kinematic variables). Determining this with the
accuracy required may be challenging is discussed further below.

  * •

The traction boundary conditions must be known. The ability to satisfy this
requirement is highly problem dependent, for example, when deforming a dog
bone sample one may accurately assume that the stress is constant along the
gauge, far away from the edge or the neck. However, for may more complex
geometries this will likely be more challenging.

  * •

The initial stress state must be known. This can be met by simply considering
the deformation from the stress free state.

  * •

The body force must be known, if this is not negligible.

  * •

The initial density must be known, if the inertial terms are not negligible.
The density at later times can be computed using the mass conservation
equation if the material is compressible.

  * •

The elastic component of the constitutive equation must be known. This is
theoretically straightforward for cases where the relationship is isotropic-
linear-elastic, as the relationship could be determined at every point in the
material by applying the elastic method developed in [23] at the initial
stages of the deformation.

  * •

The constitutive relationship must meet the requirements specified in Section
2. I.e that the Kröner decomposition is applicable and that the elastic and
plastic components of the constitutive equation are isotropic. This is
discussed further below with examples of materials and deformation regimes.

Provided these conditions are met, the system of governing equations (Eq. 36)
can then be solved for any general geometry or boundary condition. We make
several additional remarks

  * •

Solutions can be numerically obtained utilizing the algorithm outlined in
Section 5.2. This is largely based on the finite volume method developed in
[24], which can be directly applied to general geometries and boundary
conditions in the case of infinitesimal and finite deformation.

  * •

The algorithms in the current form are not directly applicable to deformations
with repeated stress eigenvalues (though they are applicable to stress-free
initial conditions). The properties of the equations in this scenario are
under investigation. Other than this restriction, the equations and algorithms
in there current from are directly applicable to arbitrary deformations.

  * •

In the common case that only the two-dimensional deformation field is obtained
using an approach like DIC, the equations can still be solved provided the
material is in a state of plane stress, plane strain, or any other state where
the out of plane $\mathbf{T}$ shear components are zero (see Section 5).

  * •

The method can be applied to finite deformations, and can be straight-
forwardly applied to cases of large plastic deformation (Section 5), or fluid
flow.

  * •

The method is not restricted to cases where the elastic strain is much less
than the plastic strain, or the elastic strain rate is much less than the
plastic strain rate.

One important constraint on the practical application is that the deformation
field data must be sufficiently accurate to determine the principal directions
of the elastic strain. In cases where the elastic strain is typically small,
such as the deformation of metals, this may pose a significant challenge. To
quantify the effect of an error in the observed deformation $\mathbf{F}$,
which we denote $\mathbf{F}^{\prime}$, one must (i) understand the error
introduced in $\mathbf{q}$ and (ii) understand how the error in $\mathbf{q}$
translates to error in the stress. The latter is more complicated as it
involves non-local effects throughout the domain, and will likely be
investigated in future work. The constraints imposed by the former can be
understood in a simple way by considering the elastic and plastic limits of
the equation. First, we consider the elastic limit where
$\\{\mathbf{q}^{(i)}\\}$ are computed directly from
$\mathbf{B}^{e}+\mathbf{B}^{e\prime}$, where $\mathbf{B}^{e\prime}$ is the
error. Assuming $\mathbf{B}^{e\prime}<<\mathbf{B}^{e}$, we have that the error
in $\\{\mathbf{q}^{(i)}\\}$ will scale with $|\mathbf{B}_{0}^{e\prime}|$ due
to the theory of matrix perturbations (the subscript $0$ refers to the
deviatoric component of the $\mathbf{B}^{e}$). The error will be small if
$\mathbf{B}_{0}^{e\prime}<<\mathbf{B}_{0}^{e}$. Second, we consider a
deformation step in the plastic limit, where the principal directions are
computed directly from $\mathbf{D}^{p}+\mathbf{D}^{p\prime}$ where
$\mathbf{D}^{\prime}$ is the error. Again, assuming
$\mathbf{D}^{p\prime}<<\mathbf{D}^{p}$, we have that the error in
$\\{\mathbf{q}^{(i)}\\}$ will scale with $|\mathbf{D}_{0}^{p\prime}|$. The
error will be small if $\mathbf{D}_{0}^{\prime}<<\mathbf{D}_{0}$. In the case
with both elastic and plastic deformation, there will be non-linear
interaction of these errors. For example, during elastic unloading, the error
stress computed from the proceeding plastic deformation will translate to an
error in the stress. Nevertheless, the errors in the elastic and plastic
limits may prove useful guidelines until this non-linearity is more fully
understood. To summarize we require that:

$\begin{gathered}\mathbf{B}_{0}^{e\prime}<<\mathbf{B}_{0}^{e},\\\
\mathbf{D}_{0}^{p\prime}<<\mathbf{D}_{0}^{p}.\end{gathered}$ (70)

The accuracy of typical DIC strain data typically ranges between $1\times
10^{-3}$ and $1\times 10^{-5}$. Hence, the approach may be viable given
current levels of accuracy - although this clearly needs investigation using
experimental data and/or noisy synthetic detests. Here, depending on the noise
level, it may or may not be necessary to implement more specialized
discretizations, data smoothing methods, or other techniques. For discussion
and results related to the noise sensitivity of Eq. 1, the reader should refer
to [23].

We note that the approach will be valid for rate-dependent extensions of
isotropic elasto-plasticity, as all fundamental assumptions remain valid.
Furthermore, material dependence on temperature does not alter the underlying
assumptions. The incorporation of temperature and rate dependent deformation
makes the approach suited to elucidating high strain rate deformations via
full field deformation measurements - resolving issues that arise due to
heterogeneous deformation [44, 45]. While constitutive equations for numerous
materials are formulated utilizing the Kröner decomposition, the equations for
some elasto-visco-plastic materials are formulated using an evolution equation
instead. However, these alternate formulations are closely related and can
often be used to describe the same phenomena (e.g. see [46] and references
therein), so the approach presented may also be in many of these cases.

##  5 Case study: necking of a uniaxial tensile specimen

We demonstrate the utility of the governing equations by applying them to a
strain field output by a finite element simulation (solving the forward
problem). This ensures the assumptions are correct and allows us to validate
the stress computed, by comparing it to the stress output by the finite
element simulation. We consider a tensile test conducted on a two-dimensional
sample, inspired by typical tensile tests of ductile sheet metal where DIC is
readily available. In a real experiment, the sample will typically deform with
a relatively constant gauge cross-section, a diffuse neck will then develop,
then a localized neck will form and the sample will fracture [27]. As we
conduct only two-dimensional simulation, we exclude simulation of localized
necking. Assessing our algorithm on this problem will demonstrate the approach
on a problem where (i) the effects of finite deformation are significant, and
(ii) there are both elasto-plastically and elastically deforming regions.
Furthermore, the problem is of practical importance for obtaining the
constitutive equations for a sheet metal.

To validate the inverse problem approach, we follow the same broad methodology
explained in [23]. First, we run a forward problem computation using a
commercial finite element solver which outputs the stress and displacement
fields (Section 5.1). Second, we export the displacement field to use as input
to the inverse problem. We use this, along with the boundary conditions, to
calculate the stress-field using the inverse problem equations (Section 5.2).
Third, we compare the stress field output from the forward problem, with that
output from the inverse problem to assess the accuracy (Section 5.3).

###  5.1 The forward problem

We have a two-dimensional specimen subject to plane stress that is rectangular
in the undeformed frame with dimensions $l_{X}$ and $l_{Y}$ in the $X$ and $Y$
directions (Fig. 2). We have $l_{X}=6l_{Y}$. $(X,Y)$ corresponds the the
location in undeformed frame, and $(x,y)$ corresponds to the position in the
deformed frame. We have boundaries at $X=0$, $X=l_{X}$, $Y=0$ and $Y=l_{Y}$.

Figure 2: Domain and boundary conditions for forward and inverse problem. The
inverse problem domain is chosen as a subset of the forward problem domain
(highlighted in red). Green shows the boundary conditions: displacement for
the forward problem and traction for the inverse problem.

We specify the displacement on the $X=0$ and $X=l_{x}$ boundaries so that the
distance between the boundaries increases as a function of time. To ensure
that the sample necks somewhere near the center of the gauge, we specify that
the $y$ displacement is zero on these boundaries. We specify zero traction on
the $Y=0$ and $Y=l_{Y}$ boundaries.

We use the ABAQUS/expicit (2017) finite element solver for the forward problem
and specify the following boundary conditions:

$\begin{gathered}x(X=0,y)=0,\qquad y(X=0,y)=0,\\\ x(X=l_{x},y)=0,\qquad
y(X=l_{x},y)=\xi(t),\\\ \mathbf{t}(x,y=0)=0,\\\ \mathbf{t}(x,y=l_{y})=0.\\\
\end{gathered}$ (71)

Here, $\xi(t)$ is is a cubic function specified so that $d\xi(t=0)/dt=0$. We
specify a homogeneous elasto-plastic constitutive equation. The material has a
Young’s Modulus $E=100$ GPa and a Poisson’s ratio of $\nu=0.3$. For the
plastic properties, we specify isotropic hardening with a Von-Mesis yield
surface and the yield stress as a pairwise linear function of the total
plastic strain, $\varepsilon_{t}=\sqrt{2/3}|\mathbf{E}^{p}|$. This is shown
later in the article in Fig. 10, and was chosen to be similar to the copper
material investigated in [24]. More complex effects such as texturing of the
material [47] or damage nucleation [48] are not incorporated into the
simulation. The simulation is time dependent, with non-linear geometry so that
include finite strains and rotations are accounted for. We do not introduce
any geometrical imperfections to trigger necking, and instead rely upon
elastic waves propagating through the material to trigger the necking
instability. We use CPS4R elements which correspond to two-dimensional square
(in the undeformed frame), plane stress elements with nodes at each corner. We
have 40 elements in the $X$ direction and $240$ elements in the $Y$ direction.
We specify the simulation to run for 1 s with automatic time incrimination.
However, at approximately 0.6 s, which we refer to as $t_{f}$, the deformation
becomes highly localized in individual elements and the simulation is no
longer valid. Hence, we do not report results after this point. This
approximately corresponds to the point in time when the localized neck would
form. The strain output from this forward problem is shown in Fig. 3 for two
different points in time, this is used as input for the inverse problem. We
output the displacement of each node and the stress in each element at 2000
evenly spaced time intervals.

Figure 3: Strain output from forward problem at (a) $t=(11/12)t_{f}$ and (b)
$t=t_{f}$. Here we use the definition
$\varepsilon=\mathrm{ln}(1+\sqrt{2/3}|\mathbf{E}|)$, for this problem as this
is approximately equivalent to accumulated plastic strain.

###  5.2 Inverse problem

Here, we only consider a subset of the domain considered for the forward
problem. This allows us to avoid complications arising from the singularities
at the corners of the sample, for which the numerical method has not been
developed to address. The modified domain is shown in Fig. 2 where the length
in the $X$ direction is $l^{\prime}_{X}$ (1/8 of the domain is removed from
each side). We redefine $X$ such that $X=0$ on the left hand side boundary. We
have the traction on the $X=0$ and $X=l^{\prime}_{X}$ boundaries, obtained
directly from the forward problem. When a subset of the characteristics are
leaving the domain, we choose the component of the traction to specify such
that there is no reflection (see the discussion on boundary conditions in
[24]). We choose the direction of information propagation to be rightward for
the approximately horizontal characteristics, and upward for the approximately
vertical characteristics (Fig. 9). The reader should refer to [23, 24] for a
discussion of the characteristics. For the particular deformation we consider,
the inertial terms in force balance equation will be negligible, hence, we do
not include them in the inverse problem computation. We also use the elastic
constants specified in the forward problem directly (for a real experiment
they would have to be determined using one of the approaches discussed in
Section 2.2).

Figure 4: The domain and computational grid used for the inverse problem at
time $t=t_{f}$. Each cell in the figure is four finite volumes in the
computation.

We use the same approach to account for the variation in thickness as
discussed in [24]. The thickness $h$ is defined as a function of position and
time:

$h=\bar{h}(\mathbf{x},t).$ (72)

As the elastic volume change will be negligible for this problem, we assume
$\mathrm{det}(F_{ij})=1$. This gives $h/h_{0}=1/\mathrm{det}(F_{\alpha\beta})$
at every point in the sample throughout deformation where $h_{0}$ is the
initial thickness. Note we use Greek letters ($\alpha,\beta,...$) to indicate
the index when it takes values in the set $\\{1,2\\}$, and Latin letters
($i,j,...$) to indicate the index when it takes values in the set
$\\{1,2,3\\}$. Therefore, $\mathrm{det}(F_{\alpha\beta})$ corresponds to the
determinant of the two by two matrix in the $(x,y)$ plane. We use a modified
form of of the governing equations to account for varying thickness using the
variable $\mathbf{T}^{\prime}$ where $\mathbf{T}^{\prime}=h\mathbf{T}$:

$\displaystyle T_{ij}^{\prime}=hT_{ij}$ (73a) $\displaystyle\frac{\partial
T^{\prime}_{\alpha\beta}}{\partial x_{\beta}}=0,$ (73b)
$\displaystyle\dot{q}^{(\alpha)}_{\gamma}=\frac{q^{(\beta)}_{\gamma}q^{(\beta)}_{\delta}D_{\delta\zeta}q^{(\alpha)}_{\zeta}(\lambda^{e(\beta)}_{\mathbf{B}}+\lambda^{e(\alpha)}_{\mathbf{B}})}{\lambda^{e(\beta)}_{\mathbf{B}}-\lambda^{e(\alpha)}_{\mathbf{B}}}+W_{\gamma\delta}q^{(\alpha)}_{\delta}\quad\mathrm{where}\quad\alpha\neq\beta,$
(73c) $\displaystyle q^{(\alpha)}_{3}=0,$ (73d) $\displaystyle
T_{\alpha\beta}^{\prime}=\sum_{\gamma}\sigma^{\prime(\gamma)}q_{\alpha}^{(\gamma)}q_{\beta}^{(\gamma)},\qquad
B_{ij}^{e}=\sum_{k}\lambda^{e(k)}_{\mathbf{B}}q_{i}^{(k)}q_{j}^{(k)},$ (73e)
$\displaystyle B_{ij}^{e}=\bar{B^{e}}_{ij}(T_{ij}).$ (73f)

Here we implement an algorithm that is first order accurate in time. We use
$j$ to specify the increment in time and $n_{t}$ for the total number of
increments in time. We use the subscript to denote the variable evaluated at
that increment in time, e.g. $\mathbf{T}_{j}$ corresponds to the stress at the
$j$ increment. Instead of computing the rate of change of the $\mathbf{q}$
values directly, we compute $\mathbf{A}$ using Eq. 27. We obtain the principal
directions at the time step $j$ from $\mathbf{A}$:

$\mathbf{A}_{j}=\mathbf{F}^{*}\mathbf{B}^{e}_{j-1}\mathbf{F}^{*T}$ (74)

where $\mathbf{F}^{*}=\mathbf{F}_{j}\mathbf{F}_{j-1}^{-1}$. This is done for
each material point.

For this problem, $\mathbf{F}$ is computed for each element at each time step
using the displacements from the nodes at the corners of the element.
Specifically, linear regression is used to fit the displacement gradient.
$\mathbf{F}^{*}$, $\mathbf{B}$ and the thickness $h$ can then be computed from
$\mathbf{F}$. Once the principal directions $\\{\mathbf{q}_{j}^{(i)}\\}$ are
obtained from the computed $\mathbf{A}_{j}$, the finite volume method
developed in the [23] is used to obtain the stress as a function of position.
This stress is then used to compute the elastic strain, and then we can
proceed to time step $j+1$. The algorithm used is as follows:

$j\leftarrow 0$;

$\mathbf{T}_{j}\leftarrow 0$;

$j\leftarrow 1$;

compute $\\{\mathbf{q}^{(i)}_{j}\\}$ by solving
$\mathbf{D}_{j}\mathbf{q}^{(i)}_{j}=\lambda^{(i)}\mathbf{q}^{(i)}_{j}$ for all
$\mathbf{X}$;

while  _$j\leq n_{t}$_ do

compute $\mathbf{T}^{\prime}_{j}$ from $\\{\mathbf{q}_{j}\\}$ and boundary
conditions for all $\mathbf{X}$ using the finite volume method;

compute $\mathbf{T}_{j}$ from $\mathbf{T}^{\prime}_{j}$ using Eq. 73a and
$\mathbf{F}_{j}$ for all $\mathbf{X}$ using Eq. 72;

compute $\mathbf{B}^{e}_{j}$ from $\mathbf{T}_{j}$ for all $\mathbf{X}$ using
Eq. 37;

compute $\mathbf{A}_{j+1}$ from $\mathbf{B}^{e}_{j}$ and $\mathbf{F}$ for all
$\mathbf{X}$ using Eq. 74;

compute $\\{\mathbf{q}^{(i)}_{j+1}\\}$ by solving
$\mathbf{A}_{j+1}\mathbf{q}^{(i)}_{j+1}=\lambda^{(i)}\mathbf{q}^{(i)}_{j+1}$
for all $\mathbf{X}$;

$j\leftarrow j+1$;

end while

Algorithm 1 Compute $\mathbf{T}$ as a function of position and time

Deformation data is output from the forward problem for 2000 time steps,
however, as mentioned previously the deformation becomes highly localized in
individual elements and is not valid at approximately step 1220. Hence, we do
not report values at later steps. Furthermore, there is essentially zero
deformation in the first few time steps, as $d\xi(t)/dt=0$ at $t=0$. The
algorithm is not robust when there is close zero strain because computation
for the transmission and reflection coefficients result in singular matrices.
Hence, we begin the simulation at time step 20, and end at 1220. Hence, we
have have the number of steps $n_{t}=1200$, corresponding to each set of
output data form the forward problem. We also consider the problem with larger
time steps ($n_{t}=240,120,24$ and $12$) in order to evaluate how the error
depends on this value. The algorithm, including the finite volume method, was
implemented in Python and utilizes the Scipy sparse linear algebra package
[49]. The computation with 1200 time steps took approximately 4 hours on a on
a 2019 Macbook Pro with a 2.6 GHz Intel core i7 processor (the computational
time was approximately proportional to $n_{t}$).

###  5.3 Results

We compare the stress fields computed from the inverse problem with those from
the forward problem in Figs. 5 and 6. As can be seen, the stress distributions
are in close agreement, indicating close to exact results and demonstrating
the validity of the approach.

Figure 5: Comparison of the stress field output by the forward problem and
inverse problem at time $t=(11/12)t_{f}$ and $t=t_{f}$. Figure 6: Line
comparison of the stress output by the forward problem and inverse problem at
time $t=(11/12)t_{f}$ (a1,b1,c1,d1) and $t=t_{f}$ (a2,b2,c2,d2). (a) $T_{xx}$
stress component along $Y=l_{Y}/2$, (b) $T_{xy}$ stress component on a
$X=(7/18)l_{X}$ line that passes through the neck, but not at the center (so
the shear is non-zero), (c) $T_{yy}$ component along component along
$Y=l_{Y}/2$, (d) $T_{xy}$ stress component on a $Y=l_{Y}/4$ line.

We quantify the error as a function of position using
$|\Delta\mathbf{T}|/\mathrm{mean}_{x}({\mathbf{T}})$, where $\Delta\mathbf{T}$
is the difference between the stress computed from the forward problem and the
stress computed from the inverse problem, and the $\mathrm{mean}_{x}$ is taken
over all positions (but only at the time considered). The $\mathbf{T}$ used in
the denominator is the $\mathbf{T}$ from the forward problem. We show the
error between the distributions in Fig. 7. As can be seen in Fig. 7a, the
error increases toward the right hand side. This is because the boundary
condition is specified on the left hand side and as the solution is computed
along characteristics going from left to right, the numerical error
accumulates. One can also see that the error increases near the boundary to
the right of the neck. This is due to the characteristic lines not being
perfectly parallel to the boundary, and is particularly apparent for large
time-steps (the characteristic lines are discussed further later in this
section). Finally, we note that we see the larger errors in the center of the
neck. This is likely because $\\{\mathbf{q}^{(i)}\\}$ is changing more rapidly
with time, and the finite time step will cause a larger error.

Figure 7: Error as a function of position when compared to the forward
problem, plotted at different times for solutions computed using different
temporal resolutions.

We also quantify the error as a function of time using the normalized mean
absolute error (NMAE):

$\mathrm{NMAE}(\mathbf{T})=\frac{\mathrm{mean}_{x}(|\Delta\mathbf{T}|)}{\mathrm{mean}_{x}(|\mathbf{T}|)}.$
(75)

Note that elements at the boundary are not included in the error computation.
There is moderate error accumulation at the boundary elements when using the
current numerical approach, but this can practically can be removed during
post processing when using the data to investigate real materials or
deformations. Hence, the authors considered it more representative not to
include this in the overall error metric (if included the boundary error would
largely control the metric). The reader can observe the magnitudes of these
errors in Fig. 6b2,c1, and can also refer to [24] for a more extensive
discussion on this issue.

The NMAE error value obtained at the end of the simulation was 2.1$\times
10^{-3}$ at the time when the simulation was terminated (this was the point in
time with the maximum error). The error is plotted as a function of time for
different sizes of time-step (Fig. 8). As can be seen in Fig. 8a, the error
decreases as the number of temporal increments are increased. There are
several possible sources for error: the temporal resolution, the spatial
resolution, approximation of constant volume when computing the thickness,
inertial terms in the force balance equation, and error in the forward problem
computation. As the error is decreasing with increasing $n_{t}$, this
indicates that the temporal resolution is the primary source of error. This is
particularly true at higher values of $t$ where necking occurs. It is likely
that with further reduction of $n_{t}$, the error would decrease further.

Figure 8: NMAE error. (a) as a function of the resolution at two different
times. It can be seen that increasing the number of steps $n_{t}$ is likely to
improve accuracy further. (b) Error as a function of time for five different
values of $n_{t}$. It can be seen that error increases substantially near the
end where necking and significant rotations start to occur.

The characteristic lines, lines which are aligned with the principal
directions of the stress at each point in space, are critical for the
numerical computation (see [23, 24]). As the governing equations are
hyperbolic, propagation of information along these lines is used to compute
the solution at each time step. For visualization purposes, the characteristic
lines corresponding to the stress computation at specific times are computed
using an ODE solver (Fig. 9). As the computation is two-dimensional, there are
two sets of characteristic lines corresponding to each principal direction.
The horizontal characteristic lines are shown, and the arrows indicate the
arbitrarily chosen direction of information propagation used in the
computational procedure. As can be seen from Fig. 9c, when the time step is
large there is substantial error in the characteristic lines. I.e. they should
be aligned with the zero traction boundary at the boundary, as there will be
no shear traction applied. However, this error diminishes as the time step is
made smaller.

Figure 9: Spatial characteristic lines obtained when solving the inverse
problem. These are obtained using an ODE solver and the principal directions.
There is an additional set of characteristic lines at right angles with
information propagating approximately in the upward direction.

We envision that this approach could be used to obtain the stress-strain
relationship for a material at high strains (i.e. by plotting the stress and
strain for an element in the neck). We show this plot for the material element
subject to the highest stress in Fig. 10. This is compared with the
relationship specified in the forward problem, and one can see they are a
close match. When plotting the stress-strain relationship for a material
element near the $X=0$ boundary, we see the elastic unloading that occurs at
the same time the neck forms in the center of the sample. We contrast these
curves with the stress strain relationship that would be obtained using the
conventional method used to analyze tensile specimens: measuring $\Delta l$
along the gauge, calculating the engineering strain, computing the engineering
stress from the force, and calculating the true stress and strain assuming a
constant cross-section. The methods are both accurate for strain values up to
necking ($\varepsilon=0.26$). However, information at higher strains is not
available. The approach presented in this paper gives accurate stress-strain
behavior up to $\varepsilon=0.90$, this corresponds to an increase in length
of 142%. Other approaches may be applied to estimate the stress in the neck,
as discussed in the introduction, however, more restrictive assumptions must
be made regarding the constitutive equations, and significant errors still
arise as the solutions are inexact. For example, the Bridgman correction
(applicable to cylindrical geometries), gives errors on the order of 10% for
strains on the order of $70-80\%$. In contrast, the approach developed here is
theoretically exact and does not need to be adapted for specific material
hardening laws - so could be a promising approach. Of course, before
application to real materials, the sensitivity to noise must be investigated,
along with the validity of additional assumptions such as constant elastic
modulus, and the lack of texturing of the sample. These results simply provide
a proof of principal demonstration.

Figure 10: Stress-strain curve for material elements obtained via the inverse
problem, compared to the specified relationship in the forward problem. We
plot the Von Mesis stress $\sigma=\sqrt{3/2}|\mathbf{T}_{0}|$ and strain
measure $\varepsilon=\mathrm{ln}(1+\sqrt{2/3}|\mathbf{E}|)$. In addition, we
plot the stress-strain curve that would be obtained using a conventional
analysis that assumes a uniform cross section (the dotted line corresponds to
when the assumption breaks down). Also note the lines directly overlay each
other until necking and some are not visible.

##  6 Conclusion

The method developed in [23], previously valid for isotropic elastic or
plastic deformation, is generalized to elasto-plastic deformation. An
additional assumption, that the elastic component of the constitutive equation
is known, is used to achieve this. This assumption appears to be reasonable
for materials with simple linear elastic elasticity and more complex plastic
deformation, such as metals. This generalization significantly broadens the
range of applicability of this approach as there will frequently be both
elastic and plastic regions during a complex deformation. Unlike the prior
formulation, the governing partial differential equations are non-linear and
time dependent. In future work it will be important to validate the governing
equations by applying them to a broader range of problems such as those in
three-dimensions, those with non-negligible accelerations, and those with more
complex geometries. Though we do not anticipate any issues, it is possible
that unexpected challenges may arise in these scenarios. Furthermore, a
detailed investigation into noise will be required before application to
experimental data. The approach approach could be applied to other localized
deformations such as shear band formation to enable measurement of the
constitutive equations in these deformation regimes [50]. Furthermore, the
approach may be able to be applied to heterogeneous material microstructures
to be used as input into data driven microstructural optimizations [51] or
material knowledge systems [52]. This paper presents one generalization of the
governing equations, however, other generalizations to phenomena such as
anisotropic elasticity or crystal plasticity may also be possible.

The example used to demonstrate the method, is itself important (as discussed
in the introduction). There are a wide range of methods to approximately solve
the problem, but this is the first approach (to the best of our knowledge)
which gives a a theoretically exact solution consistent with force balance and
an observed full-field deformation.

## References

  * [1] Bing Pan, Kemao Qian, Huimin Xie, and Anand Asundi.  Two-dimensional digital image correlation for in-plane displacement and strain measurement: A review.  Measurement Science and Technology, 20(6), 2009.
  * [2] J. Kang, D. S. Wilkinson, J. D. Embury, and M. Jain.  Microscopic Strain Mapping Using Scanning Electron Microscopy Topography Image Correlation at Large Strain.  The Journal of Strain Analysis for Engineering Design, 40(6):559–570, 2005.
  * [3] Dingshun Yan, C. Cem Tasan, and Dierk Raabe.  High resolution in situ mapping of microstrain and microstructure evolution reveals damage resistance criteria in dual phase steels.  Acta Materialia, 96(August):399–409, 2015.
  * [4] N. Lenoir, M. Bornert, J. Desrues, P. Bésuelle, and G. Viggiani.  Volumetric digital image correlation applied to x-ray microtomography images from triaxial compression tests on argillaceous rock.  Strain, 43(3):193–205, 2007.
  * [5] Stéphane Roux, François Hild, Philippe Viot, and Dominique Bernard.  Three-dimensional image correlation from X-ray computed tomography of solid foam.  Composites Part A: Applied Science and Manufacturing, 39(8):1253–1265, 2008.
  * [6] E. Tudisco, S. A. Hall, E. M. Charalampidou, N. Kardjilov, A. Hilger, and H. Sone.  Full-field Measurements of Strain Localisation in Sandstone by Neutron Tomography and 3D-Volumetric Digital Image Correlation.  Physics Procedia, 69(October 2014):509–515, 2015.
  * [7] H. D. Harris and D. M. Bakker.  A soil stress transducer for measuring in situ soil stresses.  Soil and Tillage Research, 29(1):35–48, 1994.
  * [8] K Ramesh.  Digital photoelasticity.  Springer, Berlin, 2000.
  * [9] T Vermeij, E Plancher, and C C Tasan.  Preventing damage and redeposition during focused ion beam milling: the “umbrella” method.  Ultramicroscopy, 186:35–41, 2018.
  * [10] Michel Grédiac, Evelyne Toussaint, and Fabrice Pierron.  Special virtual fields for the direct determination of material parameters with the virtual fields method. 1-Principle and definition.  international journal of solids and structures, 39:15, 2002.
  * [11] Eric Florentin and Gilles Lubineau.  Using constitutive equation gap method for identification of elastic material parameters: technical insights and illustrations.  International Journal on Interactive Design and Manufacturing (IJIDeM), 5(4):227–234, 2011.
  * [12] Fabrice Pierron and Michel Grédiac.  The virtual fields method: extracting constitutive mechanical parameters from full-field deformation measurements.  Springer Science & Business Media, 2012.
  * [13] Aleksander Marek, Frances M Davis, and Fabrice Pierron.  Sensitivity-based virtual fields for the non-linear virtual fields method.  Computational mechanics, 60(3):409–431, 2017.
  * [14] Laurent Crouzeix, Jean Noël Périé, Francis Collombet, and Bernard Douchin.  An orthotropic variant of the equilibrium gap method applied to the analysis of a biaxial test on a composite material.  Composites Part A: Applied Science and Manufacturing, 40(11):1732–1740, 2009.
  * [15] Eric Florentin and Gilles Lubineau.  Identification of the parameters of an elastic material model using the constitutive equation gap method.  Computational Mechanics, 46(4):521–531, 2010.
  * [16] Julien Réthoré.  A fully integrated noise robust strategy for the identification of constitutive laws from digital images.  International Journal for Numerical Methods in Engineering, 84(6):1885–1891, 2010.
  * [17] M. Z. Siddiqui, S. Z. Khan, M. A. Khan, M. Shahzad, K. A. Khan, S. Nisar, and D. Noman.  A Projected Finite Element Update Method for Inverse Identification of Material Constitutive Parameters in transversely Isotropic Laminates.  Experimental Mechanics, 57(5):755–772, 2017.
  * [18] Romain Viala, Vincent Placet, and Scott Cogan.  Identification of the anisotropic elastic and damping properties of complex shape composite parts using an inverse method based on finite element model updating and 3D velocity fields measurements (FEMU-3DVF): Application to bio-based composite violin sou.  Composites Part A: Applied Science and Manufacturing, 106:91–103, 2018.
  * [19] J. C. Gelin and O. Ghouati.  An inverse method for determining viscoplastic properties of aluminium alloys.  Journal of Materials Processing Tech., 45(1-4):435–440, 1994.
  * [20] Adrien Leygue, Michel Coret, Julien Réthoré, Laurent Stainier, and Erwan Verron.  Data-based derivation of material response.  Computer Methods in Applied Mechanics and Engineering, 331:184–196, 2018.
  * [21] Harshad M. Paranjape, Kenneth I. Aycock, Craig Bonsignore, Jason D. Weaver, Brent A. Craven, and Thomas W. Duerig.  A probabilistic approach with built-in uncertainty quantification for the calibration of a superelastic constitutive model from full-field strain data.  Computational Materials Science, 192(July 2020):110357, 2021.
  * [22] Sevan Goenezen, Baik Jin Kim, Maulik Kotecha, Ping Luo, and Mohammad R. Hematiyan.  Mechanics Based Tomography (MBT): Validation using experimental data.  Journal of the Mechanics and Physics of Solids, 146, 2021.
  * [23] Benjamin C. Cameron and Cemal Cem Tasan.  Full-field stress computation from measured deformation fields: a hyperbolic formulation.  Journal of mechanics and physics of solids, 147, 2021.
  * [24] Benjamin C.; Cameron and C. C. Tasan.  On discontinuities when computing the stress-field from the strain: a finite volume discretization (in submission).
  * [25] P W Bridgman.  Studies in Large Plastic Flow and Fracture.  McGraw-Hill, New York - London, 1952.
  * [26] Z L Zhang, M Hauge, J Odegård, and C Thaulow.  Determining material true stressÐstrain curve from tensile specimens with rectangular cross-section.  International Journal of Solids and Structures, 36(9919), 1998.
  * [27] Amit K Ghosh.  Strain localization in the diffuse neck in sheet metal.  Metallurgical Transactions, 5(7):1607–1616, 1974.
  * [28] B. C. Cameron, M. Koyama, and C. C. Tasan.  Phase stability effects on hydrogen embrittlement resistance in martensite–reverted austenite steels.  Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science, 50(1):29–34, 2019.
  * [29] Jörgen Kajberg and Göran Lindkvist.  Characterisation of materials subjected to large strains by inverse modelling based on in-plane displacement fields.  International Journal of Solids and Structures, 41(13):3439–3459, 2004.
  * [30] J-H Kim, A Serpantié, F Barlat, F Pierron, and M-G Lee.  Characterization of the post-necking strain hardening behavior using the virtual fields method.  International Journal of Solids and Structures, 50(24):3829–3842, 2013.
  * [31] Sam Coppieters, Sousi Sumita, Daisaku Yanaga, Kristof Denys, Dimitri Debruyne, and Toshihiko Kuwabara.  Identification of post-necking strain hardening behavior of pure titanium sheet.  In Residual Stress, Thermomechanics & Infrared Imaging, Hybrid Techniques and Inverse Problems, Volume 9, pages 59–64. Springer, 2016\.
  * [32] K S Zhano and Z H Li.  Numerical analysis of the stress-strain curve and fracture initiation for ductile material.  Engineering fracture mechanics, 49(2):235–241, 1994.
  * [33] Shengwen Tu, Xiaobo Ren, Jianying He, and Zhiliang Zhang.  Stress–strain curves of metallic materials and post-necking strain hardening characterization: A review.  Fatigue and Fracture of Engineering Materials and Structures, 43(1):3–19, 2020.
  * [34] Hung-yang Yeh and Jung-ho Cheng.  NDE of metal damage : ultrasonics with a damage mechanics model.  International Journal of Solids and Structures, 40:7285–7298, 2003\.
  * [35] R M Cleveland and A K Ghosh.  Inelastic effects on springback in metals.  International Journal of Plasticity, 18:769–785, 2002.
  * [36] Ch. Feichter, Z Major, and R W Lang.  Methods for Measuring Biaxial Deformation on Rubber and Polypropylene Specimens.  In E E Gdoutos, editor, Experimental Analysis of Nano and Engineering Materials and Structures, pages 273–274, Dordrecht, 2007. Springer Netherlands.
  * [37] Auricchio Ferdinando and Robert L Taylor.  Shape-memory alloys: modelling and numerical simulations of the finite-strain superelastic behavior.  Computer Methods in Applied Mechanics and Engineering, 143:175–194, 1997.
  * [38] I. U. Ojalvo.  Gradients for large structural models with repeated frequencies.  SAE Technical Papers, 1986.
  * [39] M. I. Friswell.  The derivatives of repeated eigenvalues and their associated eigenvectors.  Journal of Vibration and Acoustics, Transactions of the ASME, 118(3):390–397, 1996.
  * [40] Klaus Jürgen Bathe and Haluk Ozdemir.  Elastic-plastic large deformation static and dynamic analysis.  Computers and Structures, 6(2):81–92, 1976.
  * [41] Y J Liu, Z Liu, Y Jiang, G W Wang, Y Yang, and L C Zhang.  Gradient in microstructure and mechanical property of selective laser melted AlSi10Mg.  Journal of Alloys and Compounds, 735:1414–1421, 2018.
  * [42] H. Chalal, S. Avril, F. Pierron, and F. Meraghni.  Experimental identification of a nonlinear model for composites using the grid technique coupled to the virtual fields method.  Composites Part A: Applied Science and Manufacturing, 37(2):315–325, 2006.
  * [43] R. J. Adrian.  Twenty years of particle image velocimetry.  Experiments in Fluids, 39(2):159–169, 2005.
  * [44] H. Kolsky.  An Investigation of the Mechanical Properties of Materials at very High Rates of Loading.  Proceedings of the Physical Society. Section B, pages 676–700, 1949\.
  * [45] Jonathan B. Estrada, Carlos Barajas, David L. Henann, Eric Johnsen, and Christian Franck.  High strain-rate soft material characterization via inertial cavitation.  Journal of the Mechanics and Physics of Solids, 112:291–317, 2018\.
  * [46] M. B. Rubin and A. L. Yarin.  On the relationship between phenomenological models for elastic-viscoplastic metals and polymeric liquids.  Journal of Non-Newtonian Fluid Mechanics, 50(1):79–88, 1993.
  * [47] Surya R. Kalidindi, C. A. Bronkorst, and Lallit Anand.  Crystallographic texture evolution in bulk deformation processing of FCC metals.  Journal of mechanics and physics of solids, 40(3):537–569, 1990\.
  * [48] J. Kováčik.  Correlation between Young’s modulus and porosity in porous materials.  Journal of Materials Science Letters, 18(13):1007–1010, 1999.
  * [49] Eric Jones, Travis Oliphant, Pearu Peterson, and others.  Open source scientific tools for Python, 2001.
  * [50] J Réthoré, F. Hild, and S. Roux.  Shear-band capturing using a multiscale extended digital image correlation technique.  Computational methods in applied mechanics and engineering, 196:5016–5030, 2007.
  * [51] Benjamin C. Cameron and C. Cem Tasan.  Microstructural damage sensitivity prediction using spatial statistics.  Scientific reports, 9:2774, 2019.
  * [52] Surya R. Kalidindi.  Microstructure Informatics.  Informatics for Materials Science and Engineering: Data-Driven Discovery for Accelerated Experimentation and Application, pages 443–466, 2013\.

##  Appendix A Appendix: Infinitesimal deformation governing equations

Here we give a derivation of the governing elasto-plastic deformation
equations in the simplified case of infinitesimal strain. This argument
follows parallel to that presented in Section 2.3. We also note that the same
result can be obtained by simplifying the more general finite deformation
equations (see Section 2.4).

The deformation is considered from time $t=\tau$ to time $t=\tau+\Delta t$. It
is assumed that the Cauchy stress $\mathbf{T}$ is known at time $\tau$ and we
wish to determine the stress at time $\tau+\Delta t$. We assume the standard
additive decomposition for $\mathbf{E}$

$\mathbf{E}=\mathbf{E}^{e}+\mathbf{E}^{p},$ (76)

where $\mathbf{E}^{e}$ is the elastic component and $\mathbf{E}^{p}$ is the
plastic component. Subscripts are introduced to denote quantities evaluated at
these times, e.g. $\mathbf{T}_{\tau}=\bar{\mathbf{T}}(\tau)$ and
$\mathbf{T}_{\tau+\Delta t}=\bar{\mathbf{T}}(\tau+\Delta t)$. We define
$\mathbf{E}^{*}$ to correspond to the deformation from $\tau$ to ${\tau+\Delta
t}$, i.e we have

$\mathbf{E}_{\tau+\Delta t}=\mathbf{E}_{\tau}+\mathbf{E}^{*}.$ (77)

We define $\mathbf{E}^{p*}$ in a similar way

$\mathbf{E}^{p}_{\tau+\Delta t}=\mathbf{E}^{p}_{\tau}+\mathbf{E}^{p*}.$ (78)

Substituting these expressions into Eq. 76 gives

$\mathbf{E}^{*}+\mathbf{E}^{e}_{\tau}+\mathbf{E}^{p}_{\tau}=\mathbf{E}^{e}_{\tau+\Delta
t}+\mathbf{E}^{p*}+\mathbf{E}^{p}_{\tau},$ (79)
$\mathbf{E}^{*}+\mathbf{E}^{e}_{\tau}=\mathbf{E}^{e}_{\tau+\Delta
t}+\mathbf{E}^{p*}.$ (80)

We express the spectral decomposition of $\mathbf{T}$ as

$\mathbf{T}\mathbf{q}^{(i)}=\sigma^{(i)}\mathbf{q}^{(i)}\qquad\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},$
(81)

where $\\{\mathbf{q}^{(i)}|i\in 1,2,3\\}$ are the eigenvectors of the stress
and $\\{\sigma^{(i)}|i\in 1,2,3\\}$ are the eigenvalues. We also assume the
eigenvalues do not repeat. Note that $\mathbf{E}^{e}$, $\dot{\mathbf{E}}^{p}$
and $\mathbf{T}$ have the same $\\{\mathbf{q}^{(i)}\\}$ (see Section 2.1).
Hence, we have

$\mathbf{E}_{\tau+\Delta t}^{e}\mathbf{q}_{\tau+\Delta
t}^{(i)}=\lambda^{(i)}\mathbf{q}_{\tau+\Delta t}^{(i)},$ (82)
$\dot{\mathbf{E}}_{\tau+\Delta t}^{p}\mathbf{q}^{(i)}_{\tau+\Delta
t}=\lambda^{(i)}\mathbf{q}_{\tau+\Delta t}^{(i)}.$ (83)

where $\\{\lambda^{(i)}|i\in 1,2,3\\}$ are generic eigenvalues that will be
different in Eqs. 82 and 83. In contrast, $\\{\mathbf{q}^{(i)}\\}$ always
refers to the eigenvectors of the stress and will be the same in all
equations. We consider the expression

$\mathbf{E}^{e}_{\tau+\Delta t}+\dot{\mathbf{E}}^{p}_{\tau+\Delta t}\Delta t.$
(84)

This will have the same principal directions as $\mathbf{T}_{\tau+\Delta t}$
as the sum of two tensors with the same eigenvectors also will have the same
eigenvector. We are also free to multiply any tensor by a scalar quantity.
Hence we have

$(\mathbf{E}^{e}_{\tau+\Delta t}+\dot{\mathbf{E}}^{p}_{\tau+\Delta t}\Delta
t)\mathbf{q}_{\tau+\Delta t}^{(i)}=\lambda^{(i)}\mathbf{q}_{\tau+\Delta
t}^{(i)}.$ (85)

In addition, we note that

$\mathbf{E}^{*}=\dot{\mathbf{E}}_{\tau+\Delta t}\Delta t+O({\Delta t}^{2}),$
(86) $\mathbf{E}^{p*}=\dot{\mathbf{E}}^{p}_{\tau+\Delta t}\Delta t+O({\Delta
t}^{2}).$ (87)

Hence, we have

$(\mathbf{E}^{e}_{\tau+\Delta t}+\mathbf{E}^{p*})\mathbf{q}_{\tau+\Delta
t}^{(i)}+O({\Delta t}^{2})=\lambda^{(i)}\mathbf{q}_{\tau+\Delta t}^{(i)}.$
(88)

Combining this with Eq. 80 then gives

$(\mathbf{E}^{*}+\mathbf{E}^{e}_{\tau})\mathbf{q}_{\tau+\Delta
t}^{(i)}+O({\Delta t}^{2})=\lambda^{(i)}\mathbf{q}_{\tau+\Delta t}^{(i)}.$
(89)

In the case where ${\Delta t}$ is small this gives an expression for the
principal directions of $\mathbf{T}_{\tau+\Delta t}$ in terms of the observed
quantity $\mathbf{E}^{*}$ and $\mathbf{E}^{e}_{\tau}$. The latter of which can
be computed from $\mathbf{T}_{\tau}$ using the assumption that the elastic
component of the constitutive equation is known (Section 2.2). We can use this
expression to derive a differential equation for the evolution of
$\\{\mathbf{q}^{(i)}\\}$.

For some arbitrary tensor $\mathbf{A}$ with eigenvectors
$\\{\mathbf{q}_{\mathbf{A}}^{(i)}\\}$ and eigenvalues
$\\{\lambda_{\mathbf{A}}^{(i)}\\}$ that varies with time,
$\mathbf{A}=\bar{\mathbf{A}}(t)$, we have the following general expression for
the derivative of its eigenvectors

$\dot{\mathbf{q}}^{(i)}_{\mathbf{A}}=\sum_{j,j\neq
i}\frac{\mathbf{q}_{\mathbf{A}}^{(j)}\mathbf{q}_{\mathbf{A}}^{(j)T}\dot{\mathbf{A}}\mathbf{q}_{\mathbf{A}}^{(i)}}{\lambda^{(i)}_{\mathbf{A}}-\lambda^{(j)}_{\mathbf{A}}},$
(90)

where this assumes there are no repeated eigenvalues and ensures that
$\\{\mathbf{q}_{\mathbf{A}}^{(i)}\\}$ remain normalized.

Let $\mathbf{A}=\mathbf{E}^{*}+\mathbf{E}^{e}_{\tau}$. This will have the same
eigenvectors as the stress up to error of $O(\Delta t^{2})$ (Eq. 89). Hence,
in the limit of $\Delta t\rightarrow 0$, we have
$\mathbf{q}^{(i)}_{\mathbf{A}}=\mathbf{q}^{(i)}$ and
$\dot{\mathbf{q}}^{(i)}_{\mathbf{A}}=\dot{\mathbf{q}}^{(i)}$. We will
determine an expression for $\dot{\mathbf{A}}$ and
$\\{\lambda_{\mathbf{A}}^{(i)}\\}$ to derive an expression for
$\dot{\mathbf{q}}^{(i)}$. Technically, $\mathbf{A}$ will be a function $\tau$
and $\tau+{\Delta t}$ as $\mathbf{E}^{e}_{\tau}=\bar{\mathbf{E}^{e}}(\tau)$
and $\mathbf{E}^{*}$ depends on both $\tau$ and $\tau+{\Delta t}$, i.e.
$\bar{\mathbf{E}^{*}}(\tau,\tau+{\Delta t})$. However, because we wish to
determine $\\{\mathbf{q}^{(i)}_{\tau+\Delta t}\\}$ using Eq. 89, we consider
$\tau$ fixed and allow time $t=\tau+\Delta t$ to vary. Hence,
$\bar{\mathbf{A}}(t)=\bar{\mathbf{E}^{*}}(\tau,t)+\bar{\mathbf{E}^{e}}(\tau)$.
Evaluating $\dot{\mathbf{A}}$ gives

$\displaystyle\dot{\mathbf{A}}$ $\displaystyle=\lim_{\Delta t\rightarrow
0}\frac{\bar{\mathbf{A}}(t+\Delta t)-\bar{\mathbf{A}}(t)}{\Delta t},$ (91)
$\displaystyle=\lim_{{\Delta t}\rightarrow
0}\frac{\bar{\mathbf{E}^{*}}(\tau,\tau+\Delta t)}{{\Delta t}},$

where we used the fact that $\bar{\mathbf{E}^{*}}(\tau,\tau)=0$. Substituting
in the expression for Eq. 86 gives

$\dot{\mathbf{A}}=\lim_{{\Delta t}\rightarrow 0}\frac{\dot{\mathbf{E}}{\Delta
t}+O({\Delta t}^{2})}{{\Delta t}},$ (92) $\dot{\mathbf{A}}=\dot{\mathbf{E}}.$
(93)

Evaluating $\\{\lambda_{\mathbf{A}}^{(i)}\\}$ at time $\tau=t$ gives
$\lambda_{\mathbf{A}}^{(i)}=\lambda^{e(i)}_{\mathbf{E}}$ where
$\\{\lambda_{\mathbf{E}}^{e(i)}\\}$ are the eigenvalues of $\mathbf{E}^{e}$.

Substituting these quantities into Eq. 90 gives

$\dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\dot{\mathbf{E}}\mathbf{q}^{(i)}}{\lambda^{e(i)}_{\mathbf{E}}-\lambda^{e(j)}_{\mathbf{E}}}.$
(94)

This, combined with force balance and the assumption that the elastic
component of the constitutive equation is known gives a full deterministic
system of partial differential equations

$\displaystyle\mathrm{div}{\mathbf{T}}+\mathbf{b}=\rho\ddot{\mathbf{x}},$
(95a) $\displaystyle\ \dot{\mathbf{q}}^{(i)}=\sum_{j,j\neq
i}\frac{\mathbf{q}^{(j)}\mathbf{q}^{(j)T}\dot{\mathbf{E}}\mathbf{q}^{(i)}}{\lambda^{e(i)}_{\mathbf{E}}-\lambda^{e(j)}_{\mathbf{E}}},$
(95b)
$\displaystyle\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},\qquad\mathbf{E}^{e}=\sum_{i}\lambda^{e(i)}_{\mathbf{E}}\mathbf{q}^{(i)}\mathbf{q}^{(i)T},$
(95c) $\displaystyle\mathbf{E}^{e}=\bar{\mathbf{E}^{e}}(\mathbf{T}).$ (95d)

Similarly to the finite deformation case, these equations are non-linear and
time dependent.