{'fraction_non_alphanumeric': 0.03712884922602325, 'fraction_numerical': 0.009621917405754788, 'mean_word_length': 4.5339904481251905, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deterministic calculation of elasto-plastic stress-strain behavior from arbitrary deformation fields

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.


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 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: where b is the body force,ẍ is the acceleration, and ρ is the density. Here, Q is the matrix of stress eigenvectors, i.e. the spectral decomposition of T is T = QΛ T Q T 1 . Eq. (1a) is simply force balance and Eq. (1b) constrains T to have the principal directions given by Q. This can be intuitively understood as expressing T in the coordinate system aligned with its principal directions and then specifying the shear components in this coordinate system to be zero. Assuming Q,ẍ, ρ, and 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 Q can be determined at every point in space from the measured strain (or strain rate), these equations can be directly solved for T. It is then not necessary to apply standard simulation approaches where material properties must be specified.
It is straight forward to determine 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 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 1 Notation: Bold upper-case letters to denote second order tensors (A, B,...), bold lowercase letters to denote vectors (a, b,...), unbolded letters correspond to scalars (a, b, ...), and the superscript T denote the transpose (A T , a T , ...). We use standard matrix notation, for example, we have a · b = a T b and a ⊗ b = ab T . T denotes the Cauchy stress, b denotes the body force, x donates the position of a material element in the deformed frame, and ρ 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Ȧ andÄ respectively. DivA and divA correspond to the divergence A with respect to the initial and deformed reference frame respectively. |A| = AijAij gives the magnitude of A (Einstein summation notation is used). We denote the deviatoric part of a second order tensor A as A0, and the trace as tr(A). The spectral decomposition of a symmetric tensor A is written as A = Q A Λ A Q T A where Q A has columns with the eigenvectors of A denoted {q (i) |i = 1, 2, 3} (principal directions) and Λ A is a diagonal matrix containing the eigenvalues {λ i |i = 1, 2, 3} (principal values). This decomposition is also expressed as As the spectral decomposition is used extensively for the stress, we use the special notation: T = i σ (i) q (i) q (i)T where {σ (i) |i = 1, 2, 3} are the principal stresses and {q (i) |i = 1, 2, 3}. 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 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 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].

