# On discontinuities when computing the stress-field from the strain: A finite
volume discretization

Benjamin C. Cameron  Cem C. Tasan<EMAIL_ADDRESS>
###### Abstract

This article focuses on developing a finite volume discretization for a
recently-developed and widely-applicable system of PDEs, enabling computation
of the full-field stress from a measured strain field and traction boundary
conditions. This is applicable to datasets obtained via DIC, and applies when
the material properties are unknown. Interestingly, the strain or strain-rate
discontinuities at the interfaces between finite volumes cause reflection and
transmission of the principal stresses. Hence, we analytically analyze the
solution structure in this case, enabling the development of the
discretization. In addition, analysis is conducted translating knowledge of
the traction at the boundary into boundary conditions, so that the numerical
method can be applied to a variety of loading conditions. The finite volume
discretization is validated using strain fields obtained from finite element
simulations. The developed discretization enables the system of equations to
be applied to irregular geometries, finite deformation, general boundary
conditions, and avoids instabilities that arise due to the Courant-Friedrichs-
Lewy condition, significantly increasing the range of applications for which
these PDEs can be applied.

###### keywords:

constitutive behavior , characteristics , mechanical testing , numerical
algorithms , inverse problem

††journal: International Journal of Solids and Structures

[label1]organization=Department of Civil and Environmental Engineering
Massachusetts Institute of Technology, [label2]organization=Department of
Materials Science and Engineering
Massachusetts Institute of Technology,

Developed computational method enables determination of constitutive equations
from digital image correlation (DIC).

Theoretical solution obtained for discontinuity giving principal stress
reflection.

Applies to isotropic elastic, elastoplastic, plastic and viscous deformation.

Alternative procedure to the virtual fields method and finite element model
updating.

Demonstrated computation of stress-strain relationship for strains on the
order of 100%.

##  1 Introduction

The availability and quality of full-field deformation data have substantially
improved in recent years. Image analysis algorithms, such as digital image
correlation (DIC), have become more accurate [1, 2], and there has been a
rapid increase in software and hardware availability [3]. These approaches can
be used to quantify complex plastic deformations, such as shear band
formation, necking, and high strain-rate deformation. Although extensive full-
field strain data is often available, the full-field stress typically remains
unknown. Obtaining the full-field stress would enable (i) quantification of
material properties when heterogeneous plastic flow instabilities arise; (ii)
computation of material properties as a function of position; and (iii)
quantification of stress partitioning in complex microstructures and composite
materials. This motivates solutions to the inverse problem: computing the
full-field stress from the strain when the constitutive equations remain
unknown. In contrast, the forward problem is defined as computing the stress
and strain fields from known constitutive relations and boundary conditions.

Inverse problem approaches target solutions consistent with force balance, the
observed deformation field, and an assumed constitutive equation form.
However, these constraints typically overdetermine the problem, so various
optimization procedures are developed to minimize error. For example, one can
minimize the difference between the assumed and computed deformation field
using finite element (FE) model updating [4, 5, 6], or error in the force-
balance equation using the virtual fields method [7, 8]. Other approaches
minimize the difference between two stress fields, one consistent with force
equilibrium, and one consistent with the assumed constitutive equation form
and observed strain field. These include the constitutive equation gap method
[9, 10], and the data-driven method recently introduced [11, 12, 13]. These
optimization-based methods can arrive at accurate solutions in the case of
homogeneous material properties, for example, Fletcher et al. [14] compute the
parameters of various isotropic elasto-plastic constitutive equations with
accuracies ranging between 0.07 and 5.51 %. However, accuracies are
inconsistent and depend on the specific problem, solution implementation, and
assumptions made, for example, Leygue et al. [11] arrive at an average stress-
field error of 8.5 % and maximum error greater than 50 % for a case of
homogeneous isotropic non-linear elasticity. Furthermore, these methods are
often computationally intensive (or unfeasible) due to the non-linearity of
the problem and the large number of iterations required. The problems of
accuracy and computational expense are exacerbated in heterogeneous materials,
where errors can be on the order of 10-50 % (e.g. [9, 15]). The reader can
refer to Pierron and Grédiac [16] for an extensive review of these approaches.

An alternative deterministic approach has recently been developed by the
authors [17], and subsequently by Liu [18], which allows one to compute the
stress field by solving a system of linear hyperbolic partial differential
equations. The key assumption used to derive these equations is that the
principal directions (eigenvectors) of the stress are known at each point from
the observed deformation. This is straightforward in the case of isotropic
(potentially non-linear) elasticity, or isotropic plastic deformation when the
elastic strain is negligible (see Section 2.1). Given this assumption, the
Cauchy stress, $\mathbf{T}$, will satisfy the following system of equations:

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

where $\mathbf{Q}$ is the matrix that diagonalizes $\mathbf{T}$, $\mathbf{b}$
is the body force, $\ddot{\mathbf{x}}$ is the acceleration, and $\rho$ is the
density111Notation: We use bold uppercase letters to denote second order
tensors ($\mathbf{A},\mathbf{B}$,…), bold lowercase letters to denote vectors
($\mathbf{a},\mathbf{b}$,…), unbolded letters to denote scalars ($a$, $b$,
$...$), and a superscript $\top$ to denote the transpose ($\mathbf{A}^{\top}$,
$\mathbf{a}^{\top}$, $...$). We use standard matrix notation, for example, we
have $\mathbf{a}\cdot\mathbf{b}=\mathbf{a}^{\top}\mathbf{b}$ and
$\mathbf{a}\otimes\mathbf{b}=\mathbf{a}\mathbf{b}^{\top}$. An over bar denotes
functional dependence, e.g. $\mathbf{A}=\bar{\mathbf{A}}(\mathbf{b})$. 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 of $\mathbf{A}$ with respect to the initial and deformed reference
frame respectively. $\mathrm{Grad}(\mathbf{A})$ and
$\mathrm{grad}(\mathbf{A)}$ correspond to the gradient of $\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}}^{\top}}$
where $\mathbf{Q_{\mathbf{A}}}$ has columns with the eigenvectors of
$\mathbf{A}$ denoted $\\{\mathbf{q}_{\mathbf{A}}^{(i)}|i=1,2,3\\}$ (principal
directions), and $\mathbf{\Lambda_{A}}$ is a diagonal matrix containing the
eigenvalues $\\{\lambda_{\mathbf{A}}^{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)\top}$.
As the spectral decomposition is used extensively for the stress we use
special notation:
$\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)\top}$ where
$\\{\sigma^{(i)}|i=1,2,3\\}$ are the principal stresses and
$\\{\mathbf{q}^{(i)}|i=1,2,3\\}$ are the principal directions. [17]. These
variables are generally a function of position and time, however, these
arguments will often be suppressed to avoid clutter (e.g.
$\mathbf{T}=\bar{\mathbf{T}}(\mathbf{x},t)$,
$\mathbf{Q}=\bar{\mathbf{Q}}(\mathbf{x},t)$, etc.). Eq. 1a is force balance
and Eq. 1b arises from the constraint on the principal directions of
$\mathbf{T}$. The basic workflow envisioned is (i) a deformation experiment is
conducted, (ii) the full-field kinematics are extracted using a technique such
as DIC, (iii) the principal directions of stress are extracted from the
kinematic fields, (iv) the stress field is obtained by solving the equations,
(v) useful information, such as material properties, are extracted in post-
processing. Step (iii) generally requires isotropy assumptions to determine
the principal directions of the stress and traction boundary conditions are
required to solve the equations. As discussed in Cameron and Tasan [17], this
system of equations can be applied to a wide range of materials and
deformation regimes including isotropic elastic solids, granular materials,
plastic solids consistent with co-directionality, a subset of other plastic
solids, Newtonian fluids, and a subset of non-Newtonian fluids. Furthermore,
these equations have been generalized to apply to isotropic elasto-plastic
deformation when the elastic strain is non-negligible (Cameron and Tasan,
submitted). One limitation is that these equations cannot be directly applied
to deformations with repeated stress eigenvalues. Addressing this will likely
be the subject of future work.

Before this approach can be widely applied, a robust numerical procedure must
be developed to solve the system of equations. While FE discretizations are
well suited to parabolic or elliptic systems of partial differential
equations, such as static linear elasticity, they are not well suited to
hyperbolic systems, such as Eq. 1. Hyperbolic systems describe wave-like or
transport phenomena, specifically, information propagates from the boundary
along characteristic curves throughout the domain. The variable(s) at each
point are influenced only by the variable(s) in the domain of dependence,
which is controlled by the characteristic curves and the direction of
information propagation. This imposes restrictions on the discretization, such
as the Courant–Friedrichs–Lewy condition ([20, chap. 6]), and often gives rise
to numerical instabilities if not properly accounted for.

Finite difference discretizations are commonly used for hyperbolic systems of
equations with a rectangular computational mesh (e.g. [21]). Indeed, in
Cameron and Tasan [17], a finite difference discretization was successfully
applied to solve Eq. 1 for a rectangular domain subject to infinitesimal
deformation. However, difficulties arise when applying finite difference
discretizations to solve Eq. 1 with an irregular computational mesh or domain
[20]. This is not just problematic for computing the stress field in
irregularly shaped objects, but also rectangular domains subject to finite
heterogeneous deformation. Additional problems arise because the direction of
the characteristic curves with respect to the chosen coordinates may change
throughout the domain. This is similar to the steady state advection equation
(e.g. [22]), but is in contrast to other systems such as elastic wave
propagation (e.g. Karal and Keller [23]). This means that the implicit and
explicit time stepping techniques, used for the latter case to resolve
instabilities, cannot generally be applied to Eq. 1. The finite difference
discretization applied in Cameron and Tasan [17] relied on the characteristic
curve orientations that were specific to that deformation, and becomes
unstable when applied to other problems.

Liu [18] developed a custom numerical procedure based on the method of
characteristics to solve the system of equations [18]. This addresses issues
of instability, and Liu demonstrates solutions in singularly-connected and
multiply-connected domains. However, the method is limited to cases with
geometric boundaries that are aligned with the characteristic curves of the
equation, imposing restrictions on the traction boundary conditions that may
be applied. Accuracy and computation time of this proposed approach needs to
be quantitatively assessed, as well.

Finite volume (FV) methods are well suited to hyperbolic systems with
irregular domains and unstructured meshes [24]. This research field is
relatively mature, and the discretization approach is regularly applied to
linear and non-linear problems in various applications [25, 26, 27, 28]. In
these applications, one utilizes the time dependence to give a stable solution
procedure. However, Eq. 1 is not time dependent, and there is no readily
available FV discretization that can be directly applied, hence, this
motivates developing such a method.

In order to develop a FV discretization, it is important to understand the
structure of the equation in the case of a discontinuity. This includes
continuous problems because the domain will be discretized into volume
elements with variables and coefficients that discontinuously change. Once the
behavior of the governing equation in the presence of these discontinuities is
deduced, a discretization can be formulated. Interestingly, reflection
phenomena arise at the discontinuities as the equations are second order and
hyperbolic.

Unlike when solving the forward problem, one cannot directly apply the known
traction boundary conditions when solving Eq. 1. For example, at a location
where all the characteristic lines are leaving the domain, the traction is
entirely determined as part of the solution [17]. Instead of directly
specifying the traction, one must translate knowledge of the traction and
characteristics into the boundary conditions. In this article we analyze the
general case and determine the equations that govern this.

This paper is structured as follows. In Section 2 we present analysis
addressing the characteristic curves, discontinuities, reflection, and
boundary conditions. In Section 3 we give analytical examples of an idealized
shear band and the Riemann problem, a common theoretical problem for
hyperbolic partial differential equations. Finally, in Section 4 we use the
results to develop a FV discretization which is validated against simulated
data.

##  2 Theory

###  2.1 Background

We review the governing equations and core assumptions which are presented in
Cameron and Tasan [17] or alternatively Liu [18]. Consider the simplest case
of infinitesimal isotropic elastic deformation that is potentially non-linear.
The spectral decomposition for the strain, $\mathbf{E}$, and the Cauchy
stress, $\mathbf{T}$, is given by:

$\mathbf{E}=\mathbf{Q_{E}\Lambda_{E}Q_{E}^{\top}},\qquad\mathbf{T}=\mathbf{Q_{T}\Lambda_{T}Q_{T}^{\top}}.$
(2)

The key assumption is that $\mathbf{Q}_{\mathbf{T}}=\mathbf{Q}_{\mathbf{E}}$,
which corresponds to the two tensors sharing the same eigenvectors. This is
consistent with isotropic elasticity including non-linear hyperelastic models.
We can therefore deduce the following linear relationship for $\mathbf{T}$

$[\mathbf{Q^{\top}_{E}TQ_{E}}]_{ij}=0\qquad\mathrm{for}\quad i\neq j.$ (3)

Intuitively, this corresponds to observing the stress tensor using a
coordinate system aligned with the principal directions of the strain and then
specifying that the shear components of the stress in this coordinate system
are zero. We assume that $\mathbf{Q}_{\mathbf{E}}$ and other kinematic
variables are known from experimental full-field deformation measurement,
parameters such as $\mathbf{b}$ and $\rho$ are known, and the boundary
conditions are specified. $\mathbf{T}$ will be the only unknown. When Eq. 3 is
combined with force balance, this gives a total of six equations and six
unknowns (Eq. 1). Hence, the solution is uniquely specified (with some
exceptions discussed in Section 2.2).

The traction boundary conditions are required over a subset of the domain
boundary - the areas of the boundary with incoming characteristic curves
(discussed later in this subsection). These boundary conditions are known to
be zero at free boundaries, and other boundary information can be obtained
through careful experimental design - either measuring the traction at the
boundaries or designing the problem in such a way that accurate assumptions
can be made. In the future, it may also be possible to calculate the boundary
conditions through optimization procedures and imposing various consistency
requirements on the solution.

The prior argument used to determine the principal direction of stress can be
applied to plastic deformation if one makes two key assumptions. First, that
$\mathbf{\dot{E}}^{e}<<\mathbf{\dot{E}}^{p}$. This may be practical in cases
where material elements are undergoing large strains (much larger than the
yield strain), although issues arise in heterogeneous deformations where
different regions are elastically and elasto-plastically deforming. Second,
that the principal directions of stress are aligned with
$\dot{\mathbf{E}}\approx\mathbf{\dot{E}}^{p}$. This is commonly used in
constitutive equations modeling isotropic plastic deformation, such as those
based on the co-directionality hypotheses. If these assumptions are met, one
can determine the principal directions directly from $\dot{\mathbf{E}}$.

These equations can be easily generalized to the case of finite deformation.
Before proceeding, we state the notation used. The position of a point in the
reference body is given by $\mathbf{X}$. The smooth one to one mapping from
the reference body to the deformed body is given by
$\mathbf{x}=\bar{\mathbf{x}}(\mathbf{X},t)$. The deformation gradient,
velocity, and velocity gradient are then given by

$\mathbf{F}=\mathrm{Grad}(\mathbf{x}),\qquad\mathbf{v}=\dot{\mathbf{x}},\qquad\mathbf{L}=\mathrm{grad}(\mathbf{v})=\dot{\mathbf{F}}\mathbf{F^{-1}}.$
(4)

The right and left Cauchy-Green deformation tensors are

$\mathbf{C}=\mathbf{U^{2}}=\mathbf{F^{\top}F},\qquad\mathbf{B}=\mathbf{V^{2}}=\mathbf{FF^{\top}}.$
(5)

$\mathbf{L}$ can be decomposed into its symmetric and skew parts

$\mathbf{D}=\frac{1}{2}[\mathbf{L}+\mathbf{L^{\top}}],\qquad\mathbf{W}=\frac{1}{2}[\mathbf{L}-\mathbf{L^{\top}}].$
(6)

In the case of isotropic elastic finite deformation, the principal directions
of stress align with the left Cauchy-Green deformation tensor $\mathbf{B}$.
And in the case of isotropic finite plastic deformation, where the elastic
strain is negligible, the principal directions align with $\mathbf{D}$. Again,
these assumptions are consistent with established constitutive equations for
large deformation elasticity and plasticity (e.g. [29]). Furthermore, the
alignment with $\mathbf{\dot{E}}$ or $\mathbf{D}$ is consistent with viscous
deformation in infinitesimal or finite deformation case respectively, so the
same approach can be applied to experimental fluid mechanics. Thus we have

$\displaystyle{}_{ij}=0\quad\textrm{for}\quad i\neq j$
$\displaystyle\textrm{infinitesimal elastic deformation},$ (7)
$\displaystyle[\mathbf{Q^{\top}_{\dot{E}}TQ_{\dot{E}}}]_{ij}=0\quad\textrm{for}\quad
i\neq j$ $\displaystyle\textrm{infinitesimal plastic and/or viscous
deformation},$
$\displaystyle[\mathbf{Q^{\top}_{B}TQ_{B}}]_{ij}=0\quad\textrm{for}\quad i\neq
j$ $\displaystyle\textrm{finite elastic deformation},$
$\displaystyle[\mathbf{Q^{\top}_{D}TQ_{D}}]_{ij}=0\quad\textrm{for}\quad i\neq
j$ $\displaystyle\textrm{finite plastic and/or viscous deformation}.$

These are all special cases of Eq. 1b. This system equations is linear,
hyperbolic, and variable coefficient as $\mathbf{Q}$ will vary as a function
of position. The boundary conditions will be discussed in Section 2.4. Note
that the equation and variables presented correspond to the deformed reference
frame. Note also that for the remainder of this paper we will assume
$\mathbf{Q}$, or equivalently the eigenvectors
$\\{\mathbf{q}^{(i)}|i=1,2,3\\}$, are known and we will not specify how they
are obtained, except for the case studies presented in Section 4.2. This
system of equations does not apply to cases with repeated stress eigenvalues
because $\mathbf{Q}$ will not be unique. For the remainder of this paper, we
will assume that the eigenvalues are unique, though this limitation will
likely be addressed in future work.

The two-dimensional form of this equation is important because full-field
deformation data is often only obtainable in two dimensions. Here, it can be
shown that Eq. 7 becomes equivalent to:

$\frac{T_{xx}-T_{yy}}{T_{xy}}=\frac{E_{xx}-E_{yy}}{E_{xy}}$ (8)

for the infintesimal elastic case where $\mathbf{T}$ shares eigenvectors with
$\mathbf{E}$ (equivalent expressions exist in the the plastic deformation
and/or finite deformation cases). This can be rearanged to a linear equation
in $\mathbf{T}$:

$T_{xx}-T_{yy}-cT_{xy}=0,$ (9)

with the known variable:

$c=\frac{E_{xx}-E_{yy}}{E_{xy}}.$ (10)

Eq. 9 is valid in the case of plane stress, or any other stress state where
the out-of-plane shear stress components are zero (such as an isotropic
material subject to plane strain).

We give a brief introduction to the characteristic curves to give the reader
intuition that is important for the remainder of the paper. Consider Fig. 1
which depicts a deformation in two dimensions where the principal directions
of the stress are known. Also, we assume the traction boundary conditions are
known. From knowledge of the principal directions, one can form a grid of
curved lines aligned with the principal directions at each point. These
curves, referred to as characteristic curves, will always intersect at right
angles because the principal directions of $\mathbf{T}$ are orthogonal.
Consider the material element shown in Fig. 1b. The boundaries of the element
are comprised of these characteristic curves, hence, the boundaries are
aligned with the principal directions of $\mathbf{T}$ and no shear force can
be exerted across the boundary, only normal force. If we know the magnitude of
the stresses on the left and bottom sides ($\sigma_{1}$ and $\sigma_{2}$), we
can use force balance in the $x$ and $y$ direction to compute the stresses on
the right and top sides ($\sigma_{3}$ and $\sigma_{4}$). To compute the stress
field throughout the entire domain, we can start on the lower left corner of
the domain, and proceed by moving upward and rightward. In the limit of small
material elements, the direction of computation will be along the
characteristic curves and one can think of this as information propagating
along these characteristics. One interesting feature is that the boundary
conditions on the top and right side will be computed as part of the solution,
rather than specified in the problem (as is the case with standard FE
simulations).

Figure 1: (a) A quasi-static deformation in two dimensions with the boundary
conditions specified on one side, and a traction free boundary on the top and
bottom. The dotted lines show the characteristic curves, which are aligned
with the principal directions of stress at each point (see text). (b) A
material element with boundaries defined by the characteristic curves. One can
use force balance to compute two of the stresses from two other stresses. When
this computation is conducted, information propagates from left to right and
bottom to top along the characteristic curves.

###  2.2 Characteristic curves and principal directions

Before discussing discontinuities and reflections, it is important to
mathematically show that the characteristic curves are indeed aligned with the
principal directions of the stress. Previously, this has simply been argued on
an intuitive basis [18, 17].

To show this, we carry out the following procedure. First, we combine the
system of differential equations and algebraic equations into one system of
differential equations. Second, we diagonalize this system and determine the
eigenvalues, which in turn specify the characteristic directions. Third, we
compare these characteristic directions to the principal directions of the
stress to show they are the same. We present equations and arguments in two
dimensions as the formulas are more concise and the relevant physics is the
same as the three-dimensional case.

We start with the following equations

$\begin{gathered}\frac{\partial T_{xx}}{\partial x}+\frac{\partial
T_{xy}}{\partial y}=0,\\\ \frac{\partial T_{xy}}{\partial x}+\frac{\partial
T_{yy}}{\partial y}=0,\\\ \end{gathered}$ (11)

and

$T_{xx}-T_{yy}-cT_{xy}=0.$ (12)

One can substitute for $T_{xx}$ in Eq. 11 using Eq. 12

$\begin{gathered}\frac{\partial T_{yy}}{\partial x}+c\frac{\partial
T_{xy}}{\partial x}+\frac{\partial T_{xy}}{\partial y}+T_{xy}\frac{\partial
c}{\partial x}=0,\\\ \frac{\partial T_{xy}}{\partial x}+\frac{\partial
T_{yy}}{\partial y}=0.\end{gathered}$ (13)