Preliminaries
We use standard notation to specify the kinematics of large deformation. We have a smooth one to one mapping from the reference body X to the deformed body is given by x =x(X, t). This can be used to specify the deformation gradient F = Grad(x) (2) and the velocity gradient L =ḞF −1 .
(3) This can be decomposed into symmetric and asymmetric parts The left Cauchy-Green tensor is given by For elasto-plastic deformation we take the standard Kröner-Lee multiplicative decomposition Similarly we have L e =Ḟ e F e−1 , In standard theories of isotropic elasto-plasticity we have that T has the same eigenvectors as the elastic left Cauchy deformation tensor B e = F e F eT . We also have that the stress (Cauchy stress) has the same principal directions as F e D p F eT .
The key challenge in determining the principal directions of T in the case of elasto-plastic deformation is that one can only directly obtain B and D from the observed deformation. B e , F e and D p are unknowns. This is not a problem in the simplified case of elastic deformation as one can simply take B e = B. Similarly in the case of plastic deformation where the elastic strain is negligible one can take F e D p F eT = D. Note that this translates to infinitesimal deformation if one replaces (B, D, B e , F e D p F eT ) with the infinitesimal deformation measures (E,Ė, E e ,Ė p ).

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 B e =B e (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 elastoplastic 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.

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 = τ to t = τ + ∆t. It is assumed that the Cauchy stress T is known at time τ and we wish to determine the stress at time τ + ∆t based on the observed deformation ( Fig. 1). We introduce subscripts to denote quantities evaluated at these times, e.g. T τ =T(τ ) and T τ +∆t =T(τ + ∆t). This increment of deformation is characterized by the deformation gradient F * : We define F p * in the similar way: This is the relationship depicted in the triangle in Fig. 1. We express the spectral decomposition of T as where {q (i) |i ∈ 1, 2, 3} are the eigenvectors and {σ (i) |i ∈ 1, 2, 3} are the eigenvalues such that: (observed) (known from , up to rotation) principal direction of undeformed reference frame Figure 1: The references frames and deformation gradients used to map between them. The sample is initially deformed with F τ = F e τ F p τ , then an additional deformation F * is observed which is related to the unobserved F p * and F e τ +∆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.
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 T (see Section 1), we have: Here, we use the convention that {q (i) } always corresponds to the eigenvector of the stress, however, {λ (i) } are just generic eigenvalues. If we wish to refer to the eigenvectors and eigenvalues of some other specific tensor A we will add a subscript so we have {q We also take the eigenvectors to be normalized. Hence, we have q (i)T q (j) = δ ij where δ 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 F e D p F eT we have where the {λ (i) } are different from Eq. (16) but {q (i) } are still the eigenvectors of the stress. Consider the following expression: This will have the same eigenvectors as T τ +∆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: We manipulate this expression in order to find an expression for the eigenvectors in terms of observed quantities. We have: where, I + ∆tL p τ +∆t + ∆tL pT τ +∆t = (I + ∆tL p τ +∆t )(I + ∆tL p τ +∆t ) T − ∆t 2 L p τ +∆t L pT τ +∆t This last expression follows from Substituting Eq. (22) into Eq. (21) gives: From Eq. (12) we have F e τ +∆t F p * = F * F e τ , hence: This shows that when ∆t is small we can determine the eigenvectors for T τ +∆t from the observed F * and the quantity F e τ F eT τ which known from T τ and the functionB e (T). We can use this expression to derive a differential equation for the evolution of {q (i) }.

For some arbitrary tensor A with eigenvectors {q
we have the following general expression for the derivative of its eigenvectors: where this assumes there are no repeated eigenvalues and ensures that {q This will have the same eigenvectors as the stress up to error of O(∆t 2 ) (Eq. (25)). Hence, in the limit of ∆t → 0, we have q A =q (i) . We will determine an expression forȦ and {λ (i) A } to derive an expression foṙ q (i) . Technically, A will be a function τ and τ + ∆t as F e τ is evaluated at τ and F * depends on both τ and τ + ∆t. However, because A is used to determine {q (i) } at times just after τ , we consider τ fixed and allow t = τ + ∆t to vary. Hence we have A =Ā(t) andȦ = dA/dt. We then havė where this substitution arises from a similar argument to that made in Eq. (23). Substituting for B e and dropping the subscripts τ we have:Ȧ Evaluating {λ B } are the eigenvalues of B e . Substituting these quantities into Eq. (26) gives:q When combined with force balance, the assumption that B e can be determined by T, and that 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 D and spin W. Using Eq. (16) we have: Substituting back into Eq. (32) gives: Note that for i = j, q (j)T Wq (i) = 0 as W is antisymmetric. Hence, we can replace the second sum in Eq. (34) with a sum over all j. Then the term j q (j) q (j)T becomes the identity I. This gives: To summarize these results, we now have a complete set of partial differential equations 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

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, D ≈Ė. We also have that λ (i)e B will be close to 1, we express this as λ e(i) We neglect the O( ) term in the numerator as it is small compared to 2. We note that λ e(i) As the denominator of this sum expression will be O( ), the second Wq (i) term will be negligible in comparison. This gives the following small strain formulation: 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 T =T(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 B e = B, Eq. (32) becomes: Note the following relationshipḂ We now haveq This is simply the same as the expression for the rate of change of the eigenvector of B (see Eq. (26)). Hence, choosing q to be an eigenvector of B will solve the equation and we can replace Eq. (35) with 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 F e → I indicating no elastic strain. Here all λ e Bi → 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 Since the q (i) vectors are orthogonal, i.e. q (i)T q (j) = 0 for all i = j, one can see that q (i) are the eigenvectors of 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.

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 {σ (k) } denote the set of principal stresses that are repeated over some region of (x, t). The corresponding principal directions, {q (k) }, are not unique and any linear combination of {q (k) } is also a principal direction. Note that this also corresponds to the repeated principal values of B e , {λ e(k) B }, as the functionB e (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 {q (i) } will be non-unique. To resolve this ambiguity we additionally require that each q (i) is continuously varying and differentiable with respect to time. This restriction gives rise to unique {q [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 {q However, as the material is unstressed at time t = 0, we have that B t0 = I. Therefore {q So {q (i) t0 } will be consistent with L e t0 + L eT t0 .
Second, we have that {q (i) t0 } will be consistent with F e t0 D p t0 F eT t0 . However, as F e t0 = I, the stress will be aligned with D p . Due to the definition of D p , {q (i) t0 } will also be consistent with L p t0 + L pT t0 .
As {q t0 } is consistent with both L e t0 +L eT t0 and L p t0 +L pT t0 , we have that {q given: Hence we have D t0 q Therefore, the initial conditions for {q (i) } can be obtained from the observed quantity 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 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 t, one can use the equations explained in [24] to specify these components.
where n is the domain boundary normal, t is the known traction, and {d (j) * } is a subset of size n e of an arbitrarily chosen orthonormal basis {d (i) }. These n e equations can be used to solve for the unknown principal stresses at the boundary.

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ẍ = 0 and b = 0. Where needed, incorporation of these terms into the linearization can be carried out without significant challenges.
We introduce a small quantity , where << 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 In general we treat {σ (i) } and {q (i) } as our independent variables, though we will often refer to T and T [0] where convenient with Furthermore, we simply refer toḞ, {q (i) } and L with no subscripts as second derivatives will be small in comparison to other terms. We have in where div [0] is the divergence with respect Next we consider the linearized form of Eq. (36b). Substituting Eq. (32) into Eq. (50b) gives: 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: We also define F * = ∂x/∂x [0] . We have: The referential form of the force balance equation gives Consider the summation component of this equation: We note that q (i) [1] will be O( ) since Eq. (56) requiresq (i) to be perpendicular to q (i) (with O( ) error). Hence, this simplifies to Hence, Eq. (61) becomes: Then substituting for t [1] L using Eq. (58b) we have: Using Eq. (53) gives: This gives the linearized form of force balance.
We denote t * as the first Piola traction vector with respect to the x [0] reference frame such that: We now consider the limit as → 0 to obtain the linearized system: in 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

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 (x =x(X, t) is sufficient as one can compute 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 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 F, which we denote F , one must (i) understand the error introduced in q and (ii) understand how the error in 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 {q (i) } are computed directly from B e + B e , where B e is the error. Assuming B e << B e , we have that the error in {q (i) } will scale with |B e 0 | due to the theory of matrix perturbations (the subscript 0 refers to the deviatoric component of the B e ). The error will be small if B e 0 << B e 0 . Second, we consider a deformation step in the plastic limit, where the principal directions are computed directly from D p + D p where D is the error. Again, assuming D p << D p , we have that the error in {q (i) } will scale with |D p 0 |. The error will be small if D 0 << 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: The accuracy of typical DIC strain data typically ranges between 1 × 10 −3 and 1 × 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.

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).

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 .
forward problem, undeformed forward problem, deformed plastic region inverese problem, undeformed inverse problem, deformed 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.
Here, ξ(t) is is a cubic function specified so that dξ(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 ν = 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, ε t = 2/3|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.

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 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 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). 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 =h(x, t). (72) As the elastic volume change will be negligible for this problem, we assume det(F ij ) = 1. This gives h/h 0 = 1/det(F αβ ) at every point in the sample throughout deformation where h 0 is the initial thickness. Note we use Greek letters (α, β, ...) 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, det(F αβ ) 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 T where T = hT: 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. T j corresponds to the stress at the j increment. Instead of computing the rate of change of the q values directly, we compute A using Eq. (27). We obtain the principal directions at the time step j from A: where F * = F j F −1 j−1 . This is done for each material point. For this problem, 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. F * , B and the thickness h can then be computed from F. Once the principal directions {q (i) j } are obtained from the computed 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: Algorithm 1: Compute T as a function of position and time j ← 0; compute T j from {q j } and boundary conditions for all X using the finite volume method; compute T j from T j using Eq. (73a) and F j for all X using Eq. (72); compute B e j from T j for all X using Eq. (37); compute A j+1 from B e j and F for all X using Eq. (74); compute {q 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ξ(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 ).

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.
We quantify the error as a function of position using |∆T|/mean x (T), where ∆T is the difference between the stress computed from the forward problem and the stress computed from the inverse problem, and the mean x is taken over all positions (but only at the time considered). The T used in the denominator is the 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 {q (i) } is changing more rapidly with time, and the finite time step will cause a larger error.
We also quantify the error as a function of time using the normalized mean absolute error (NMAE): 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. t = t f t = (11/12)t f 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 .
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.
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.
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 ∆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 (ε = 0.26). However, information at higher strains is not available. The approach presented in this paper gives accurate stress-strain behavior up to ε = 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 (a) t = (11/12)t f , n t =1200 (b) t = t f , n t =1200 (c) t = t f , n t = 12  t/t f n t = 1200 n t = 240 n t = 120 n t = 24 n t = 12 assumptions such as constant elastic modulus, and the lack of texturing of the sample. These results simply provide a proof of principal demonstration.

Conclusion
The method developed in [23], previously valid for isotropic elastic or plastic deformation, is generalized to elastoplastic 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 (a) t = (11/12)t f , n t =1200 (b) t = t f , n t =1200 (c) t = t f , n t = 24 specified relationship in forward problem from inverse problem (near x = 0 boundary) from inverse problem (location of max ) from conventional method 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 σ = 3/2|T 0 | and strain measure ε = ln(1 + 2/3|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. 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.