This can be written in matrix form:

$\frac{\partial\mathbf{u}}{\partial
x}+\mathbf{A}\frac{\partial\mathbf{u}}{\partial y}+\mathbf{B}\mathbf{u}=0,$
(14)

where,

$\mathbf{u}=\begin{bmatrix}T_{yy}\\\
T_{xy}\end{bmatrix},\qquad\mathbf{A}=\begin{bmatrix}0&1\\\
1&c\end{bmatrix},\qquad\mathbf{B}=\begin{bmatrix}0&0\\\ 0&\frac{\partial
c}{\partial x}\end{bmatrix}.$ (15)

According to standard theory on hyperbolic partial differential equations
(e.g. see [30]), the eigenvalues of $\mathbf{A}$,
$\\{\lambda_{\mathbf{A}}^{(\alpha)}|\alpha=1,2\\}$, determine the
characteristic curves of the equation, and the eigenvectors of $\mathbf{A}$,
$\\{\mathbf{q}_{\mathbf{A}}^{(\alpha)}|\alpha=1,2\\}$, give the information
that propagates along these characteristics. Note that we use Greek letters to
indicate the index when it takes values in the set $\\{1,2\\}$. One may
express the unknown $\mathbf{u}$ as a linear combination of the eigenvectors
of $\mathbf{A}$

$\mathbf{u}=a^{(1)}\mathbf{q}_{\mathbf{A}}^{(1)}+a^{(2)}\mathbf{q}_{\mathbf{A}}^{(2)}.$
(16)

Here, the $a^{(1)}\mathbf{q}_{\mathbf{A}}^{(1)}$ component of the stress will
propagate along the characteristic determined by $\lambda_{\mathbf{A}}^{(1)}$.
Similarly, the $a^{(2)}\mathbf{q}_{\mathbf{A}}^{(2)}$ component will propagate
along the characteristic determined by $\lambda_{\mathbf{A}}^{(2)}$. We have

$\begin{gathered}\lambda^{(1)}_{\mathbf{A}}=(c-\sqrt{c^{2}+4})/2,\qquad\lambda^{(2)}_{\mathbf{A}}=(c+\sqrt{c^{2}+4})/2,\\\
\mathbf{q}_{\mathbf{A}}^{(1)}=\begin{bmatrix}-(c+\sqrt{c^{2}+4})/2\\\
1\end{bmatrix},\qquad\mathbf{q}_{\mathbf{A}}^{(2)}=\begin{bmatrix}(-c+\sqrt{c^{2}+4})/2\\\
1\end{bmatrix}\end{gathered}.$ (17)

The vectors $\mathbf{v}^{(\alpha)}$ corresponding to the tangents to the
characteristic curves can then be computed using the equation defining the
$\alpha$ characteristic curve, $dx=\lambda^{(\alpha)}_{\mathbf{A}}dy$. This
gives

$\mathbf{v}^{(\alpha)}=\begin{bmatrix}\lambda_{\mathbf{A}}^{(\alpha)}\\\
1\end{bmatrix}.$ (18)

Hence we have

$\mathbf{v}^{(1)}=\begin{bmatrix}(c-\sqrt{c^{2}+4})/2\\\
1\end{bmatrix},\qquad\mathbf{v}^{(2)}=\begin{bmatrix}(c+\sqrt{c^{2}+4})/2\\\
1\end{bmatrix},$ (19)

where these vectors are defined up to a multiplicative constant. This gives
the spatial direction along which information propagates throughout the
material.

We now consider the components of stress that propagate along each
characteristic curve, which can be computed from
$a^{(1)}\mathbf{q}_{\mathbf{A}}^{(1)}$ and
$a^{(2)}\mathbf{q}_{\mathbf{A}}^{(2)}$ respectively. We label these
$\mathbf{T}^{(1)}$ and $\mathbf{T}^{(2)}$. Due to Eq. 16 and the definition of
$\mathbf{u}$ we have

$\mathbf{T}=\mathbf{T}^{(1)}+\mathbf{T}^{(2)},$ (20)

where $\mathbf{T}$ is the stress tensor in matrix form (in contrast to the
vector form $\mathbf{u}$). First we consider the stress matrix corresponding
to the component of $\mathbf{u}$ that propagates along $\alpha=1$
characteristic curve, given by $a^{(1)}\mathbf{q}_{\mathbf{A}}^{(1)}$. Using
the expression for $\mathbf{q}_{\mathbf{A}}^{(1)}$ we have
$T_{yy}^{(1)}=-a^{(1)}(c+\sqrt{c^{2}+4})/2$, and $T_{xy}^{(1)}=a^{(1)}$. Eq.
12 gives the $T_{xx}^{(1)}$ component. Hence, we have

$\begin{gathered}\mathbf{T}^{(1)}=a^{(1)}\begin{bmatrix}(c-\sqrt{c^{2}+4})/2&1\\\
1&-(c+\sqrt{c^{2}+4})/2\end{bmatrix}.\end{gathered}$ (21)

This is a tensor with matrix rank one, as one of the eigenvalues are zero
(this should not be confused with a vector, or first-order tensor, which is
often referred to as a rank one tensor). The non-zero eigenvalue and
corresponding eigenvector of $\mathbf{T}^{(1)}$ is

$\begin{gathered}\lambda_{\mathbf{T}(1)}=-a^{(1)}\sqrt{c^{2}+4},\\\
\mathbf{q}_{\mathbf{T}(1)}=\begin{bmatrix}(c-\sqrt{c^{2}+4})/2\\\
1\end{bmatrix},\end{gathered}$ (22)

where the $(1)$ is in the subscript position to indicate that it corresponds
to the $\mathbf{T}^{(1)}$ matrix. One can compare this expression with Eq. 19
to conclude that the eigenvector of the stress matrix $\mathbf{T}^{(1)}$ is
the same as $\mathbf{v}^{(1)}$, and gives the characteristic direction for the
$\mathbf{T}^{(1)}$ component of the stress.

One can repeat the same exercise for the $\alpha=2$ characteristic curve

$\begin{gathered}\mathbf{T}^{(2)}=a^{(2)}\begin{bmatrix}(c+\sqrt{c^{2}+4})/2&1\\\
1&-(c-\sqrt{c^{2}+4})/2\end{bmatrix},\end{gathered}$ (23)
$\begin{gathered}\lambda_{\mathbf{T}(2)}=-a^{(2)}\sqrt{c^{2}+4},\\\
\mathbf{q}_{\mathbf{T}(2)}=\begin{bmatrix}-(c-\sqrt{c^{2}+4})/2\\\
1\end{bmatrix}.\end{gathered}$ (24)

Again, this shows that the principal direction of the stress is coaxial with
the respective characteristic curve.

One can substitute the spectral decompositions for each
$\mathbf{T}^{(\alpha)}$ into Eq. 20, noting that each $\mathbf{T}^{(\alpha)}$
is matrix rank one

$\mathbf{T}=\lambda_{\mathbf{T}(1)}\mathbf{q}_{\mathbf{T}(1)}\mathbf{q}^{\top}_{\mathbf{T}(1)}+\lambda_{\mathbf{T}(2)}\mathbf{q}_{\mathbf{T}(2)}\mathbf{q}^{\top}_{\mathbf{T}(2)}.$
(25)

The spectral decomposition for $\mathbf{T}$ in two dimensions can be expressed
as

$\mathbf{T}=\sigma^{(1)}\mathbf{q}^{(1)}\mathbf{q}^{(1)\top}+\sigma^{(2)}\mathbf{q}^{(2)}\mathbf{q}^{(2)\top},$
(26)

where $\\{\sigma^{(\alpha)}\\}$ are the principal stresses and
$\\{\mathbf{q}^{(\alpha)}\\}$ (without the subscript) are the principal
directions.

Comparing Eqs. 25 and 26, one can take
$\lambda_{\mathbf{T}(\alpha)}=\sigma^{(\alpha)}$ and
$\mathbf{q}_{\mathbf{T}(\alpha)}=\mathbf{q}^{(\alpha)}$. Hence, the stress can
simply be decomposed into two matrix-rank one tensors using the spectral
decomposition, and each tensor will propagate along the direction of its non-
zero principal direction. We can interpret this as the principal stresses
simply propagating along their respective principal directions. This is the
key result of the argument just made and allows one to make a clear physical
interpretation of the characteristic curves.

We note that the same argument can be made in three dimensions, leading to the
same conclusion: that the principal directions of stress are aligned with the
respective characteristics of the governing equation. Hence, it will
frequently be convenient to decompose $\mathbf{T}$ into its eigenvalues and
eigenvectors

$\mathbf{T}=\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)\top}.$ (27)

We make several remarks:

  * 1.

One must assign the direction of information propagation along each
characteristic curve when solving the system of equations. This gives
$2^{2}=4$ possibilities in two dimensions and $2^{3}=8$ possibilities in three
dimensions, one binary choice for each characteristic curve. The choice of the
direction of information propagation is intrinsically linked to which boundary
conditions are specified, as the principal stresses can only be specified
along the characteristics entering the domain.

  * 2.

The stress at a particular point is only determined by the stress in the
domain of dependence (Fig. 2). The stress at a particular point influences
stresses only in its domain of influence. We refer to upstream and downstream
as positions along characteristic curves reverse to the direction of
information propagation and along the direction of information propagation
respectively. Points in the domain of dependence can be reached by following
any combination of characteristic lines upstream from the point of interest.
Similarly, the points in the domain of influence can be found by following any
combination of characteristics in the downstream direction.

  * 3.

The characteristic curves are always perpendicular to one another when they
intersect. In two dimensions, the variable coefficient $c$ controls the
orientation of the characteristic curves, however, in three dimensions, three
independent variables are required. We will often use $\\{\mathbf{q}^{(i)}\\}$
in two or three dimensions to specify the orientation of the characteristic
lines. In some parts of the paper, it will be more concise to refer to the
principal direction orientation angle with respect to the $x$ axis, $\theta$.
The other principal direction will be oriented at $\theta+\pi/2$. We have the
expression

$\theta=\tan^{-1}\left(\frac{1}{c+\sqrt{1+c^{2}}}\right)$ (28)

  * 4.

There are several forms of the governing equation which may be useful in
different applications. For example, one equation form can be obtained in two
dimensions using the stress Airy function

$\frac{\partial\psi}{\partial y^{2}}-\frac{\partial\psi}{\partial
x^{2}}+c\frac{\partial\psi}{\partial x\partial y}=0.$ (29)

This can be decomposed into a system of first-order equations

$\begin{gathered}\frac{\partial\psi}{\partial
y}+\frac{1}{2}(c+\sqrt{c^{2}+4})\frac{\partial\psi}{\partial x}=0,\\\
\frac{\partial\psi}{\partial
y}+\frac{1}{2}(c-\sqrt{c^{2}+4})\frac{\partial\psi}{\partial
x}=0,\end{gathered}$ (30)

where we assume there is no body force or acceleration term. This particular
form may be useful in some applications as there is no source term.

  * 5.

While in many cases the solution is deterministic, non-uniqueness arises when
the characteristic curves form closed loops. For example, the case of a
rotationally symmetric expansion depicted in Fig. 2c. This will likely the
subject of future work.

Figure 2: (a) The domain of dependence and influence in two dimensions where
the characteristic directions are varying with position. (b) The domain of
dependence and influence in three dimensions where the characteristic
directions are constant. (c) An example where the solution is non-unique as
the characteristic lines form closed loops. This example shows a rotationally
symmetric problem where a hole exerts a constant pressure on the domain
boundary.

###  2.3 Discontinuities

Next, we consider a simple discontinuity in the principal directions of the
stress. This could arise from discontinuities in the observed strain or strain
rate field, however, it also arises in the finite volume discretization
process even when the strain field is continuous. Hence, it is necessary to
understand prior to developing the discretization procedure. We will show how
discontinuities give rise to reflection and transmission phenomena. For
simplicity, we assume the material is in static equilibrium and that the
problem is two dimensional (we comment at the end of this subsection on the
three-dimensional case). We also note that one can solve these problems
trivially by balancing traction across the interface and utilizing Eq. 9,
however, it is insightful to show how the same solution arises directly from
the system of partial differential equations. We will then utilize the simpler
traction balance approach for the remainder of the paper.

First, consider the problem depicted in Fig. 3a. Here, we assume the stress is
known for $x<0$, and must be computed for $x>0$. The coefficient $c$ is given
by

$c=\left.\begin{cases}c_{L},&\text{if }x\leq 0,\\\ c_{R},&\text{if }x>0,\\\
\end{cases}\right.$ (31)

where the subscripts $L$ and $R$ correspond to left and right respectively.
The unit vectors $\\{\mathbf{q}^{(\alpha)}_{L}\\}$ and
$\\{\mathbf{q}^{(\alpha)}_{R}\\}$ are determined by $c_{L}$ and $c_{R}$
respectively using Eqs. 22 and 24. Here the solution will not depend on $y$,
hence, we eliminate the $y$ derivatives from Eq. 13 giving

$\begin{gathered}\frac{\partial T_{yy}}{\partial x}+c\frac{\partial
T_{xy}}{\partial x}+T_{xy}\frac{\partial c}{\partial x}=0,\\\ \frac{\partial
T_{xy}}{\partial x}=0.\end{gathered}$ (32)

This directly gives:

$T_{xyL}=T_{xyR}=T_{xy},$ (33)

and

$\frac{\partial T_{yy}}{\partial x}+T_{xy}\frac{\partial c}{\partial x}=0.$
(34)

We have the expression for $\frac{\partial c}{\partial x}$

$\frac{\partial c}{\partial x}=(c_{R}-c_{L})\delta(x),$ (35)

where $\delta(x)$ is the Dirac delta function. Substituting this into Eq. 34
then integrating we have:

$\begin{gathered}\frac{\partial T_{yy}}{\partial
x}=-T_{xy}(c_{R}-c_{L})\delta(x),\\\ T_{yyR}-T_{yyL}=-T_{xy}(c_{R}-c_{L}),\\\
T_{yyR}=T_{xy}(c_{L}-c_{R})+T_{yyL}.\end{gathered}$ (36)

Finally, by calculating $T_{xxL}$ and $T_{xxR}$ using Eq. 12 one can obtain:

$T_{xxL}=T_{xxR}.$ (37)

Figure 3: Solving the equations over a discontinuity for different
characteristic directions. On the left and right hand side of the
discontinuity the principal directions, or characteristics, have different
orientations. (a) The stress on one side of the discontinuity is known. (b)
One principal stress is known on each side. This latter case leads to
reflection and transmission of the principal stresses.

Together, Eqs. 33, 36 and 37 give expressions for the stress where $x>0$ in
terms of the stress where $x\leq 0$. This solution is simply equivalent to
balancing the traction across the discontinuity and specifying the principal
directions.

One can repeat the same exercise for a $y=\mathrm{constant}$ discontinuity,
where force balance plus the constraint on the principal directions is
sufficient to determine the unknown stress. We can express the general formula
relating the stress on each side of a discontinuity as

$\mathbf{T}_{A}\mathbf{n}=\mathbf{T}_{B}\mathbf{n},$ (38)

where $\mathbf{n}$ is the vector normal to the discontinuity, and A and B
correspond to two regions on each side of the discontinuity. In two or three
dimensions, this can be expressed in the following equations relating
$\\{\sigma^{(i)}_{A},\sigma^{(i)}_{B}\\}$

$\sum_{i}\sigma_{A}^{(i)}\mathbf{q}_{A}^{(i)}\mathbf{q}^{(i)\top}_{A}\mathbf{n}=\sum_{i}\sigma_{B}^{(i)}\mathbf{q}_{B}^{(i)}\mathbf{q}^{(i)\top}_{B}\mathbf{n}.$
(39)

It may not be true that the stress is fully determined on one side of the
discontinuity and unknown on the other: it may be partially determined on each
side. For the general case in two dimensions, two characteristics are incoming
to the boundary, and two are outgoing from the boundary. In three dimensions,
there are three characteristics incoming and three outgoing. Given that
$\\{\sigma^{(i)}\\}$ along the incoming characteristics is assumed known,
$\\{\sigma^{(i)}\\}$ along the outgoing characteristics can be computed using
Eq. 39.

We briefly consider the case where characteristics are parallel to the
discontinuity. In two dimensions, where one characteristic curve on one side
of the discontinuity is parallel to the discontinuity, two independent force
balance equations will remain. Therefore, there should be two characteristics
outgoing from the discontinuity and one incoming so the solution is unique. In
three dimensions, there are several possibilities which depend on the number
of characteristics parallel to the discontinuity. Eq. 39 reduces to a
different number of equations depending on the scenario. For example, in the
case of two characteristics on each side being parallel to the boundary, Eq.
39 reduces to one equation, so there should be one outgoing characteristic.

One can understand this analysis in terms of the ”reflection” and
”transmission” of the principal stresses. In general, reflection and
transmission phenomena arise whenever the system is governed by a set of
partial differential equations that are linear, hyperbolic, second order or
higher, variable coefficient, and discontinuous in the variable
coefficient(s). To give the reader intuition for this, we consider the analogy
of a sound wave. The mathematics governing this refection is very similar to
the mathematics governing the behavior of the principal stresses at the
discontinuity. Specifically, a discontinuity in the principal directions of
stress is analogous to the discontinuity in sound wave speed, as the wave
speed is essentially a ”direction” in the space-time domain. In addition, the
principal stresses themselves are analogous to the amplitudes of sound waves
propagating in different directions. When a sound wave encounters a wave speed
discontinuity, a component reflects off the discontinuity and another
transmits through it, and the amplitudes of each component are determined by
the respective wave speeds and Newton’s second law. In the principal stress
case, the principal stress propagating along the direction incoming toward the
discontinuity is converted into two principal stresses propagating along
directions outgoing from the discontinuity. We can say that a component of the
principal stress is ”reflected” and a component is ”transmitted”. (The reader
can refer to the example in Section 3.2 which may give further intuition).

Consider a principal stress $\sigma^{(1)}_{L}$ on an incoming characteristic,
$\mathbf{q}^{(1)}_{L}$, as shown in Fig. 3b, where we assume
$\sigma^{(2)}_{R}=0$. The incoming principal stress, $\sigma^{(1)}_{L}$, will
have a reflected component along $\mathbf{q}^{(2)}_{L}$ of magnitude
$\sigma^{(2)}_{L}=r_{r}\sigma^{(1)}_{L}$, and a transmitted component along
$\mathbf{q}^{(2)}_{R}$ of magnitude $\sigma^{(2)}_{R}=r_{t}\sigma^{(1)}_{L}$.
To calculate $r_{r}$ and $r_{t}$ we take $\sigma^{(1)}_{L}=1$ and use Eq. 39.
One arrives two equations and two unknowns that can be solved

$r_{t}\mathbf{q}_{R}^{(1)}\mathbf{q}^{(1)\top}_{R}\mathbf{n}-r_{r}\mathbf{q}^{(2)}_{L}\mathbf{q}^{(2)\top}_{L}\mathbf{n}=\mathbf{q}^{(1)}_{L}\mathbf{q}^{(1)\top}_{L}\mathbf{n}.$
(40)

Hence, $r_{r}$ and $r_{t}$ are determined by the principal directions
$\\{\mathbf{q}^{(\alpha)}_{L},\mathbf{q}^{(\alpha)}_{R}\\}$, and boundary
normal $\mathbf{n}$. Substituting for
$\\{\mathbf{q}^{(\alpha)}_{L},\mathbf{q}^{(\alpha)}_{R}\\}$, and taking
$\mathbf{n}$ to be in the $x$ direction gives

$r_{t}=\frac{\cos(\theta_{L})}{\cos(\theta_{L}-\theta_{R})\cos(\theta_{R})},\qquad
r_{r}=\frac{\tan(\theta_{L}-\theta_{R})}{\tan(\theta_{L})},$ (41)

where we express the answer using $\theta_{L}$ and $\theta_{R}$ defined in
Section 2.2). Finally, when considering the case of two incoming principal
stresses, one can simply apply the principal of superposition: the
$\sigma^{(i)}$ on an outgoing characteristic can be expressed as the sum of
values each arising from the reflection or transmission of an incoming
characteristic. This same approach can be applied directly to the three-
dimensional case.

###  2.4 Boundary conditions

This section discusses how knowledge of the traction at the boundary can be
translated to specify the boundary condition. Specifically, we must specify
the principal stresses along the characteristics entering the domain, and not
specify the principal stresses along characteristics leaving the domain.

Here, we utilize the results from the prior subsections showing the
characteristics are aligned with the principal directions of the stress, and
solving the equations at a discontinuity simply corresponds to balancing the
traction. In general, $\\{\sigma^{(i)}\\}$ must be specified on the
characteristics entering the domain, $\\{\mathbf{q}^{(ei)}\\}$, but cannot be
specified on those leaving the domain, $\\{\mathbf{q}^{(li)}\\}$. If
$\mathbf{n}$ is the unit outward normal of the domain boundary, we have
$\mathbf{q}^{(li)\top}\mathbf{n}>0$ and $\mathbf{q}^{(ei)\top}\mathbf{n}<0$.
We have the additional special case of characteristics aligned with the
boundary, $\\{\mathbf{q}^{(ai)}\\}$, where
$\mathbf{q}^{(ai)\top}\mathbf{n}=0$. We have
$\\{\mathbf{q}^{(i)}\\}=\\{\mathbf{q}^{(ei)}\\}\cup\\{\mathbf{q}^{(li)}\\}\cup\\{\mathbf{q}^{(ai)}\\}$.
There will be $n$ characteristics at each point, where $n=3$ in three
dimensions and $n=2$ in two dimensions. Let $n^{e}$ denote the number of
characteristics entering the domain, $n^{l}$ denote the number of
characteristics leaving the domain, and $n^{a}$ denote the number of
characteristics aligned with the domain boundary such that
$n^{e}+n^{l}+n^{a}=n$.

In general, we will utilize the equation

$\sum_{i}\sigma^{(i)}\mathbf{q}^{(i)}\mathbf{q}^{(i)\top}\mathbf{n}=\mathbf{t},$
(42)

where $\mathbf{t}$ is the traction at the boundary. In the straightforward
case where all characteristics are entering the domain ($n^{e}=n$), Eq. 42
gives three equations and three unknowns which can be solved for
$\\{\sigma^{(i)}\\}$. Similarly, in the case where all characteristics are
leaving the domain ($n^{l}=n$), $\mathbf{t}$ is entirely determined from the
values in the domain so no boundary traction can be specified.

In the case where $n^{a}\neq 0$, the term(s) containing
$\mathbf{q}^{(ai)\top}\mathbf{n}$ in Eq. 42 become zero eliminating the
unknown(s) $\sigma^{(ai)}$ from the equation. Hence, the system of equations
will become over-determined for a general $\mathbf{t}$. In order to resolve
this issue, one must specify a $\mathbf{t}$ such that
$\mathbf{q}^{(ai)\top}\mathbf{t}=0$. Physically, this corresponds to the
traction on a surface having no shear component when a principal direction of
the stress is orthogonal to the surface.

For any value of $n^{e}$, one can specify $n^{e}$ components of the traction.
In the case where $0<n^{e}<n$, one can also choose which components of the
traction are specified. Let $\\{\mathbf{d}^{(j)}\\}$ be any orthonormal basis
spanning the two- or three-dimensional space. One may choose a subset of
$\\{\mathbf{d}^{(j)}\\}$ of size $n^{e}$ in which to specify the traction
components. We denote this $\\{\mathbf{d}^{(j)*}\\}$. Taking the dot product
of Eq. 42 with $\mathbf{d}^{(j)*}$ gives:

$\sum_{i}\sigma^{(i)}(\mathbf{d}^{(j)*\top}\mathbf{q}^{(i)})(\mathbf{q}^{(i)\top}\mathbf{n})=\mathbf{d}^{(j)*\top}\mathbf{t}=t_{j}.$
(43)

This is a scalar equation involving $\\{\sigma^{(i)}\\}$. There will be
$n^{e}$ equations and $n^{e}$ unknowns, so the equations can be solved. One is
free to select any $n^{e}$ orthonormal vectors $\\{\mathbf{d}^{(i)*}\\}$, with
the only restriction being that $\mathbf{d}^{(i)}\neq\pm\mathbf{q}^{(lj)}$,
otherwise the equations for $\sigma^{(ei)}$ will be singular. Freedom to
select $\\{\mathbf{d}^{(i)*}\\}$ allows one to:

  * 1.

specify boundary conditions that are non-reflective by choosing
$\\{\mathbf{d}^{(j)*}\\}=\\{\mathbf{q}^{(ei)}\\}$

  * 2.

specify the orientation of $\mathbf{t}$ by choosing particular $t_{j}=0$

  * 3.

specify the magnitude of $\mathbf{t}$ in certain directions by choosing
particular $t_{j}$ and $\mathbf{d}^{(j)*}$.

Figure 4: Determining the boundary conditions when one characteristic is
entering the domain, $\mathbf{q}^{(e)}$, and another is leaving,
$\mathbf{q}^{(l)}$. One can specify the traction components in chosen
$\mathbf{d^{(1)*}}$ directions.

We consider an example in two dimensions of one characteristic entering the
domain and one characteristic leaving the domain. Here, one can specify one
component of the traction vector, whereas the other component is determined by
the principal stress along the leaving characteristic (Fig. 4). Let
$\mathbf{d}^{(1)*}$ be a unit vector corresponding to the direction where one
chooses to specify the traction, and $\mathbf{d}^{(2)}$ correspond to the
direction where $\mathbf{t}$ will be determined. Let
$t_{1}=\mathbf{d}^{(1)*\top}\mathbf{t}$ and
$t_{2}=\mathbf{d}^{(2)\top}\mathbf{t}$. Let $\mathbf{q}^{(e)}$ correspond to
the characteristic entering the domain and $\mathbf{q}^{(l)}$ correspond to
the characteristic leaving the domain. One must determine $\sigma^{(e)}$,
where $\sigma^{(l)}$ is known. We have

$\sigma^{(1)}\mathbf{q}^{(1)}\mathbf{q}^{(1)\top}\mathbf{n}+\sigma^{(e)}\mathbf{q}^{(e)}\mathbf{q}^{(e)\top}\mathbf{n}=\mathbf{t}.$
(44)

We take the dot product of this equation with $\mathbf{d}^{(1)*}$ and
$\mathbf{d}^{(2)}$ respectively giving

$\begin{gathered}\sigma^{(l)}\mathbf{d}^{(1)*\top}\mathbf{q}^{(l)}\mathbf{q}^{(l)\top}\mathbf{n}+\sigma^{(e)}\mathbf{d}^{(1)*\top}\mathbf{q}^{(e)}\mathbf{q}^{(e)\top}\mathbf{n}=t_{1},\end{gathered}$
(45)
$\sigma^{(l)}\mathbf{d}^{(2)\top}\mathbf{q}^{(l)}\mathbf{q}^{(l)\top}\mathbf{n}+\sigma^{(e)}\mathbf{d}^{(2)}\mathbf{q}^{(e)}\mathbf{q}^{(e)\top}\mathbf{n}=t_{2}.$
(46)

rearranging Eq. 45 gives:

$\sigma^{(e)}=\frac{t_{1}-\sigma^{(l)}\mathbf{d}^{(l)*\top}\mathbf{q}^{(l)}\mathbf{q}^{(l)\top}\mathbf{n}}{\mathbf{d}^{(l)*\top}\mathbf{q}^{(e)}\mathbf{q}^{(e)\top}\mathbf{n}}.$
(47)

This gives the principal component for the entering characteristic. One can
interpret this as being comprised of two terms, one due to the specified
component of traction along the boundary, and another due to reflection off
the boundary from the leaving characteristic. We also obtain the expression
for the unknown component of the traction by substituting this into Eq. 46

$t_{2}=\sigma^{(l)}\mathbf{d}^{(2)\top}\mathbf{q}^{(l)}\mathbf{q}^{(l)\top}\mathbf{n}+\frac{(t_{1}-\sigma^{(l)}\mathbf{d}^{(1)*\top}\mathbf{q}^{(l)}\mathbf{q}^{(l)\top}\mathbf{n})\mathbf{d}^{(2)\top}\mathbf{q}^{(e)}}{\mathbf{d}^{(1)*\top}\mathbf{q}^{(e)}}.$
(48)

One has freedom to select $\mathbf{d}^{(1)*}$ except for the restriction
mentioned above. To eliminate reflection at the boundary, one may specify
$\mathbf{d}^{(1)*}=\mathbf{q}^{(e)}$ and $\mathbf{d}^{(2)}=\mathbf{q}^{(l)}$.
This gives

$\begin{gathered}\sigma^{(e)}=\frac{t_{1}}{\mathbf{q}^{(e)\top}\mathbf{n}},\\\
t_{2}=\sigma^{(l)}\mathbf{q}^{(l)\top}\mathbf{n}.\end{gathered}$ (49)

We remark on a few boundaries of physical relevance. First, consider a free
boundary that will have zero traction. Here, one characteristic is orthogonal
to the boundary, and the other characteristic(s) will align with it. In the
case where the orthogonal characteristic is entering the domain, one can
specify that the normal component of the traction is zero. Second, consider a
domain boundary that is in contact with a frictionless rigid surface. Here,
there will also be one characteristic orthogonal to the boundary and the other
characteristic(s) will be aligned with it. In the case where one or two
characteristics are entering the domain, one may specify the
$\mathbf{d}^{(i)*}$ to be orthogonal to $\mathbf{n}$ and the corresponding
traction in these direction(s) to be zero, hence, the traction vector will be
normal to the boundary with a magnitude that is determined instead of
specified. Third, consider the case where $\mathbf{t}$ is known (or assumed).
One can directly apply the equations discussed, however, in the case of one
characteristic entering and one leaving, one must choose the direction(s)
$\\{\mathbf{d}^{(i)*}\\}$ along which to specify the component of
$\mathbf{t}$.

##  3 Analytical examples

###  3.1 Riemann problem

We now consider the solution to the Riemann problem in two dimensions, which
is commonly used to perform computations for various systems of hyperbolic
partial differential equations.

The problem is depicted in Fig. 5. We have a discontinuous boundary traction
and a discontinuous coefficient $c$ that both change at $x=0$. We have
characteristic directions $\\{\mathbf{q}^{(\alpha)}_{L}\\}$ on the left side
($x\leq 0$), and $\\{\mathbf{q}^{(\alpha)}_{R}\\}$ on the right side ($x>0$).
Similarly the traction $\mathbf{t}_{L}$ is applied for $x\leq 0$, and
$\mathbf{t}_{R}$ is applied for $x>0$. We specify the direction of information
propagation such that the information propagates in the positive $y$
direction. We also assume the material is in static equilibrium. The stress
will have four different states, one in each region A, B, C and D. Regions A-B
boundary will be aligned with the $\mathbf{q}^{(2)}_{L}$ characteristic and
the C-D boundary will be aligned with the $\mathbf{q}^{(1)}_{R}$
characteristic. The normal vector to the domain boundary is denoted
$\mathbf{n}^{*}$ and the normal vector to the B-C boundary is denoted
$\mathbf{n}$. We have eight unknowns:
$\\{\sigma^{(1)}_{A},\sigma^{(2)}_{A},\sigma^{(1)}_{B},\sigma^{(2)}_{B},\sigma^{(1)}_{C},\sigma^{(2)}_{C},\sigma^{(1)}_{D},\sigma^{(2)}_{D}\\}$.

Figure 5: Schematic of the Riemann problem for the governing equations.

One can determine
$\\{\sigma^{(1)}_{A},\sigma^{(2)}_{A},\sigma^{(1)}_{D},\sigma^{(2)}_{D}\\}$
using force balance and the constraint on the principal directions as
discussed in Section 2.4

$\begin{gathered}\sigma^{(1)}_{A}\mathbf{q}_{L}^{(1)}\mathbf{q}_{L}^{(1)\top}\mathbf{n}^{*}+\sigma^{(2)}_{A}\mathbf{q}^{(2)}_{L}\mathbf{q}_{L}^{(2)\top}\mathbf{n}^{*}=\mathbf{t}_{L},\\\
\sigma_{D}^{(1)}\mathbf{q}_{R}^{(1)}\mathbf{q}_{R}^{(1)\top}\mathbf{n}^{*}+\sigma_{D}^{(2)}\mathbf{q}_{R}^{(2)}\mathbf{q}_{R}^{(2)\top}\mathbf{n}^{*}=\mathbf{t}_{R}.\end{gathered}$
(50)

One can directly solve these equations as they are two linear equations with
two unknowns. Next, consider the A-B boundary. There is no change in
$\\{\mathbf{q}^{(\alpha)}_{L}\\}$, hence, we have

$\sigma_{A}^{(1)}=\sigma_{B}^{(1)}.$ (51)

Similarly, for the C-D boundary, we have

$\sigma_{C}^{(2)}=\sigma_{D}^{(2)}.$ (52)

Finally, force balance is used along the B-C boundary to obtain
$\sigma_{B}^{(2)}$ and $\sigma_{C}^{(1)}$ (see Section 2.3)

$\sigma^{(1)}_{B}\mathbf{q}^{(1)}_{L}\mathbf{q}^{(1)\top}_{L}\mathbf{n}+\sigma^{(2)}_{B}\mathbf{q}^{(2)}_{L}\mathbf{q}^{(2)\top}_{L}\mathbf{n}=\sigma^{(1)}_{C}\mathbf{q}^{(1)}_{R}\mathbf{q}^{(1)\top}_{R}\mathbf{n}+\sigma^{(2)}_{C}\mathbf{q}^{(2)}_{R}\mathbf{q}^{(2)\top}_{R}\mathbf{n}.$
(53)

The solution of this equation can be expressed using the $\theta$ notation
(Section 2.2) to give expressions for the unknown $\sigma_{B}^{(2)}$ and
$\sigma_{C}^{(1)}$

$\sigma_{B}^{(2)}=\frac{\sigma_{A}^{(1)}\,\sin\left(2\,\theta_{L}-\theta_{R}\right)-\sigma_{A}^{(1)}\,\sin\left(\theta_{R}\right)+2\,\sigma_{D}^{(2)}\,\sin\left(\theta_{R}\right)}{\sin\left(2\,\theta_{L}-\theta_{R}\right)+\sin\left(\theta_{R}\right)},$
(54)

and

$\sigma_{C}^{(1)}=-\frac{2\,\sigma_{A}^{(1)}\,\cos\left(\theta_{L}\right)-\sigma_{D}^{(2)}\,\cos\left(\theta_{L}\right)+\sigma_{D}^{(2)}\,\cos\left(\theta_{L}-2\,\theta_{R}\right)}{\cos\left(\theta_{L}-2\,\theta_{R}\right)+\cos\left(\theta_{L}\right)}.$
(55)

Eqs. 50, 51, 52 and 53 fully determine the $\sigma^{(\alpha)}$ in every region
and solves the problem.

###  3.2 Idealized shear band

We next consider the analytical example of an idealized shear band (Fig. 6a).
Here, a sample is under tension with a uniform horizontal traction of
magnitude $t$ on the constant $x$ boundaries (taking $x$ to be in the
horizontal direction and $y$ to be in the vertical direction). The shear band
is oriented at $45^{o}$ to the horizontal, i.e.
$\mathbf{n}=[-1/\sqrt{2},-1/\sqrt{2}]^{\top}$ (Fig. 6a). The shear band also
has an associated angle $\beta$. We have regions A, B, D, and E outside the
shear band and region C inside the shear band. We assume that the principal
directions of the stress are aligned with the $x$ and $y$ axis outside the
shear band so we have

$\begin{split}\mathbf{q}^{(1)}_{A}=\mathbf{q}^{(1)}_{B}=\begin{bmatrix}0\\\
-1\end{bmatrix},\qquad\mathbf{q}^{(2)}_{A}=\mathbf{q}^{(2)}_{B}=\begin{bmatrix}1\\\
0\end{bmatrix},\\\
\qquad\mathbf{q}^{(1)}_{D}=\mathbf{q}^{(1)}_{E}=\begin{bmatrix}1\\\
0\end{bmatrix},\qquad\mathbf{q}^{(2)}_{D}=\mathbf{q}^{(2)}_{E}=\begin{bmatrix}0\\\
1\end{bmatrix}.\end{split}$ (56)

The signs are chosen so that the stress can be computed in all regions
starting from the boundary condition on the left side. We also number
$\\{\mathbf{q}^{(i)}\\}$ such that $\mathbf{q}^{(2)}$ is a $90^{o}$
counterclockwise rotation from $\mathbf{q}^{(1)}$.

Figure 6: An idealized shear band which demonstrates the reflection and
transmission of the principal stress. (a) The different regions, geometry and
boundary conditions. (b) The reflection of the principal stress off the shear
band, and the propagation through the shear band. The transmission coefficient
gives the stresses in the shear band, and the reflection coefficient enables
computation of the boundary conditions on the upper and lower side. The red
lines are the characteristics.

We assume that the stress state in the shear band is constant. Hence, due to
the free boundary on the top and bottom of the shear band, we can deduce that

$\mathbf{q}^{(1)}_{C}=\begin{bmatrix}-\sin(\beta)\\\
-\cos(\beta)\end{bmatrix},\qquad\mathbf{q}^{(2)}_{C}=\begin{bmatrix}\cos(\beta)\\\
-\sin(\beta)\end{bmatrix}.$ (57)

Further, we have that $\sigma_{C}^{(1)}=0$ because of the traction-free
boundary on the top and bottom of the shear band (assuming the stress in the
shear band is constant).

Beginning at the boundary on the left side, we have

$\sigma_{A}^{(2)}=\sigma_{B}^{(2)}=t.$ (58)

$\sigma_{A}^{(1)}$ will be governed by the traction applied to the upper
boundary of region A, which we assume to be zero. Hence we have

$\sigma_{A}^{(1)}=0.$ (59)

Consider the reflection of the $\sigma_{B}^{(2)}$ principal stress at the B-C
interface. $\mathbf{q}_{B}^{(2)}$ and $\mathbf{q}_{C}^{(1)}$ will be incoming
to the interface with known principal stresses, and that
$\mathbf{q}_{B}^{(1)}$ and $\mathbf{q}_{C}^{(2)}$ will be outgoing from the
interface with unknown principal stresses. We solve the following equations to
obtain these stresses

$\sigma^{(1)}_{B}\mathbf{q}^{(1)}_{B}\mathbf{q}^{(1)\top}_{B}\mathbf{n}+\sigma^{(2)}_{B}\mathbf{q}^{(2)}_{B}\mathbf{q}^{(2)\top}_{B}\mathbf{n}=\sigma^{(1)}_{C}\mathbf{q}^{(1)}_{C}\mathbf{q}^{(1)\top}_{C}\mathbf{n}+\sigma^{(2)}_{C}\mathbf{q}^{(2)}_{C}\mathbf{q}^{(2)\top}_{C}\mathbf{n}.$
(60)

This gives

$\sigma_{C}^{(2)}=\frac{t}{\cos^{2}(\beta)-\sin(\beta)\cos(\beta)}$ (61)

and

$\sigma_{B}^{(1)}=t\frac{\sin^{2}(\beta)-\sin(\beta)\cos(\beta)}{\cos^{2}(\beta)-\sin(\beta)\cos(\beta)}.$
(62)

We then have the stress state in the shear band as a function $t$ and $\beta$.
One can repeat the same exercise for regions D and E, or use symmetry
arguments to obtain the stress states. From this we obtain

$\sigma_{D}^{(2)}=\sigma_{B}^{(1)},\qquad\sigma_{D}^{(1)}=\sigma_{E}^{(1)}=\sigma_{A}^{(2)}=\sigma_{B}^{(2)},\qquad\sigma_{E}^{(2)}=\sigma_{A}^{(1)}=0.$
(63)

We determine that the traction applied on the bottom of region B is
$-\sigma_{B}^{(1)}\hat{\mathbf{e}}_{y}$, the traction at the top of region D
is $\sigma_{B}^{(1)}\hat{\mathbf{e}}_{y}$, and the traction on the right of
region E is $t\hat{\mathbf{e}}_{x}$. One can check that this ensures the net
moment applied to the sample is zero.

##  4 Finite volume formulation

Here, we develop a FV discretization that (i) is able to be applied to an
irregular geometry with an irregular grid; (ii) allows the direction of
information propagation to be in any direction in any dimension; (iii) is
stable and robust; (iv) can be applied with general boundary conditions; and
(v) handles discontinuities in the principal directions of stress.

We present the FV discretization in two dimensions because this is of
immediate practical importance (due to the prevalence of DIC data). However,
the generalization to three dimensions is straightforward.

###  4.1 Discretization

First, we consider the conservation law that applies to each volume
element222We use the term volume generically, and in this case, the volume
will be two dimensional. $\Omega$. We can simply integrate
$\mathrm{div}(\mathbf{T})=0$ over the volume then apply the divergence theorem
to obtain

$\begin{gathered}\int_{\partial\Omega}\mathbf{T}\mathbf{n}dS=0,\end{gathered}$
(64)

where $\partial\Omega$ is the boundary and $dS$ is an infinitesimal boundary
element. We specify the geometry in the deformed frame and use four sided
polyhedral elements with the positions of the nodes specified at each corner
(Fig. 7). We assume a constant traction $\mathbf{t}^{(j)}$ on each of the four
boundaries of the element. Hence, Eq. 64 becomes:

$\sum_{j,j=1}^{4}\mathbf{t}^{(j)}s_{j}=0,$ (65)

where $s_{j}$ is the length of the $j$th side of the polyhedral element. Here
we have assumed a fully two-dimensional deformation where the thickness does
not change, although there are straightforward modifications that can be made
to account for varying thickness.

Figure 7: The FV discretization. The volumes are irregularly shaped with four
nodes each. The force on each boundary must be computed.

The variable coefficient $c$ is specified to be constant inside each volume
element and to be the average value of two elements on the element boundary.
Note that we use the term element boundary to refer to the boundary between
two FV elements and domain boundary to refer to the edge of the entire domain.
Each element boundary will have two characteristics determined by $c$ on the
boundary, denoted $\\{\mathbf{q}^{(\alpha)}_{B}\\}$. Each element boundary
will also have an associated normal vector $\mathbf{n}$. We denote the cell in
the direction of $\mathbf{n}$ with a subscript $D$ and the cell in the
direction of $-\mathbf{n}$ with a subscript $C$ (Fig. 7). One must assign the
signs of $\mathbf{q}^{(\alpha)}$ such that they approximate a continuously
varying orientation of each characteristic, i.e. there are no abrupt reversals
in direction. One can arbitrarily choose which direction information should
propagate along for each characteristic, depending on which boundary
conditions one wishes to specify.

To develop a FV method, it is necessary to define a flux function on the
element boundaries. Note that the flux is simply the traction on the boundary
(which can be interpreted as momentum flux). Here we use a standard first
order upwinding scheme with characteristic decomposition. This ensures the
Courant–Friedrichs–Lewy condition is met, the consistency condition is met,
and that the method is stable. $\sigma^{(\alpha)}_{B}$ is specified using

$\sigma^{(\alpha)}_{B}=\left.\begin{cases}\sigma^{(\alpha)}_{C},&\text{if
}\mathbf{q}^{(\alpha)\top}_{B}\mathbf{n}>0,\\\ \sigma^{(\alpha)}_{D},&\text{if
}\mathbf{q}^{(\alpha)\top}_{B}\mathbf{n}<0.\\\ \end{cases}\right.$ (66)

One can express the stress at the boundary utilizing the spectral
decomposition. Hence, the traction along the $j$th boundary is

$\mathbf{t}^{(j)}=\sigma_{B}^{(1)}\mathbf{q}_{B}^{(1)}\mathbf{q}_{B}^{(1)\top}\mathbf{n}+\sigma_{B}^{(2)}\mathbf{q}_{B}^{(2)}\mathbf{q}_{B}^{(2)\top}\mathbf{n}.$
(67)

In the case where $\mathbf{q}^{(\alpha)\top}_{B}\mathbf{n}=0$ for some
$\alpha$, $\sigma^{(\alpha)}_{B}$ does not need to be computed as the
contribution to $\mathbf{t}$ will be zero.

The equations developed in Section 2.4 are utilized to specify the boundary
conditions. We distinguish two traction vectors at each point on the domain
boundary: the traction $\mathbf{t}^{*}$ that is given and that determines
$\sigma^{(\alpha)}$ on characteristics entering the domain, and the traction
$\mathbf{t}^{(j)}$ that arises due to the combination of entering and leaving
characteristics. For a real problem, there will be slight differences between
the two tractions due to error. When both characteristics are entering the
domain $\mathbf{t}^{(j)}=\mathbf{t}^{*}$. When both characteristics are
leaving, the traction is determined as part of the solution and
$\mathbf{t}^{*}$ does not play a role. For one characteristic entering and one
leaving, we can specify no domain boundary reflection and utilize Eq. 49.
Hence we have

$\mathbf{t}^{(j)}=\left.\begin{cases}\mathbf{t}^{*}&\text{if
}\mathbf{q}^{(1)\top}\mathbf{n}<0\text{ and
}\mathbf{q}^{(2)\top}\mathbf{n}<0,\\\
\sigma^{(1)}\mathbf{q}^{(1)}\mathbf{q}^{(1)\top}\mathbf{n}+\sigma^{(2)}\mathbf{q}^{(2)}\mathbf{q}^{(2)\top}\mathbf{n}&\text{if
}\mathbf{q}^{(1)\top}\mathbf{n}>0\text{ and
}\mathbf{q}^{(2)\top}\mathbf{n}>0,\\\
\mathbf{q}^{(1)}\mathbf{q}^{(1)\top}\mathbf{t}^{*}+\sigma^{(2)}\mathbf{q}^{(1)}\mathbf{q}^{(2)\top}\mathbf{n}&\text{if
}\mathbf{q}^{(1)\top}\mathbf{n}<0\text{ and
}\mathbf{q}^{(2)\top}\mathbf{n}>0,\\\
\mathbf{q}^{(2)}\mathbf{q}^{(2)\top}\mathbf{t}^{*}+\sigma^{(1)}\mathbf{q}^{(1)}\mathbf{q}^{(1)\top}\mathbf{n}&\text{if
}\mathbf{q}^{(1)\top}\mathbf{n}>0\text{ and
}\mathbf{q}^{(2)\top}\mathbf{n}<0.\\\ \end{cases}\right.$ (68)

where $\\{\mathbf{q}^{(\alpha)}\\}$ refer to the principal directions of the
stress in the element adjacent to the domain boundary.

For the cases considered by the authors, it was not necessary to specify the
case where $\mathbf{q}^{(\alpha)\top}\mathbf{n}=0$ on the domain boundary
because these vectors are numerically specified using floating point numbers,
i.e. the dot product $\mathbf{q}^{(i)\top}\mathbf{n}$ will rarely be exactly
zero. However, if such a case did arise one should specify that the traction
component parallel to the boundary should be zero.

Eqs. 67 and 68 are sufficient to determine the $\mathbf{t}^{(j)}$, which are
then substituted into Eq. 65. This provides two linear equations for each
finite volume. Each finite volume also has two unknown principal stresses. The
equations and unknowns can be assembled into a linear system of 2$n$ equations
and 2$n$ unknowns, $\mathbf{A}\mathbf{x}=\mathbf{b}$, where $n$ is the number
of finite volumes. This can then be solved using standard linear algebra
algorithms. Once $\\{\sigma^{(\alpha)}\\}$ has been determined for each
element, $\mathbf{T}$ can be determined using the known principal directions.
A python implementation of this code, showing how the matrix is assembled and
solved, is provided in the supplementary material.

###  4.2 Validation

The developed FV discretization is designed for problems with varying material
properties, discontinuous strain fields, and irregular geometries. We first
consider a problem developed to assess the performance of the discretization
in the case of a discontinuity. Specifically, we consider an infintesimal
elastic deformation with discontinuous material properties. In Section 4.2.2
we briefly consider the algorithm applied to two problems from other articles
to assess performance in the case of continuously varying material properties
and an irregular computational mesh. The following procedure is used in each
case: (i) stress and displacement/strain field is computed using a FE
simulation implemented with a commercial FE solver, (ii) the kinematic field
and boundary conditions are exported, (iii) the strain field is used to
compute the principal directions of stress, (iv) the stress field is computed
as a function of position using the developed FV discretization, (v) the
stress field computed principal directions is compared to that output by the
commercial FE solver and the error is quantified. In this last step we use the
normalized mean absolute error (NMAE):

$\mathrm{NMAE}(\Delta\mathbf{T})=\frac{\mathrm{mean}(|\Delta\mathbf{T}|)}{\mathrm{mean}(|\mathbf{T}|)},$
(69)

where $\Delta\mathbf{T}$ is the difference between the stresses in a
particular element, $\mathbf{T}$ is the stress computed from the FE
simulation, and the mean is taken spatially over all elements in the
computation.

####  4.2.1 Strain discontinuity

Consider the two-dimensional infintesimal-strain linear-elastic problem
depicted in Fig. 8. The domain is square with two regions, one with an elastic
modulus $C_{1}=50$ GPa, and another with an elastic modulus of $C_{2}=200$
GPa. The Poisson’s ratio $\nu$ is $0.3$ for both regions. A constant tensile
traction of 100 MPa is applied on the right and left boundaries while the
upper and lower boundaries are traction free.

Figure 8: The problem investigated in order to validate the algorithm in the
case of a discontinuity. There are two regions of the material with different
elastic moduli $C$. A horizontal traction is applied on the left and right
boundaries.

Abaqus standard (2020) is used to conduct the FE simulations. The mesh is a
regular array of square CPS4R shell elements. Simulations are conducted at a
variety of resolutions, the finest being 4000$\times$4000 elements. The
element boundaries coincide with the boundaries between the two domains of
differing elastic moduli. Fig. 9 shows the output strain field.

The principal directions are determined at each point by computing the
eigenvectors of the strain output from the FE simulation. In order to assign
the direction of information propagation for each volume, the following
procedure is used: (i) the direction of information propagation is chosen for
one ”seed” volume, (ii) the direction of information propagation is chosen to
be as similar as possible in adjacent cells (the dot product of vectors
corresponding to the direction of information propagation is set as close as
possible to 1), (iii) the prior step is repeated until the direction of
information propagation is set over the entire domain. For this problem, in
addition to the discontinuity present in the strain field, there is a
discontinuity in the orientation of the principal directions allowing us to
assess the algorithm in this case (Fig. 10).

Figure 9: The strain field output from the FE simulation and input into the FV
simulation. Figure 10: The orientation of the principal directions of the
stress and strain as a function of position. In particular, this shows the
jump between the regions with two different elastic moduli.

The same computational mesh is used for the FV procedure as the FE procedure
in order to avoid unnecessary interpolation issues. A linear system of
equations is assembled and solved using a MATLAB’s direct sparse matrix
solver. At the finest resolution, solving the sparse system of equations took
663 s on an Intel Xeon 6226 2.70GHz CPU. Note that there were memory
requirements on the order of 100 GB. Once the principal stresses are
determined, the full stress tensor is computed at every point using the
principal directions.

The stress field obtained is compared to that obtained using the FE method in
Fig. 11, with the error plotted as a function of position in Fig. 12. There
was an NMAE of 8.1$\times 10^{-4}$ at the finest resolution. The solution
steadily converged to the true solution as the number of volumes $n$ was
increased (Fig. 13).

Figure 11: A comparison of the two stress fields: a field obtained via a FE
simulation with specified material properties and the field obtained via FV
simulation using the strain field but with no known material properties.
Figure 12: The normalized absolute error $|\Delta\mathbf{T}|/|\mathbf{T}|$
where $\Delta\mathbf{T}$ is the difference between the stress computed with
the FE simulation and the stress computed using the FV method and the
denominator is the stress computed via the FE method at each point. Figure 13:
The normalized mean absolute error (NMAE) as a function of the number of
volumes used in the FV simulation. This shows the output is converging to the
true solution as the computational mesh is refined. Figure 14: The stress
field computed via the FE simulation compared with the FV simulation with
different numbers of volumes $n$. (a) A comparison with no discontinuity in
the strain with $x=1.5=3l_{x}/4$. (b) A comparison with a discontinuity in the
strain and principal directions ($\Delta\theta=8.9^{o}$). The FE curve is
taken with $n=1.6\times 10^{7}$, however, there is no visible discrepancy
between this curves and curves generated with a much smaller number of
elements.

The error is primarily localized at the discontinuity (Figs. 12 and 14) due to
issues interpolating $\mathbf{q}$ at the boundary. As the principal directions
are approximately aligned with the discontinuity in some regions, small errors
accumulate as the stress propagates along characteristics adjacent to the
boundary. One can see some streaks in Fig. 12 which corresponds to erroneous
principal stresses propagating away from the discontinuity along
characteristic lines. Nevertheless, we emphasize that these errors diminish as
the computational mesh is refined, and the current FV discretization appears
to converge to the true solution (Fig. 13). For example, Fig. 14 shows
significant error close to the discontinuity at $n=4\times 10^{4}$, but a much
smaller error at $n=1.6\times 10^{7}$.

One limitation of the developed FV procedure in the case of a discontinuity is
that a particularly fine computational mesh is required: orders of magnitude
more elements/volumes than a FE simulation for comparable accuracy. In
addition, issues such as numerical diffusion may arise in other cases, as this
depends on the orientation of the characteristic lines with respect to the
element boundaries. The developed method is only first order accurate. This
may be acceptable for many cases, however, substantial improvements to
computational time, memory requirements and accuracy would be possible with a
higher order discretization within the current FV framework.

####  4.2.2 Other validation examples

We next consider the validation example described in detail in Cameron and
Tasan [17] and briefly summarized here. We have a two-dimensional rectangular
sample with continuously varying material properties; specifically, yield
stress and hardening coefficient. Technically, this should be interpreted as
non-linear elasticity simulation because only one strain increment is
simulated (see paper for details). A constant traction is applied at both ends
of the sample so that the entire specimen is in tension. The simulation is
conducted using 400$\times$100 square elements. A much lower spatial
resolution is required compared to Section 4.2.1. For this problem, the NMAE
was 9.51$\times 10^{-4}$. This was higher than the discretization presented in
Cameron and Tasan [17] which has a NMAE of 2.52$\times 10^{-4}$. Although that
method is more accurate for this example, it is generally unstable and cannot
be widely applied as frequently will violate the Courant–Friedrichs–Lewy
condition. We also not that the specific error is of secondary importance
compared to stability because the FV discretization can reach arbitrary
accuracies through refinement of the numerical grid (see second validation
example and Cameron and Tasan [19]). The results are shown along two line
profiles in Fig. 15.

Figure 15: Comparison of the stress computed from the FE simulation described
in Cameron and Tasan [17] (solving the forward problem), the stress computed
using the finite difference method presented in that same paper (solving the
inverse problem), and the stress computed using the FV method presented in
this paper (solving the inverse problem). (a) $T_{xx}$ plotted as a function
of $x$ with $y=l_{y}/2$. (b) $T_{xy}$ plotted as a function of $y$ with
$x=3l_{x}/4$.

We consider a further example obtained from Cameron and Tasan (submitted) to
validate the method for an irregular geometry with an irregular mesh. This
example is presented comprehensively in that article and is only briefly
summarized here. There is a a two-dimensional elasto-plastic finite-strain
temporal simulation of the diffuse necking in a flat tensile specimen. The
simulation domain is a flat rectangular geometry with a length to width ratio
of 4.5 and comprised of 7200 shell elements. The local accumulated plastic
strain reaches values of 89.9 %, and the square undeformed elements become
significantly deformed, allowing the method to be tested for irregular
geometries. The FV discussed in this paper was implemented at each time step,
however, $\mathbf{Q}$ is determined using a more complex method discussed in
Cameron and Tasan (submitted) that incorporates elasto-plasticity. The max
NMAE error (computed at one increment in time but averaged over all space) is
2.1$\times 10^{-3}$, however, this is attributed to uncertainty in
$\mathbf{Q}$ rather than the FV discretization because of the sensitivity of
the error to the time step used.

Note that a substantial portion of the error arises from elements directly
adjacent to the boundary in all validation cases. For example, in the second
validation example, the NMAE is 4.39$\times 10^{-4}$ when elements on the
boundary are excluded (approximately half of the value computed when the
boundaries are included). This issue remains when alternate $\mathbf{d}^{*}$
vectors are used to specify the traction, or when using the finite difference
discretization approach presented in Cameron and Tasan [17]. Furthermore, this
issue has been observed in all problems considered by the authors. The
phenomenon likely arises because the characteristics at the boundary should be
theoretically aligned with the boundary in many cases, whereas the
characteristics of the volume element adjacent to the boundary are not. We
observe that errors accumulate along the characteristic that is approximately
aligned with the boundary (e.g. see Fig. 16). Fortunately, from a practical
perspective, one can simply ignore the information obtained close to the
boundary.

Figure 16: The error in the stress computed for the first validation example
plotted in a small sub-region of the domain. This is used to show the error
accumulation on the boundary. The left and upper boundaries correspond to
domain boundary. The directions of information propagation are in the downward
and rightward directions.

##  5 Conclusion

Analysis of the governing equations was conducted prior to developing the
finite volume discretization. Key results from this analysis include:

  * 1.

It is mathematically shown that the characteristic directions are aligned with
the principal directions of the stress. This leads to the interpretation of
the propagation of principal stress components throughout the material.

  * 2.

It is shown how a principal direction discontinuity gives rise to reflection
and transmission phenomena. This is mathematically similar to the classic
second order wave equation in variable wave-speed media.

  * 3.

The correspondence between known traction at the boundary and the specific
boundary conditions is elucidated for different configurations of
characteristics entering and leaving the domain.

  * 4.

The solution to the Riemann problem and an idealized shear band is given.
These examples give the reader intuition for how reflection and transmission
of stress apply to a physical problem.

This analysis enabled the development of a robust finite volume
discretization. This is necessary for the broader applicability of the stress
computation approach presented in Cameron and Tasan [17], Liu [18], Cameron
and Tasan (submitted) and the present paper. Unlike the method used in Cameron
and Tasan [17], this method is stable and works for general geometries, large
deformation, and more complex characteristic structures where information may
propagate in any direction. Unlike the method proposed by Liu [18], this
method generalizes to any boundary conditions and is robust to
discontinuities. The finite volume discretization has been validated on three
examples: one infinitesimal strain problem with discontinuous material
properties and strain field, one infinitesimal strain problem with material
properties and strain field, and one finite strain problem with with
homogeneous material properties. There are two limitations of this analysis
and discretization: (i) it does not apply to deformations with repeated stress
eigenvalues, (ii) it does not apply in cases where the characteristic curves
form closed loops. These limitations will likely be addressed in future work.

## References

  * Pan et al. [2009] B. Pan, K. Qian, H. Xie, A. Asundi,  Two-dimensional digital image correlation for in-plane displacement and strain measurement: A review,  Measurement Science and Technology 20 (2009). doi:10.1088/0957-0233/20/6/062001.
  * Hild and Roux [2006] F. Hild, S. Roux,  Digital image correlation: From displacement measurement to identification of elastic properties - A review,  Strain 42 (2006) 69–80. doi:10.1111/j.1475-1305.2006.00258.x.
  * Blaber et al. [2015] J. Blaber, B. Adair, A. Antoniou,  Ncorr: Open-Source 2D Digital Image Correlation Matlab Software,  Experimental Mechanics 55 (2015) 1105–1122. doi:10.1007/s11340-015-0009-1.
  * Réthoré [2010] J. Réthoré,  A fully integrated noise robust strategy for the identification of constitutive laws from digital images,  International Journal for Numerical Methods in Engineering (2010) 1885–1891. doi:10.1002/nme.
  * Siddiqui et al. [2017] M. Z. Siddiqui, S. Z. Khan, M. A. Khan, M. Shahzad, K. A. Khan, S. Nisar, D. Noman,  A Projected Finite Element Update Method for Inverse Identification of Material Constitutive Parameters in transversely Isotropic Laminates,  Experimental Mechanics 57 (2017) 755–772. doi:10.1007/s11340-017-0269-z.
  * Viala et al. [2018] R. Viala, V. Placet, S. 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 (2018) 91–103. doi:10.1016/j.compositesa.2017.12.018.
  * Crouzeix et al. [2009] L. Crouzeix, J. N. Périé, F. Collombet, B. 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 (2009) 1732–1740. doi:10.1016/j.compositesa.2009.08.016.
  * Pierron and Grédiac [2012] F. Pierron, M. Grédiac, The Virtual Fields Method: Extracting Constitutive Mechanical Parameters from Full-Field Deformation Measurements, Springer Science & Business Media, 2012\.
  * Florentin and Lubineau [2010] E. Florentin, G. Lubineau,  Identification of the parameters of an elastic material model using the constitutive equation gap method,  Computational Mechanics 46 (2010) 521–531. doi:10.1007/s00466-010-0496-y.
  * Florentin and Lubineau [2011] E. Florentin, G. 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 (2011) 227–234. doi:10.1007/s12008-011-0129-5.
  * Leygue et al. [2018] A. Leygue, M. Coret, J. Réthoré, L. Stainier, E. Verron,  Data-based derivation of material response,  Computer Methods in Applied Mechanics and Engineering 331 (2018) 184–196.
  * Stainier et al. [2019] L. Stainier, A. Leygue, M. Ortiz,  Model-free data-driven methods in mechanics: Material data identification and solvers,  Computational Mechanics 64 (2019) 381–393. doi:10.1007/s00466-019-01731-1.
  * Langlois et al. [2022] R. Langlois, M. Coret, J. Réthoré,  Non‐parametric stress field estimation for history‐dependent materials: Application to ductile material exhibiting Piobert–Lüders localization bands,  Strain 58 (2022). URL: https://onlinelibrary.wiley.com/doi/10.1111/str.12410. doi:10.1111/str.12410.
  * Fletcher et al. [2021] L. Fletcher, F. Davis, S. Dreuilhe, A. Marek, F. Pierron,  High strain rate elasto-plasticity identification using the image-based inertial impact (IBII) test part 1: Error quantification,  Strain 57 (2021). doi:10.1111/str.12375.
  * Rahmani et al. [2014] B. Rahmani, I. Villemure, M. Levesque,  Regularized virtual fields method for mechanical properties identification of composite materials,  Computer Methods in Applied Mechanics and Engineering 278 (2014) 543–566. doi:10.1016/j.cma.2014.05.010.
  * Pierron and Grédiac [2021] F. Pierron, M. Grédiac,  Towards Material Testing 2.0. A review of test design for identification of constitutive parameters from full-field measurements,  Strain 57 (2021). doi:10.1111/str.12370.
  * Cameron and Tasan [2021] B. C. Cameron, C. Tasan,  Full-field stress computation from measured deformation fields: A hyperbolic formulation,  Journal of the Mechanics and Physics of Solids 147 (2021) 104186. doi:10.1016/j.jmps.2020.104186.
  * Liu [2021] C. Liu,  Nonuniform Stress Field Determination Based on Deformation Measurement,  Journal of Applied Mechanics 88 (2021) 071005. doi:10.1115/1.4050535.
  * Cameron and Tasan [2022] B. C. Cameron, C. C. Tasan, Partial differential equations to determine elasto-plastic stress-strain behavior from measured kinematic fields (in submission), 2022\.
  * Strang [2007] G. Strang, Computational Science and Engineering, Wellesley-Cambridge Press, Wellesley, MA, 02382 USA, 2007.
  * Alford et al. [1974] R. M. Alford, K. R. Kelly, D. M. Boore,  Accuracy of finite-difference modeling of the acoustic wave equation,  Geophysics 39 (1974) 834–842.
  * Molenkamp [1968] C. R. Molenkamp,  Accuracy of finite-difference methods applied to the advection equation,  Journal of Applied Meteorology and Climatology 7 (1968) 160–167.
  * Karal and Keller [1959] F. C. Karal, Jr., J. B. Keller,  Elastic wave propagation in homogeneous and inhomogeneous media,  The Journal of the acoustical society of america 31 (1959) 694–705.
  * Wang [2002] Z. J. Wang,  Spectral (finite) volume method for conservation laws on unstructured grids. basic formulation: Basic formulation,  Journal of computational physics 178 (2002) 210–251.
  * Chai et al. [1994] J. C. Chai, H. S. Lee, S. V. Patankar,  Finite volume method for radiation heat transfer,  Journal of thermophysics and heat transfer 8 (1994) 419–425.
  * Jasak [1996] H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows., Ph.D. thesis, Imperial College London, 1996.
  * Mingham and Causon [1998] C. G. Mingham, D. M. Causon,  High-resolution finite-volume method for shallow water flows,  Journal of Hydraulic Engineering 124 (1998) 605–614.
  * Hermeline [2000] F. Hermeline,  A finite volume method for the approximation of diffusion operators on distorted meshes,  Journal of computational Physics 160 (2000) 481–499.
  * Mooney [1940] M. Mooney,  A theory of large elastic deformation,  Journal of Applied Physics 11 (1940) 582–592. doi:10.1063/1.1712836.
  * LeVeque [2002] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, volume 31, Cambridge university press, 2002.