## Abstract

In this investigation, we demonstrate a technique that, under certain circumstances, will determine stresses associated with a nonuniform deformation field without knowing the detailed constitutive behavior of the deforming material. This technique is based on (1) a detailed deformation measurement of a domain (currently in 2D) and (2) the observation that for isotropic materials, the strain and the stress, which form the so-called work-conjugate pair, are co-axial, or their eigenvectors share the same directions. The particular measures for strain and stress chosen in this study are the Lagrangian (or Green-Lagrangian) strain and the second Piola–Kirchhoff stress. The deformation measurement provides the field of the principal stretch orientation θ_{λ} and since the Lagrangian strain and the second Piola–Kirchhoff stress are co-axial, the principal stress orientation θ_{s} of the second Piola–Kirchhoff stress is then determined. The Cauchy stress is related to the second Piola–Kirchhoff stress through the deformation gradient tensor, which can be measured experimentally. We then show that the principal stress orientation θ_{σ} of the Cauchy stress is the sum of the principal stretch orientation θ_{λ} and the local rigid-body rotation θ_{q}, which is determinable by the deformation gradient through polar decomposition. Such a relationship is valid for finite deformations. With the principal stress orientation θ_{σ} known, the equation of equilibrium, now in terms of the two principal stresses, σ_{1} and σ_{2}, and θ_{σ}, can be solved numerically with appropriate traction boundary conditions. The stresses determined using this technique obviously satisfy the equation of equilibrium, in contrast to those obtained from a constitutive model with input from deformation measurement. The technique and the associated numerical scheme are verified and validated through two virtual test cases representative of the simply-connected and multiply-connected domains, where exact solutions are available. The technique is then applied to an experimental case of nonuniform deformation of a polyvinyl chloride (PVC) sheet with a circular hole and subject to uniaxial tension. In this case, the associated stress field is also determined through a constitutive model of hyperelasticity, the generalized neo-Hookean (GNH) model, calibrated for the particular PVC sheet. Limitations and restrictions of the technique and the associated numerical scheme, as well as possible extensions will be discussed.

## 1 Introduction

Associated with a nonuniform deformation field is usually a nonuniform stress field, either as the driving force to produce the deformation or as the resistance of the material against externally imposed deformation. For characterizing and quantifying deformation and failure process of a structure or a material, both the deformation and the stress fields are desired. However, there is a tremendous disparity in our ability of experimentally measuring deformation and experimentally measuring stress. Advances in experimental techniques in the past several decades have enabled us to quantitatively measure the deformation field to a very high accuracy, especially in the case of two-dimensional (2D) deformation. On the other hand, not much can be said about the stresses associated with the nonuniform deformation field, even though every detail of the deformation can be obtained experimentally.

The usual scheme for inferring stresses associated with a deformation field relies on the constitutive model that connects stress to strain/deformation, and the constitutive relation itself might need to be measured and calibrated through separate independent experiment on some simple stress/deformation configurations, e.g., uniaxial tension or compression. The problem with this approach is that the accuracy and fidelity of a particular constitutive model, as well as the parameters associated with it, are often questionable. Moreover, the stress field determined in such way does not necessarily satisfy the equilibrium equation and the traction boundary conditions.

*σ*=

*F*/

*A*, where

*F*is the applied force and

*A*is the cross-sectional area of the cylinder. As all these requirements are satisfied, the experimenter declares or determines that the stress, as a tensor, within the gage section is given by

*x*

_{3}is along the axial direction of the cylinder. Reaching such a conclusion is based on two considerations. The first consideration is the belief that the stress associated with a homogeneous deformation is also homogeneous or uniform. The uniformity of the deformation has been checked and ensured during the uniaxial tension test. The second consideration, which is more important and convincing, is that the stress given in such manner satisfies the equation of equilibrium, as well as the traction boundary conditions, not only over the top and the bottom surfaces of the cylinder, but also over the lateral surface of the cylinder.

For any testing configuration that is slightly more complicated than uniaxial stress, measuring stresses becomes hopeless, even more so if the deformation field associated with the testing configuration is nonuniform. Nevertheless, if one accepts the fact that stresses cannot be directly measured experimentally, at least in the usual sense, but only be determined as suggested in the uniaxial tension case discussed above, we may seek to determine the stresses associated with a nonuniform deformation field, at least for some limited scopes.

Numerical simulation approaches, like the finite element analysis (FEA), have become a mature tool for assessing engineering design and structure/material performance [1–3]. The FEA method seeks numerical/approximate solution to the full set of partial differential equations (PDEs) involving kinematics, kinetics, and constitutive relations, together with the appropriate boundary and initial conditions, of the deforming domain. In principle, FEA does not require experimental input for completing the simulations, except some calibration tests for the determination of parameters in the constitutive model for specific materials. The predictions from a FEA simulation are often compared with the experimental measurement for validation of the simulation results, and the deformation field is used for such comparison.

In recent years, various techniques were developed, e.g., the virtual field technique [4–6], by taking advantage of the availability of the detailed deformation field, to quantitatively determine the material parameters in a constitutive model. By maintaining the deformation/strain field, obtained from experimental measurement, and by assuming certain form of the constitutive model, the material parameters in the constitutive model are iteratively updated such that the equation of equilibrium and the boundary conditions are satisfied. It is the objective of the present study to identify the requirement and the scope, where it is possible to obtain nonuniform stress field based on the deformation measurement without knowing the detailed constitutive behavior of the deforming body.

In the present investigation, we demonstrate a technique, together with the numerical scheme, for determining the stresses associated with a nonuniform deformation field. In addition to the detailed deformation measurement of a 2D domain, this technique is based on an observation that for *isotropic* materials, the strain and the stress, which form the so-called work-conjugate pair, are co-axial, or their eigenvectors share the same directions. Therefore, the second limitation of the technique is that we only deal with isotropic materials. Note that “isotropy” is indeed a constitutive description of a material, but the technique we are discussing in this study demands no further information regarding the material’s constitutive behavior.

This paper is organized as the following. In Sec. 2, analysis of finite deformation is briefly summarized. Here, we emphasize all the quantities that can be experimentally measured in a 2D deformation. The connection between the principal stretch orientation and the principal stress orientation for a 2D deformation is discussed in Sec. 3. In Sec. 4, we discuss the equation of equilibrium in terms of the Cauchy principal stresses and the principal stress orientation. We show that when the field of the principal stress orientation is known, the equation of equilibrium can be solved for the two principal stresses with appropriate traction boundary conditions. Some of the unique characteristics of the partial differential equations are discussed. The numerical scheme for solving the PDEs discussed in Sec. 4, is the topic of Sec. 5. In Sec. 6, we present two virtual testing cases to verify and validate the technique and the numerical scheme discussed in Secs. 4 and 5. The two virtual testing cases include the Brazilian disk compression (Sec. 6.1), and an infinite sheet with a circular hole and subject to remote tension (Sec. 6.2). An experimental case, the uniaxial tension of a polyvinyl chloride (PVC) sheet with a circular hole is presented in Sec. 7, where the stress field determined using the technique and the numerical scheme developed in the present study is also compared with the prediction using a constitutive model of hyperelasticity, the GNH model, calibrated for the particular PVC sheet. Concluding remarks are made in Sec. 8.

When the author is ready to submit this article to a journal, we learned the new publication by Cameron and Tasan [7], where they addressed the same questions considered in this study. Their approach is based on what they referred to as the “alignment assumption,” where the principal directions of stress align with that of strain. Our approach relies on the mathematical notion of tensor-valued isotropic functions. The subsequent treatment of the equation of equilibrium is also different, in Ref. [7], the Cartesian form of the equation is maintained, while in the present study, we recast the equation into the curvilinear system formed by the two families of characteristic curves and by doing so, the numerical solution to the equation of equilibrium is believed to be simpler.

## 2 Finite Deformation

**, in the undeformed state, to the new position**

*x***in the deformed state. In the present investigation, we focus on planar deformation. Consider a continuous planar object, which occupies a closed two-dimensional (2D) region $Rref$. The motion of the 2D object can be written as**

*y**y*

_{α}is the position vector at the current moment during the motion of a particle that is located at

*x*

_{α}in the reference region $Rref$. We can also write

*u*

_{α}(

*x*

_{1},

*x*

_{2}) is the displacement at the current moment of the particle located at

*x*

_{α}in the reference region. In defining the notion of motions, we require that $y^\alpha (x1,x2)$ is twice differentiable continuous over region $Rref$, the mapping $y^\alpha $ is one-on-one, and the Jacobian $J(x1,x2)=det{y^\alpha ,\beta (x1,x2)}>0$ for all $x\alpha \u2208Rref$ (

*α*,

*β*= 1, 2).

*b*

_{α}is a vector and

*F*

_{αβ}is a tensor dependent on location (

*x*

_{1},

*x*

_{2}). Since $J(x1,x2)=det{y^\alpha ,\beta (x1,x2)}=det{F\alpha \beta}$, necessarily $det{F\alpha \beta}>0$. Meanwhile, the deformation gradient

*F*

_{αβ}can be related to the displacement

*u*

_{α}through

*x*

_{1},

*x*

_{2}) and with the orientation angle

*θ*. The original length of the material fiber is

*δ*ℓ, and it becomes

*δL*after deformation. The ratio

*δL*/

*δ*ℓ represents the stretch

*λ*of the material fiber at location (

*x*

_{1},

*x*

_{2}) and along the direction

*θ*. When

*λ*> 1, the material fiber is elongated, and when

*λ*< 1, the material fiber is shortened. The stretch

*λ*can be expressed in terms of the deformation gradient

*F*

_{αβ}through

*λ*is a function of location (

*x*

_{1},

*x*

_{2}) and angle

*θ*. The extension of a material segment in the direction of

*m*

_{α}, is simply

*λ*− 1.

*x*

_{1},

*x*

_{2}) but varying the angle

*θ*, we find that at some angle the stretch reaches extrema, maximum, or minimum. The directional angle,

*θ*

_{λ}, at which the stretch

*λ*reaches extrema can be determined by

*θ*

_{λ}are the principal stretches. In planar deformation, the direction perpendicular to (

*x*

_{1},

*x*

_{2}) plane is always a principal direction. The other two principal stretches

*λ*

_{1}and

*λ*

_{2}, where

*λ*

_{1}⩾

*λ*

_{2}, are within the (

*x*

_{1},

*x*

_{2})-plane, and we assign the

*θ*

_{λ}to be measured from the

*x*

_{1}axis to the major stretch

*λ*

_{1}counterclockwise.

From experimental perspective, all the quantities discussed so far can be measured in an experiment over the sample domain $Rref$. For example, using the optical technique of digital image correlation (DIC), one can measure the displacement field *u*_{α} over the region $Rref$. Consequently, deformation gradient *F*_{αβ} can be computed through differentiation of the displacement field *u*_{α} according to Eq. (4). The stretch *λ* can be determined at any given location (*x*_{1}, *x*_{2}) and for any arbitrary angle *θ* from Eq. (5), and this will lead to the determination of the principal stretches *λ*_{1} and *λ*_{2}, and the principal stretch orientation *θ*_{λ} at any given location (*x*_{1}, *x*_{2}).

## 3 Principal Stretch Orientation and Principal Stress Orientation

In this section, we establish the relationship between the principal stretch orientation *θ*_{λ}, defined in Sec. 2, and the principal stress orientation *θ*_{σ} in the deformed body.

### 3.1 Strains in Terms of Principal Stretches and Principal Stretch Orientation.

**is an arbitrary vector. From polar decomposition [8], we have**

*b***is a proper orthogonal tensor, i.e., $Q\u22121=Q\u22a4$ and $detQ=1$. Finally, $V=QUQ\u22a4$.**

*Q**u*

_{α}(

*x*

_{1},

*x*

_{2}) (

*α*= 1, 2) can be experimentally measured, so is the deformation gradient tensor

*F*

_{αβ}. Meanwhile, the principal stretches,

*λ*

_{1}and

*λ*

_{2}, and the principal stretch orientation

*θ*

_{λ}, can be determined as well. Since the principal stretches,

*λ*

_{1}and

*λ*

_{2}, are the two eigen values of the right-stretch tensor

**, the two eigenvectors associated with**

*U**λ*

_{1}and

*λ*

_{2}are

*x*

_{1},

*x*

_{2}), the right-stretch tensor

**can be written as**

*U***(**

*R**θ*

_{λ}) is a proper orthogonal tensor, as well as a rotational tensor associated with angle

*θ*

_{λ}. Meanwhile, the inverse of the right-stretch tensor

**can be written through the principal stretches,**

*U**λ*

_{1}and

*λ*

_{2}, and the principal stretch orientation

*θ*

_{λ}as

**can be written as**

*Q***represents a rigid-body rotation, and we may express the orthogonal tensor**

*Q***as a rotation tensor**

*Q***(**

*R**θ*

_{q}), where the rotation angle

*θ*

_{q}is completely determined once the deformation gradient

*F*

_{αβ}, the principal stretches

*λ*

_{1}and

*λ*

_{2}, and the principal stretch orientation

*θ*

_{λ}are experimentally measured.

**, is defined by**

*E***has two eigen values, $(\lambda 12\u22121)/2$ and $(\lambda 22\u22121)/2$, with the two eigenvectors**

*E*

*N*_{1}and

*N*_{2}given in Eq. (9). Therefore, the principal stretch orientation

*θ*

_{λ}also defines the direction of the major principal strain, $(\lambda 12\u22121)/2$, of the Lagrangian strain

**.**

*E*To summarize, from experiment, we can measure the displacement field *u*_{α}(*x*_{1}, *x*_{2}) (*α* = 1, 2), the deformation gradient tensor ** F**, the principal stretches,

*λ*

_{1}and

*λ*

_{2}, and the principal stretch orientation

*θ*

_{λ}, and consequently, we can also determine the components involved in the polar decomposition of the deformation gradient tensor

**,**

*F***, and**

*U***.**

*Q*### 3.2 Notion of Tensor-Valued Isotropic Functions.

**and**

*A***be two second-order and symmetric tensors, and let**

*B***be a tensor-valued function such that**

*φ***is called isotropic if for every proper orthogonal tensor**

*φ***[9],**

*Q***represents a rotation, the tensor-valued isotropic function**

*Q***is invariant under arbitrary rotations. One consequence of an isotropic tensor-valued function**

*φ***that relates**

*φ***and**

*A***is that tensors**

*B***and**

*A***are co-axial, or they share the same (orthonormal) eigenvectors. Also note that tensors**

*B***and**

*A***are co-axial if and only if**

*B***=**

*AB***.**

*BA*A material is isotropic if it exhibits the same mechanical properties in all directions. In terms of stress and strain at any location within the deforming body, a material is isotropic if the relation between stress and strain at that location is independent of directions, though this relation may still be a function of locations, e.g., for inhomogeneous materials. The stress-strain relation at any point can be represented as a tensor-valued function and for isotropic materials, such tensor-valued function must be isotropic as well. However, not any pair of stress and strain is related by such a function; here, we only consider the stress and strain that form the so-called work-conjugate pair [10,11]. Some specific representations of the isotropic tensor-valued function ** φ** for certain isotropic materials can be found in Ref. [11]. Based on the notion of isotropic tensor-valued function, for isotropic materials such pair of stress tensor and strain tensor are, therefore, always co-axial.

### 3.3 Work Conjugate and the Second Piola–Kirchhoff Stress.

**is the second Piola–Kirchhoff stress**

*E***, where $S\u22a4=S$. Suppose that stress**

*S***has the principal components**

*S**s*

_{1}and

*s*

_{2}, where

*s*

_{1}⩾

*s*

_{2}, and the major principal stress

*s*

_{1}is in the direction described by the angle

*θ*

_{s}measured from

*x*

_{1}-axis to

*s*

_{1}counterclockwise. The components of

**in the Cartesian coordinate (**

*S**x*

_{1},

*x*

_{2}) can be related to the principal stresses

*s*

_{1}and

*s*

_{2}by

**(**

*R**θ*

_{s}) is the rotation tensor associated with angle

*θ*

_{s}.

*isotropic*materials, the principal stress orientation

*θ*

_{s}of the second Piola–Kirchhoff stress

**, is the same as the principal stretch orientation**

*S**θ*

_{λ}, or the principal orientation of the Lagrangian strain

**. Since both**

*E***and**

*S***are expressed in the reference or undeformed configuration,**

*E***and the second Piola–Kirchhoff stress**

*E***for isotropic materials. As we have mentioned in Sec. 3.1 that the principal stretch orientation**

*S**θ*

_{λ}(

*x*

_{1},

*x*

_{2}) can be measured experimentally over the reference domain $Rref$, as a result, the principal stress orientation

*θ*

_{s}(

*x*

_{1},

*x*

_{2}) of the second Piola–Kirchhoff stress

**is also experimentally measurable over domain $Rref$.**

*S*### 3.4 Cauchy Stress and Principal Stress Orientation.

**. The Cauchy stress**

*σ***is related to the second Piola–Kirchhoff stress**

*σ***, through the deformation gradient tensor**

*S***, by**

*F***and the second Piola–Kirchhoff stress**

*E***, i.e.,**

*S**θ*

_{s}=

*θ*

_{λ}, the polar decomposition

**=**

*F***, and**

*QU***=**

*Q***(**

*R**θ*

_{q}). Furthermore, using Eq. (10),

**(**

*R**θ*

_{q})

**(**

*R**θ*

_{λ}) =

**(**

*R**θ*

_{q}+

*θ*

_{λ}), has been used in Eq. (20). Note that $\Lambda S^\Lambda /J$ is a diagonal tensor, so Eq. (20) provides one form of diagonalization of the Cauchy stress tensor

**for a planar or 2D deformation.**

*σ***has the principal components**

*σ**σ*

_{1}and

*σ*

_{2}, where

*σ*

_{1}⩾

*σ*

_{2}, with the principal stress orientation

*θ*

_{σ}, which is measured from the

*y*

_{1}-axis to

*σ*

_{1}counterclockwise, we will have

*θ*

_{σ}of a material particle in the deformed state equals to the principal stretch orientation

*θ*

_{λ}of the same material particle in the undeformed state plus the rotation angle

*θ*

_{q}due to the rigid-body rotation of the material particle associated with

**, i.e.,**

*Q**α*= 1, 2). Now, we have established the relationship between the principal stress orientation

*θ*

_{σ}of a material particle in the deformed state and the principal stretch orientation

*θ*

_{λ}of the same material particle in the reference state within a domain made of an isotropic material. Since in Eq. (22), both the principal stretch orientation

*θ*

_{λ}and the rotation angle

*θ*

_{q}due to local rigid-body rotation of the material particle, can be measured experimentally, the principal stress orientation

*θ*

_{σ}can be experimentally determined as well. In the following sections, we will explore a scheme for determining the nonuniform stress field

**of a domain $R$ when the field of the principal stress orientation**

*σ**θ*

_{σ}is given within domain $R$ and the traction boundary conditions are described over the boundary of domain $R$.

## 4 Equation of Equilibrium, Characteristic Lines, and Recast Partial Differential Equation

### 4.1 Equation of Equilibrium.

*y*

_{1},

*y*

_{2}), together with the traction boundary condition, is given by

**, where $\sigma \u22a4=\sigma $. Domain $R$ represents the deforming region in the current configuration and its boundary $\u2202R$ has the unit outward normal vector**

*σ***.**

*n***in the 2D region $R$,**

*σ**σ*

_{1}and

*σ*

_{2}, and with the associated principal stress orientation

*θ*

_{σ}, measured from

*y*

_{1}-axis to the major principal stress

*σ*

_{1}counterclockwise. In the Cartesian coordinate (

*y*

_{1},

*y*

_{2}), the Cauchy stress

**can be expressed in terms of its principal components**

*σ**σ*

_{1}and

*σ*

_{2}, and the principal stress orientation

*θ*

_{σ}by Eq. (21), or to be more explicit,

*θ*

_{σ}can be related to the principal stretch orientation

*θ*

_{λ}through Eq. (22), where both the principal stretch orientation

*θ*

_{λ}and the rotation angle

*θ*

_{q}due to local rigid-body rotation can be measured experimentally. Therefore, the principal stress orientation

*θ*

_{σ}is also experimentally determinable over the current domain $R$.

Such experimentally measured field *θ*_{σ}(*y*_{1}, *y*_{2}) imposes a constraint to the PDEs in Eq. (26), where we have two unknown functions, the principal stresses *σ*_{1} and *σ*_{2}, and two PDEs. With appropriate traction boundary conditions, the stress field, (*σ*_{1}, *σ*_{2}), which in general is nonuniform, can be determined through numerical solutions to the PDEs in Eq. (26) over the given domain $R$.

### 4.2 Characteristic Curves and Recast Partial Differential Equations.

*σ*

_{1}as the

*α*-lines, and the ones associated with the minor principal stress

*σ*

_{2}as the

*β*-lines. Once the field

*θ*

_{σ}(

*y*

_{1},

*y*

_{2}) is given, such two-family of characteristic curves,

*α*-lines and

*β*-lines, are completely determined.

Note that at each location in domain $R$, the *α*-line and the *β*-line passing through that location are perpendicular to each other since the product of their slopes, *dy*_{2}/*dy*_{1}, at that location equals to −1. As a result, we may prefer the “curvilinear coordinates,” formed by the *α*- and the *β*-lines, over the more conventional Cartesian coordinates (*y*_{1}, *y*_{2}). Unlike the Cartesian coordinates (*y*_{1}, *y*_{2}), where there is a well defined origin *y*_{1} = *y*_{2} = 0, along each *α*-line and at any given location the coordinate *s*_{α} is given by the length, usually, from one end of the *α*-line to that location. Similarly coordinate along the *β*-line, *s*_{β}, can be defined. At the intersection of an *α*-line and a *β*-line, the coordinate (*s*_{α}, *s*_{β}) is thus given in such manner. One final note is that similar to the Cartesian coordinate (*y*_{1}, *y*_{2}), the curvilinear coordinate (*s*_{α}, *s*_{β}) should form a right-handed system.

*s*

_{α},

*s*

_{β}), consider the derivative of an arbitrary variable

*f*along an

*α*-line and along a

*β*-line,

*s*

_{α},

*s*

_{β}), the recast PDEs in Eq. (26) become

*θ*

_{σ}(

*s*

_{α},

*s*

_{β}) in the curvilinear coordinate (

*s*

_{α},

*s*

_{β}).

Equation (28) represents the simplest form of the equation of equilibrium in a 2D deformation. Our subsequent exploration of solving the equation of equilibrium, together with some traction boundary conditions, for the principal stresses, *σ*_{1} and *σ*_{2} once the functions $\theta \sigma ,s\alpha (s\alpha ,s\beta )$ and $\theta \sigma ,s\beta (s\alpha ,s\beta )$ are given, will be based on Eq. (28).

## 5 Numerical Scheme

*σ*

_{2}is known as a function of the curvilinear coordinate (

*s*

_{α},

*s*

_{β}) along a given

*α*-line, say

*L*

_{α}, the PDE can be written as

*s*

_{α},

*s*

_{β}) ∈

*L*

_{α}. Note that the first-order linear PDE in the above equation only involves a single-derivative term, ∂

*σ*

_{1}/∂

*s*

_{α}, and thus, it is integrable along the characteristic

*α*-line and one can formally write the expression for

*σ*

_{1}as

*c*

_{α}(

*s*

_{β}) is a integration constant but may depend on

*s*

_{β}and must be determined by the appropriate boundary condition at either ends of the specific

*α*-line

*L*

_{α}.

*σ*

_{1}is known as a function of the curvilinear coordinate (

*s*

_{α},

*s*

_{β}) in the second PDE in Eq. (28) along a given

*β*-line,

*L*

_{β}, we will have

*s*

_{α},

*s*

_{β}) ∈

*L*

_{β}, and the formal expression for

*σ*

_{2}can be written as

*c*

_{β}(

*s*

_{α}) is another integration constant but may depend on the curvilinear coordinate

*s*

_{α}.

Although Eqs. (30) and (32) indicate both the existence and the form of solutions for the two principal stresses, *σ*_{1} and *σ*_{2}, under the given assumptions, for practical applications, we will not directly use Eq. (30) or (32). We will use a numerical form of Eqs. (29) and (31) to obtain the principal stresses, *σ*_{1} and *σ*_{2}, at discrete locations inside the deformed domain $R$.

*θ*

_{σ}(

*y*

_{1},

*y*

_{2}) and according to Eq. (27), we already determined the two family of characteristic curves

*α*-line $L\alpha (i)$, the coordinate is described by (

*s*

_{α},

*s*

_{β}) and in particular $0\u2a7ds\alpha \u2a7dS\alpha (i)$, where $S\alpha (i)$ is the total length of the

*α*-line $L\alpha (i)$. Similarly, coordinate

*s*

_{β}along a

*β*-line $L\beta (j)$ is such that $0\u2a7ds\beta \u2a7dS\beta (j)$, where $S\beta (j)$ is the total length of the

*β*-line $L\beta (j)$. If characteristic curves $L\alpha (i)$ and $L\beta (j)$ intersect, the intersecting point is described by the coordinates $(s\alpha (i,j),s\beta (i,j))$. In general, we have

*M*×

*N*such intersecting points over the deformed domain and it is our objective to solve the equation of equilibrium (28) and to obtain numerical values of the principal stresses,

*σ*

_{1}and

*σ*

_{2}, at every intersecting point $(s\alpha (i,j),s\beta (i,j))$. Note that in some domain configurations, e.g., multiply-connected domain, for a given pair (

*i*,

*j*), the characteristic lines $L\alpha (i)$ and $L\beta (j)$ may not intersect, and as a result, we will have less than

*M*×

*N*intersecting points.

*α*-line family, and assume that $L\alpha (1)$ is part of the domain boundary. We also assume that the traction boundary condition is such that the principal stress

*σ*

_{2}is given over $L\alpha (1)$, denoted as $\sigma 2=\sigma 2(1)(s\alpha ,s\beta )$. Note that the traction boundary condition over $L\alpha (1)$ cannot be described in terms of the principal stress

*σ*

_{1}. This is related to the restrictions of the so-called Cauchy problem and Cauchy data. Now, we have the equation for

*σ*

_{1}over the one-dimensional line $L\alpha (1)$ (from the first PDE in Eq. (28)),

*α*-line $L\alpha (1)$, i.e.,

*s*

_{α}= 0 or $s\alpha =S\alpha (1)$, the

*α*-line $L\alpha (1)$ will intersect

*β*-lines, $L\beta (1)$ and $L\beta (N)$. We assert that both $L\beta (1)$ and $L\beta (N)$ also form part of the boundary of domain $R$, where the principal stress

*σ*

_{1}is prescribed. In other words, the above equation has the boundary condition at both ends of characteristic line $L\alpha (1)$ for

*σ*

_{1}; thus, it can be numerically solved to have the value of

*σ*

_{1}at discrete locations $(s\alpha (1,j),s\beta (1,j))$ for

*j*= 1, 2, …,

*N*. The value of $\sigma 2=\sigma 2(1)(s\alpha ,s\beta )$ at those locations is already known since it is part of the boundary condition.

*i*th

*α*-line, $L\alpha (i)$, the principal stresses,

*σ*

_{1}and

*σ*

_{2}, have been obtained, and they are given at discrete locations $(s\alpha (i,j),s\beta (i,j))$ as

*i*+ 1)th

*α*-line, $L\alpha (i+1)$, the value for principal stress

*σ*

_{2}can be estimated or updated to be

*σ*

_{2}over line $L\alpha (i+1)$, denoted as $\sigma 2=\sigma 2(i+1)(s\alpha ,s\beta )$, is known and the principal stress

*σ*

_{1}satisfies the following equation over line $L\alpha (i+1)$,

*j*= 1, 2, …,

*N*.

Such a process is thus repeated for all characteristic lines $L\alpha (i)$, *i* = 1, 2, …, *M*. The principal stresses, *σ*_{1} and *σ*_{2}, are therefore obtained at all locations $(s\alpha (i,j),s\beta (i,j))$, where *i* = 1, 2, …, *M* and *j* = 1, 2, …, *N*.

We can switch the role of *α*-lines and *β*-lines and start from the first *β*-line $L\beta (1)$, where the value of the principal stress *σ*_{1} is known due to the fact that $L\beta (1)$ is part of the boundary of the domain $R$ and the traction boundary condition for *σ*_{1} is given. Then, the principal stress *σ*_{2} can be solved by the equation over the line $L\beta (1)$. Using the similar updating scheme and repeating for all the *β*-lines $L\beta (j)$ (*j* = 1, 2, …, *N*), the principal stresses, *σ*_{1} and *σ*_{2}, are therefore obtained at all locations $(s\alpha (i,j),s\beta (i,j))$, where *i* = 1, 2, …, *M* and *j* = 1, 2, …, *N*.

This numerical scheme will be verified and validated through two virtual testing cases. The factors that might affect the convergency and accuracy of the numerical results will be identified and discussed.

## 6 Verification and Validation

Two “virtual” testing cases are considered in this section to verify and to validate the technique and the numerical scheme discussed in Secs. 4 and 5. One testing case is the Brazilian disk compression, where the deforming domain is simply connected. The second testing case is an infinite sheet with a circular hole and subject to remote tension, and the deforming domain in this case is multiply connected. The numerically determined stress field is then compared with the closed-form solution for the two testing cases.

### 6.1 Virtual Test Case I: Brazilian Disk Compression.

The first testing case is the Brazilian disk compression, which has been used for observing and measuring mechanical properties, damage/cracking initiation, and crack extension in many different materials [12,13], for the reason that closed-form solution for isotropic elastic material exists.

*R*, subjects to a pair of concentrated forces

*F*

_{2}along its vertical diameter, as shown in Fig. 1(a). For such a loading configuration, the stress field can be expressed in a simple algebraic form listed below [14]:

*θ*

_{σ}can be expressed as

*θ*

_{σ}(

*x*

_{1},

*x*

_{2}) is shown in Fig. 1(b). With the principal stress orientation field

*θ*

_{σ}given, the two family of the characteristic curves can be determined through the ODEs given in Eq. (27), which can be specialized for the Brazilian disk so that the

*α*-lines are given by

*β*-lines are given by

*a*and

*b*are two parameters, where

*a*is the vertical coordinate of the intersecting point of an

*α*-line with the

*x*

_{2}-axis and

*b*is the horizontal coordinate of the intersecting point of a

*β*-line with the

*x*

_{1}-axis. Both

*α*-lines and

*β*-lines are portions of circles, and they are also shown in Fig. 1(b), where the

*α*-lines are red (or running left and right) and the

*β*-lines are blue (or runing up and down).

The Brazilian disk compression test has a characteristic feature that all three stress components vanish along the disk boundary except at the two ends of the vertical diameter. Along the traction-free boundary, since the shear stress is already zero, one of the principal stresses must be zero as well. Such principal stress can be identified by tracking the variation of the principal stress orientation *θ*_{σ} along the boundary. There are only two choices for the principal stress orientation *θ*_{σ} along any traction-free boundary, either it is pointing to the normal direction of the boundary or it is pointing to the tangential direction of the boundary. If *θ*_{σ} is in the normal direction of the boundary, then *σ*_{1} = 0. As shown in both plots in Fig. 1(b), we should have the boundary condition *σ*_{1} = 0 everywhere along the circular edge of the disk, except at the two poles and at the two poles, the minor principal stress *σ*_{2} = *σ* = −*F*_{2}/*R*.

*θ*

_{σ}in the curvilinear coordinate (

*s*

_{α},

*s*

_{β}), formed by the two family of characteristic curves, can be expressed as

*σ*

_{1}and

*σ*

_{2}.

*α*-line

*L*

_{α}(

*a*

_{1}), which is very close to the loading point near the bottom of the disk. As shown in Fig. 3(a), the minor principal stress

*σ*

_{2}can be approximated by using the Boussinesq solution or

*r*(

*θ*) is the distance from the bottom pole to any location on

*L*

_{α}(

*a*

_{1}) and

*θ*is the polar angle shown in Fig. 3(a). Once

*σ*

_{2}is known over

*L*

_{α}(

*a*

_{1}), Eq. (33) is used for solving

*σ*

_{1}over

*L*

_{α}(

*a*

_{1}), with the boundary condition

*σ*

_{1}(

*s*

_{α}= 0,

*s*

_{β}) = 0. Equations (34) and (35) are used for solving the principal stresses,

*σ*

_{1}and

*σ*

_{2}, over the subsequent

*α*-lines

*L*

_{α}(

*a*

_{i}) (

*i*= 1, 2, …,

*M*). The update of the principal stress

*σ*

_{2}from line

*L*

_{α}(

*a*

_{i}) to line

*L*

_{α}(

*a*

_{i+1}) is sketched in Fig. 3(b).

We can also switch the role of *α*- and *β*-lines and start from the first *β*-line *L*_{β}(*b*_{1}), where *σ*_{1} ≡ 0 according to the traction boundary condition. The minor principal stress *σ*_{2} over *L*_{β}(*b*_{1}) can be obtained by solving the integrable PDE with the end condition given by the Boussinesq solution. The update of the major principal stress *σ*_{1} along *α*-lines *L*_{α}(*a*) will lead to the value of *σ*_{1} along the next *β*-line *L*_{β}(*b*) similar to that shown in Fig. 3(b). This procedure is then repeated for all *β*-lines. Such a scheme is illustrated in Fig. 3(c).

Figure 4(a) presents the principal stress fields, *σ*_{1} and *σ*_{2}, from the closed-form solution. The numerical results of the two principal stresses, *σ*_{1} and *σ*_{2}, according to the proposed scheme, indicated in Fig. 3(b), are presented in Fig. 4(b). The top plot in Fig. 4(b) shows the result where 101 *α*-lines and 41 *β*-lines are used, and the maximum length of the line segment *δs*_{β}, used in the *σ*_{2} updating scheme, is about 2% of the disk radius *R*. The comparisons presented in Fig. 4(b) highlight the important role of the maximum length of the line segment *δs*_{β}, or step size, in the proposed numerical scheme. Only when the step size *δs*_{β} is small enough that the numerical solutions of the two principal stresses, *σ*_{1} and *σ*_{2}, converge to the closed-form solution, as shown in the bottom plot of Fig. 4(b), where max (*δs*_{β}) ∼ 0.004*R*. The small notches shown in Fig. 4(b) are the first and the last *α*-lines, where the concentrated force boundary condition is approximated by the Boussinesq solution.

The numerical results of the two principal stresses, *σ*_{1} and *σ*_{2}, according to the alternative scheme shown in Fig. 3(c), are presented in Fig. 4(c). In the top plot of Fig. 4(c), we only use very few data points, but we see that the numerical result is already quite good, though slightly off symmetry in *σ*_{1}. In the middle and the bottom plots of Fig. 4(c), we moderately increase the number of data points, and the numerical results become much better. Also, we see that the smallest step size is only about 2% of the disk radius *R*, but in the previous updating scheme, we have to use much smaller step size to reach the same accuracy.

The virtual testing case of Brazilian disk compression is unique in that every quantity for checking the proposed scheme can be written in closed form and these include each individual characteristic curve. As a result, it can isolate the factors that might affect the accuracy of the numerical scheme easily. We have shown that using the curvilinear coordinate formed by the two-family of characteristic curves makes the PDE and boundary condition much simpler and cleaner, since the boundary is one of the characteristic curves. The determination of the characteristic curves only requires the experimentally obtainable principal stress orientation *θ*_{σ} and based on the curvilinear coordinate system, we can easily convert solving the coupled PDEs to solving a series of integrable PDEs along one family of the characteristic curves.

There are two factors that control the accuracy and efficiency of the numerical scheme. The first factor is the approximation the point-force boundary condition by a traction along the *α*-line *L*_{α}(*a*_{1}), using the Boussinesq solution. By choosing, the *L*_{α}(*a*_{1}) very close to the loading point will reduce the approximation error. The second factor is the finite difference scheme used for updating *σ*_{2} from *L*_{α}(*a*_{i}) to *L*_{α}(*a*_{i+1}), or updating *σ*_{1} from *L*_{β}(*b*_{j}) to *L*_{β}(*b*_{j+1}), where the accuracy of the resulting stress depends on the step size. Finally, which family of characteristic curves to choose for solving the integrable PDEs also affects the accuracy and efficiency.

### 6.2 Virtual Test Case II: Infinite Sheet With a Circular Hole and Subject to Remote Tension.

In the second testing case, we consider a multiply-connected domain: an infinite sheet with a circular hole and subject to remote tension, as shown in Fig. 5, where the radius of the circular hole is *a*, and the remote tensile stress is *σ* and in the *x*_{2} direction. The sheet is assumed to be made of homogeneous, isotropic, and linearly elastic material. Once again, for infinitesimal deformation, there is no need to distinguish between the reference and the deformed configurations.

*x*

_{1},

*x*

_{2}) is listed below [15]

*a*)

*b*)

*c*)

The field of the principal stress orientation *θ*_{σ} is shown in Fig. 6(a) for the virtual testing case. For the majority of the domain $R$, *θ*_{σ} is around 90 deg, only near the edge of the circular cavity, *θ*_{σ} varies away from 90 deg. Since we already knew/measured the principal stress orientation *θ*_{σ}, the characteristic curves can be obtained by solving the two ODEs in Eq. (27), and they are presented also in Fig. 6(a). Again, we assign the *α*-lines (red or running up and down) to be associated with the principal stress *σ*_{1} and the *β*-lines (blue or running left and right) to principal stress *σ*_{2}. These characteristic curves have some characteristics, especially around the circular cavity, as shown in Fig. 6(b).

The *β*-line far below the hole is approximately running from left to right and the farther away from the hole, the flatter the *β*-line is. As the *β*-line approaches the circular hole, which is traction-free along its boundary, the middle section of the *β*-line is “attracted” toward the circular hole and bent. At some point, the *β*-line touches the edge of the circular hole and wraps around a section of the cavity, so that the lower section of the circular cavity becomes part of the characteristic *β*-line. This section ends at two points, where the polar angles are *θ* = −60 deg and −120 deg, respectively, at which the *β*-line has sharp kinks and then turns away from the hole. Further up, the *β*-line becomes two branches on either side of the cavity, and they can be thought of being semi-infinite, starting from the edge of the circular cavity and extending to infinity. Also at points that the *β*-lines touch the hole edge, these *β*-lines are all perpendicular to the hole boundary. At locations where the polar angles *θ* = 60 deg and 120 deg, the two-branched *β*-lines merges into one and again, wraps around portion of the cavity boundary. Afterward, the *β*-line departs the hole and moves upward to infinity.

Similar observations can be made of the characteristic *α*-lines, from the right side of the hole to the left side of the hole. The *α*-line also touches the hole edge and wraps around the hole boundary between the polar angles *θ* = ±60 deg, and then becomes two branches, and finally merges into one again. Note that the characteristic curves around the circular cavity boundary are closely related to how the traction boundary conditions are described, i.e., the traction boundary condition is prescribed for the quantity whose characteristic curves are perpendicular to the boundary, not the one whose characteristic curves are part of the boundary.

*H*is chosen such that

*H*/

*a*≫ 1. Along the edge of the circular cavity, we have the traction-free boundary condition. By observing the characteristic curves near the edge of the circular hole, i.e., Fig. 6(c), we can specify the boundary condition along the edge of the circular cavity to be

The contour plots of the gradient of the principal stress orientation *θ*_{σ} in the curvilinear coordinate (*s*_{α}, *s*_{β}), formed by the two family of characteristic curves, are shown in Fig. 7.

The numerical scheme of solving for the principal stresses, *σ*_{1} and *σ*_{2}, is similar to that used in the Brazilian disk compression test. But, there are some differences due to the fact that the deforming domain is multiply connected. The similarities and differences of the numerical scheme are illustrated in Fig. 8 and discussed here.

*β*-line, $L\beta (1)$, as shown in the upper plot in Fig. 8, we have the approximated boundary condition

*σ*

_{1}=

*σ*, so that the minor principal stress

*σ*

_{2}satisfies the PDE,

*σ*

_{2}is solved along the first

*β*-line, $L\beta (1)$. Now suppose that

*σ*

_{1}and

*σ*

_{2}along

*β*-line $L\beta (j)$ are known, we are seeking the solution along

*β*-line $L\beta (j+1)$. The major principal stress

*σ*

_{1}can be updated along the characteristic

*α*-lines as

*σ*

_{2}satisfies the PDE with some initial conditions,

*σ*

_{2}is solved over the

*β*-line $L\beta (j+1)$. This step is illustrated in the middle plot in Fig. 8, and it can be repeated for all characteristic

*β*-lines.

However, the virtual test of an infinite sheet with a circular hole and subject to remote tension deals with multiply-connected domain. At some point, portion of the hole boundary becomes part of the *β*-line, where the next *β*-line has two branches, as illustrated in the bottom plot in Fig. 8. As a result, when we update the principal stress *σ*_{1} along *α*-lines, only solutions (*σ*_{1}, *σ*_{2}) along a portion of the previous *β*-line are used, while the solution (*σ*_{1}, *σ*_{2}) over part of the circular boundary is not used. The solution update scheme is carried out on the left and the right sides of the circular cavity, as seen in the bottom plot in Fig. 8, until the moment when the next *β*-line becomes a single line extending to infinity on both sides. The principal stress *σ*_{1} on the left and the right sides of the *β*-line is updated according to our scheme. However, according to the traction boundary condition, along the boundary of the circular hole, which is also part of the characteristic *β*-line, we have *σ*_{1} = 0. Therefore, the principal stress *σ*_{1} along the entire *β*-line is known. Now, the scheme illustrated in the middle plot in Fig. 8 can be applied for the rest of the domain $R$.

We now compare the numerical results with the closed-form solution for the virtual test problem, where the closed-form solution for the two principal stresses is shown in Fig. 9.

We first study the effect of data density or mesh size in the region near the circular hole. We fix the far-field boundary at the distance of *H* = 4*a*. The top plot in Fig. 10 shows the results using very few data points, less than two dozens on each side, and we see that the numerical solution already captures the essential structure of the stress field. In the middle plot in Fig. 10, the data density around the circular hole increases to a moderate level, to about 50 data points on each side, and the numerical solution captures more details of the stresses along the hole boundary. When the data density around the hole increases more, as shown in the bottom plot in Fig. 10, the stress field surrounding the hole improves but no drastic improvement can be seen. This is more evident by comparing the value of the principal stress *σ*_{2} from the closed-form solution, shown in Fig. 9, near the left and the right sides of the circular hole with those from the numerical solutions.

The more important factor that affects the accuracy of the numerical solution is the parameter *H*, i.e., at which location we approximate the far-field boundary condition. In Fig. 11, we fix the data density in the region around the circular hole and change the parameter *H* from relatively close to the hole, *H* = 4*a*, to relatively farther away from the hole, *H* = 8*a*. Since the computational domain enlarges, the number of data points also increases. But we see that as we push the far-field boundary further away, the quality of the numerically obtained stress field improves significantly. This can be seen by looking at the two regions on either side of the circular hole, where the principal stress *σ*_{2} is positive. As *H* becomes larger, these two regions become closer to that in the closed-form solution.

In many experimental measurements, the experimental sample design involves multiply-connected domain, for example, specimen with an internal crack and specimen with damage characterized by cavity of certain shape. The virtual test of an infinite sheet with a circular hole and subject to remote tension represents such class of specimen designs and test configurations. By tracking the variation of the experimentally measurable principal stress orientation *θ*_{σ}, the traction boundary condition along the edge of the cavity can be identified and such boundary conditions are exact. The far-field boundary condition can be approximated through remotely applied total force. Through careful sample design and the placement of far-field boundary, error due to approximated far-field boundary condition can be reduced and the accuracy of the numerical results improved, as shown in Fig. 11. However, the characteristic curves associated with the multiply-connected domain become more complicated surrounding the cavity and denser data grid is required.

We have transformed the solving of a set of PDEs into solving a series of integrable PDEs along one family of the characteristic curves and then updating the solution along the other family of the characteristic curves. When solving an integrable PDE along a given characteristic curve, it only involves “initial” condition, or boundary condition from one side of the curve. As a result, advantage can be gained by starting the numerical scheme from both sides of the domain, top and bottom or left and right, and the final results then meet somewhere in the middle. Such an approach will minimize the error when the finite difference scheme, for updating the solution, moves around the cavity.

## 7 Uniaxial Tension of Polyvinyl Chloride Sheet With a Circular Hole

In this section, an experimental investigation of the nonuniform deformation and the associated nonuniform stress field is carried out on a PVC sheet with a circular hole and subject to uniaxial tension. We use the optical technique of DIC [16,17] to measure the nonuniform displacement field over the sample surface and then apply the technique and the numerical scheme, discussed in Secs. 3–5, to determine the nonuniform stress field in the sample. This stress field is then compared with the prediction using a hyperelasticity constitutive model, the generalized neo-Hookean (GNH) model, which has been calibrated through uniaxial tension test for the PVC sheet.

### 7.1 Uniaxial Tension Test, Overall Response, and Deformation Measurement.

The material used in this experiment, PVC is a thermoplastic polymer. The sketch of the PVC sheet specimen is shown in Fig. 12(a) and the nominal dimensions of the sample are the gage length *L* = 20.84 mm, the width *W* = 37.50 mm, the diameter of the circular hole *d* = 5.31 mm, and the thickness of the sheet *B* = 0.254 mm.

We used an Instron test machine (model 1125) to apply the tensile load to the specimen and the load cell capacity is 500 N. The loading speed was controlled by the crosshead velocity of the test machine and it was set at 1.2 mm/min, which gives the equivalent strain rate of about 1.0 × 10^{−3}s^{−1}. The overall response of one of the PVC sheet specimens is presented in Fig. 12(b), where the applied tensile load *P* is normalized by the cross-sectional area of the sheet *A* = *W B* and the gage section elongation *δ* is normalized by the original gage length *L*.

Optical DIC is used for mapping the displacement field over the PVC sheet sample surface. Such displacement measurement is shown in Fig. 13 for the selected moments indicated in Fig. 12(b). In the contour plots of Fig. 13, *u*_{1} represents the displacement in the horizontal direction and *u*_{2} in the vertical direction. The PVC sheet specimen experiences quite large deformation with the axial extension $\delta /L\u223c55%$ at moment (D), while at the same time, the cavity at the center of the sheet specimen, which starts as a circle, is elongated over 220% in the vertical direction and contracts about 8% in the lateral direction. Meanwhile, both the top and the bottom edges of the gage section remain approximately straight during the tensile process.

With the displacement field measured, deformation gradient can be computed through numerical differentiation. The deformation gradient field *F*_{αβ} (*α*, *β* = 1, 2) at moment (D) is shown in Fig. 14.

From the discussions in Sec. 2 and according to Eqs. (5) and (6), the fields of principal stretches, *λ*_{1} and *λ*_{2}, as well as the principal stretch orientation *θ*_{λ} can be determined. These fields are presented in Fig. 15 for the selected moments indicated in Fig. 12(b). For the isotropic PVC sheet, *θ*_{λ} also represents the principal stress orientation of the second Piola–Kirchhoff stress ** S** as discussed in Sec. 3.3.

The principal stretches, *λ*_{1} and *λ*_{2}, and the principal stretch orientation *θ*_{λ} will determine the right stretch tensor ** U** according to Eq. (10), and from Eq. (12), the orthogonal tensor

**can be obtained. For 2D deformation of the PVC sheet,**

*Q***represents a rigid-body rotation and can be represented by the rotation tensor**

*Q***(**

*R**θ*

_{q}), and consequently, the field of the rotation angle

*θ*

_{q}is determined. The principal stress orientation associated with the Cauchy stress

*θ*

_{σ}is related to the principal stretch orientation

*θ*

_{λ}and the rotation angle

*θ*

_{q}through Eq. (22). The field of the principal stress orientation

*θ*

_{σ}in the current or deformed configuration is shown in Fig. 16 for the selected moments indicated in Fig. 12(b).

### 7.2 Stress Field in Polyvinyl Chloride Sheet.

To determine the stress field by solving the equation of equilibrium, Eq. (28), the traction boundary conditions along the external and the cavity edges of the PVC sheet specimen are required. Along the external side edges of the sheet, we have traction-free boundary condition. By noting the principal stress orientation *θ*_{σ} along the side edges, we have the condition *σ*_{2} = 0. Note that the nominal applied load at the top and the bottom edges of the gage section is *σ*, but in the current configuration, this applied load is represented by $\lambda \xaf1\sigma $; here, $\lambda \xaf1$ is the average principal stretch *λ*_{1} along the top and the bottom edges of the gage section, as shown in Fig. 17(a).

Along the edge of the cavity, we also have traction-free boundary condition, but we need to identify which principal stress vanishes over which portion of the boundary. Figure 17(b) presents the variation of the principal stress orientation *θ*_{σ} as a function of the polar angle *θ* over the edge of the cavity. The polar angle *θ* is measured from the center of the cavity. In Fig. 17(b), the tangential direction of the cavity edge is shown as dashed line, and we see that for the most part, the principal stress orientation *θ*_{σ} tracks the tangential direction, which indicates that the minor principal stress *σ*_{2} = 0 over the portion of the cavity edge that *θ*_{σ} aligns with the tangential direction. Only between the two open circles for either *θ* < 0 or *θ* > 0, i.e., near the top and the bottom poles of the cavity, that the principal stress orientation *θ*_{σ} deviates from the tangential direction and aligns with the normal direction of the edge. Therefore, near the top and the bottom poles of the cavity, the boundary condition becomes *σ*_{1} = 0. However, the span of the polar angle at the top and the bottom poles of the cavity, where *σ*_{1} = 0, is only about 7.2 deg, so that the majority of the cavity edge is controlled by the boundary condition of *σ*_{2} = 0.

The other ingredient for solving Eq. (28) is the curvilinear coordinate system (*s*_{α}, *s*_{β}) formed by the two-family of characteristic curves obtained by solving the two ODEs in Eq. (27). For the principal stress orientation *θ*_{σ} field shown in Fig. 17(a) and at moment (D), the characteristic curves are given in Fig. 18(a). In Fig. 18(a), the red lines (running up and down) are the characteristic curves associated with the major principal stress *σ*_{1} (*α*-lines) and the blue lines (running left and right) are the ones associated with the minor principal stress *σ*_{2} (*β*-lines). One thing to note from Fig. 18(a) is that the characteristic curves (the *α*-lines) on the left and the right sides of the PVC sheet coincides with the specimen edge, where the PVC sheet specimen edge is shown as dashed line. This further confirms that the boundary condition *σ*_{2} = 0 along the lateral edges of the sheet specimen.

*θ*

_{σ}in the current configuration and the projection of the gradient into the curvilinear coordinate system (

*s*

_{α},

*s*

_{β}). Such gradient fields are shown in Fig. 18(b). In evaluating the gradient of the field shown in Fig. 17(a), it is convenient to use the following relation:

**is the deformation gradient with components shown in Fig. 14 as an example. From perspective of experiment, the data grid in the reference state is often regular, i.e., evenly spaced and aligned with the coordinate axes**

*F**x*

_{1}and

*x*

_{2}. As a result, the gradient $\u2207x\theta \sigma $ can be evaluated numerically in a straight forward fashion. In the current or deformed state, a straight line in the reference state may become curved, and therefore, Eq. (47) provides a convenient connection between gradients in the two states.

With all the ingredients ready, we follow the steps illustrated in Fig. 8 and use Eqs. (44)–(46) to determine the two principal stresses *σ*_{1} and *σ*_{2} in the PVC sheet. The computed principal stress fields (*σ*_{1}, *σ*_{2}) in PVC sheet specimen at moment (D) are presented in Fig. 19 for three different mesh sizes. We observe no significant differences when the mesh size decreases or the number of data point increases in the computation.

Finally, we would like to point out, again, that in determining the principal stress fields, *σ*_{1} and *σ*_{2}, as shown in Fig. 19, the only information required is the principal stress orientation field *θ*_{σ} shown in Fig. 17(a), together with the average stretch $\lambda \xaf1$ over the top and the bottom edges of the gage section. The principal stress orientation field *θ*_{σ} and the average stretch $\lambda \xaf1$ are the results from the deformation measurement. The principal stress fields, *σ*_{1} and *σ*_{2}, are then obtained by solving the equation of equilibrium over the domain shown in Fig. 18(a), with traction boundary conditions. Except the assumption of isotropy, no other knowledge of constitutive behavior of the deforming material gets involved.

### 7.3 Comparison With Constitutive Model Predictions.

We have mentioned in Sec. 1 that the usual approach for obtaining stresses associated with a nonuniform deformation field is to use a constitutive relation. Such a constitutive model needs to be calibrated for a particular material through independent experiment involving only simple loading and deformation configurations, or the stress state of the deforming sample in the calibration experiment is simple and can be calculated from applied force. In this section, we will follow such an approach to obtain the stresses associated with the deformation field of the PVC sheet shown in Figs. 13 and 15.

**(**

*τ***) is given by**

*x***is the deformation gradient tensor and in the present study, it is experimentally measurable. Further, the Cauchy stress**

*F***is related to the nominal stress**

*σ***by**

*τ**W*will depend on position

**only through the deformation gradient tensor**

*x***, i.e.,**

*F**W*(

**) =**

*x**W*(

**(**

*F***)).**

*x**μ*> 0,

*b*> 0, and

*n*> 0 are constants, and

*λ*

_{i}(

*i*= 1, 2, 3) are the local principal stretches, where

*λ*

_{1}⩾

*λ*

_{2}⩾

*λ*

_{3}.

*σ*

_{αβ}, according to the GNH model, are given by

*σ*

_{1}and

*σ*

_{2}are given by

*λ*

_{1}and

*λ*

_{2}are given, will be compared with the stresses determined using the technique and scheme proposed in the present study.

Note that to use Eq. (53) to predict stresses, we need to know the constants *μ*, *b*, and *n* for the PVC sheet material. Uniaxial tension test was performed using the standard dog-bone shaped samples made of the PVC sheet to calibrate the GNH model. Figure 20(a) presents the displacement fields over the tensile specimen surface during the tension process. The axial stretch *λ* is measured from the displacement field, while the nominal axial stress Σ is calculated by the applied tensile load divided by the original cross-sectional area of the tensile specimen. The variation of nominal axial stress Σ as function of the axial stretch *λ* from the uniaxial tension tests is shown in Fig. 20(b) as solid lines.

*λ*through

*μ*,

*b*, and

*n*for the PVC sheet. The fitting result is shown Fig. 20(b) as open diamonds. We see that the GNH model, with appropriate parameters determined, captures the deformation response of the PVC sheet quite well in the uniaxial tension case.

With the calibrated parameters *μ*, *b*, and *n* for the PVC sheet, we can use Eq. (53) to predict the principal stresses in the case of PVC sheet with a circular hole and subject to uniaxial tension, since the principal stretch fields *λ*_{1} and *λ*_{2} have been experimentally measured, as shown in Fig. 15. The predicted stress fields *σ*_{1} and *σ*_{2}, based on the GNH model, are presented in Fig. 21 for the selected moments indicated in Fig. 12(b).

*σ*

_{1}and

*σ*

_{2}; however in general, three components are required to describe the stress field in a planar deformation. By examining Eq. (52), we note that

**=**

*Q***(**

*R**θ*

_{q}) for planar deformation, and the property of the rotation tensor

**. Combining the above expression with Eq. (53), one will reach Eq. (21). Therefore, the principal stress orientation**

*R**θ*

_{σ}=

*θ*

_{λ}+

*θ*

_{q}indeed provides the third component of the stress field.

Corresponding to those stress fields predicted from the GNH model, shown in Fig. 21, we also determine the two principal stresses, *σ*_{1} and *σ*_{2}, using the technique proposed in the present study and based on the measurement shown in Fig. 16. These stress fields are shown in Fig. 22 for the selected moments indicated in Fig. 12(b).

We first note that the domain over which the stress fields are determined in Fig. 22 is different from that shown in Fig. 21. The domain, over which the stresses are given in Fig. 21, is the current or deformed domain. Figure 23(a) shows the sample domain in the undeformed or reference state, indicated as $Rref$. The boundaries or edges of domain $Rref$ are formed by some material particles. During the deformation, these boundaries or edges travel with the material particles and form the boundary of the current or deformed domain $R$, as shown in the left plot of Fig. 23(b). The third domain, over which the stresses are computed in Fig. 22, is referred to as the computational domain, shown in the right plot of Fig. 23(b), and according to the numerical scheme of present study, the boundary of the computational domain is composed of characteristic curves. Some of the domain boundaries are part of the characteristic curves, like the lateral edges and the boundary of the cavity. This part of boundary is the same in both Figs. 21 and 22. However, the top and the bottom edges of the current domain $R$ and the top and the bottom boundaries of the computational domain are different because the top and the bottom boundaries of the computational domain depend on the principal stress orientation *θ*_{σ} at any given moment. Since the principal stress orientation *θ*_{σ} measurement is over the current or deformed domain, the computational domain shown in Fig. 22 is a subdomain of the current domain $R$. This observation highlights the necessity of careful design of both the specimen and the loading fixture when doing the experiment, so that the computational domain can be as close to the current domain as possible, thus minimizes the degree of approximation of the traction boundary conditions.

The computed stress field in Fig. 22 and the predicted stress field in Fig. 21 have similar structures and slightly different overall magnitude, especially for the component of the major principal stress *σ*_{1}, where the predicted stress from the GNH model is higher than that of computed stress. Such a difference may be the result of the approximation of the traction boundary condition we just discussed above. The matching of the minor principal stress *σ*_{2} between the GNH model prediction and the computed stress is much better. On the other hand, we also see that the predicted stress from the GNH model does not satisfy the traction-free boundary conditions along the lateral edges and along the edge of the cavity. Along the lateral edges of the PVC sheet, *σ*_{2} is small but not zero as the boundary condition requires. Along the boundary of the cavity, *σ*_{2} obviously does not vanish. It is likely that the stress field predicted from the GNH model does not satisfy the equation of equilibrium as well.

Finally, we want to point out an interesting observation that in the computed stress field, stress concentration does not occur at the cavity edges, but slightly inward at location near the cavity boundary. In contrast, stress concentration occurs right at the cavity edge according to the stress field predicted from the GNH model, which is consistent with intuition and many closed-form solutions. The fact that stress concentration does not occur at the cavity edge in Fig. 22 might be due to the artifact of either the numerical scheme or the experiment itself. We note that the PVC sheet, or any polymer sheet, has very low bending stiffness and may be easily bent during a planar deformation causing out-of-plane motion, especially along edges of the specimen. Such out-of-plane motion, though may be exceedingly small, can influence the in-plane deformation measurement, and since the stress determination scheme proposed in this study relies on the detailed deformation field measurement, any error in deformation measurement will manifest itself in the resulting stress field.

### 7.4 Further Comments on the Traction Boundary Condition and Its Approximation.

In the experimental case of a PVC sheet with a circular hole and subject to uniaxial tension, discussed in this section, the imposed boundary conditions over the top and the bottom of the gage section are actually displacement boundary conditions. As demonstrated in Fig. 13, the lateral displacement *u*_{1} along the top and the bottom edge of the gage section is approximately zero and the vertical displacement *u*_{2} remains constant across the edges. Although experimental measurement has provided the precise description of the displacement boundary condition, but directly consider such a condition in the determination of the stresses would involve constitutive relations, which is against the premise of the present investigation.

To determine the stresses without invoking constitutive relations, we need to deal with the equation of equilibrium with traction boundary conditions, together with the additional information, the field of principal stress orientation *θ*_{σ}, provided by the experimental deformation measurement. The traction boundary condition is exact over some portion of the boundary, but approximate over others. Careful specimen and loading fixture design may reduce the effect of approximation, as we discussed in Sec. 7.3, but there is another avenue in getting better approximation of traction boundary conditions, and we will give a rough sketch of the approach here.

*β*-lines, $L\beta (1)$ and $L\beta (N)$, over which major principal stress

*σ*

_{1}needs to be described as additional boundary conditions. We formally write the boundary condition as

*f*

_{1}(

*s*

_{α},

*s*

_{β}) and

*f*

_{N}(

*s*

_{α},

*s*

_{β}) are yet to be determined, but they should satisfy the constraint

*W*is the original width of the PVC sheet sample, and $S\beta (1)$ and $S\beta (N)$ are the total lengths of the

*β*-lines, $L\beta (1)$ and $L\beta (N)$, respectively.

Furthermore, we assume that *f*_{1}(*s*_{α}, *s*_{β}) and *f*_{N}(*s*_{α}, *s*_{β}) can be represented as the sum of a series of known functions, *ξ*_{i}(*s*_{α}, *s*_{β}) (*i* = 1, 2, …, *K*), like

*i*= 1, 2, …,

*K*), following the scheme discussed in the present investigation, we may determine a stress field

*λ*

_{1}and

*λ*

_{2}, we may define an energetic measure $\Pi (ai(1),ai(N),i=1,2,\u2026,K)$. Now we may invoke some variational principles to choose the set of coefficients $ai(1)$ and $ai(N)$ (

*i*= 1, 2, …,

*K*) such that the energetic measure Π is minimized. Once the set of coefficients $ai(1)$ and $ai(N)$ (

*i*= 1, 2, …,

*K*) is obtained, the stresses within the PVC sheet is also determined.

One can see that the approach sketched here for obtaining more accurate traction boundary conditions along the characteristic lines, $L\beta (1)$ and $L\beta (N)$, is much more numerically involved. However, we have converted the stress determination problem into an optimization problem of minimizing the energetic measure Π, where the stresses are just the by-product of the process, as a result, more powerful numerical techniques, like machine learning and neural networking, can be applied.

This approach may go even further. We may directly consider the displacement boundary conditions as shown in Fig. 13. Rather than limited to the characteristic curve boundary of the computational domain, we will try to determine the tractions along the top and the bottom edge of the current domain shown in Fig. 23, where the displacement is measured experimentally. At any point along this displacement-prescribed boundary, we need to specify the value of a combination of both *σ*_{1} and *σ*_{2}, as indicated in Eq. (25), which can be expressed as a sum of the expansion of known functions with yet-to-be-determined coefficients. Through the optimization or minimization of the energetic measure Π, both the tractions along the displacement-prescribed boundary and the stresses within the current domain can be obtained.

## 8 Concluding Remarks

In this study, we explored and demonstrated a technique that will determine the stresses associated with a nonuniform deformation field without knowing the detailed constitutive behavior of the deforming material. Here, we do not attempt to deal with the complete boundary value problem in continuum mechanics, but rather we focus on some special circumstances, where if the detailed deformation field over a domain has been measured and traction boundary condition can be identified or approximated, the equation of equilibrium can be solved numerically for the two principal stresses. Since the technique relies on the full-field experimental measurement of the deformation field, in particular, the field of the principal stretch orientation over the whole domain, we are limited to planar or two-dimensional (2D) field only. Meanwhile, the technique also relies on the observation that the Lagrangian strain and the second Piola–Kirchhoff stress for isotropic materials share the same eigenvectors; therefore, the technique proposed in the present study applies to isotropic solids only.

The experimentally measured deformation field over the whole domain plays the key role in the proposed technique. The measured deformation field not only provides the principal stretch orientation *θ*_{λ} at every location but also provides the rigid-body rotation *θ*_{q} of the material particle at that location. The principal stress orientation *θ*_{σ} is the sum of these two angles, according to Eq. (22). With the principal stress orientation *θ*_{σ} field specified, the equation of equilibrium becomes numerically solvable with only the traction boundary conditions. Thus, the numerical computation of the stresses closely intertwines with the experimental measurement.

The close linkage of experimental measurement and numerical computation also poses high demand regarding the accuracy and signal/noise ratios of the experimental technique used for deformation field measurements. Any experimental measurement results contain uncertainty and noise. How sensitive of the proposed numerical scheme on the measurement uncertainty and noise has not been addressed in the present study. It is believed that the optical technique DIC used in the present investigation provides very high accuracy in measuring the displacement field, especially for finite or large deformations. The close match of the lateral edges of the PVC sheet, identified from the DIC result, and the characteristic curves, which are the results of solving the ODEs based on the measured *θ*_{σ} field, is a clear indication of such accuracy of the DIC technique.

The numerical scheme proposed and used in this investigation for solving the equation of equilibrium and computing the stresses depends on the characteristic curves. Once the principal stress orientation *θ*_{σ} field is specified, the two family of the characteristic curves are completely determined. The role of characteristic curves in solving first-order PDEs is quite important. For a single PDE involving a single unknown function, if the PDE can be recast onto the characteristic line, it becomes an ODE and is readily solvable. The numerical scheme proposed in the present study uses similar idea. The traction boundary conditions are described over the characteristic curves that form the boundary of the computational domain, as discussed in Sec. 7.3. Some of these traction boundary conditions are exact and some are approximate, and, in addition to a better sample and loading fixture design, we have given a rough sketch in Sec. 7.4 for improving the accuracy of the tractions over some of the characteristic curves. Such improvement will be explored in a future study.

With the deformation field measured from experiment and complemented with the stresses determined using the technique and the numerical scheme discussed in this study, we have a useful tool for examining different constitutive models on different materials.

The technique discussed in this paper could be extended in different directions and we briefly discuss a few of them here. The first is the situation where domain $R$ is composed of several non-overlapping subdomains, and each subdomain has a sharp boundary. Suppose that we have measured the deformation field and the principal stress orientation *θ*_{σ} field over $R$, so we should be able to determine the stresses within $R$. However, since the principal stress orientation *θ*_{σ} may suffer a jump across the subdomain boundary and the equation of equilibrium (28) involves gradient of the principal stress orientation *θ*_{σ}, the coefficients of Eq. (28) will have singularities along certain lines within $R$. The existence of a solution to Eq. (28) and how to compute the solution in the presence of singularities become an issue and should be investigated.

*σ*

_{1}and

*σ*

_{2}, and the principal stress orientation

*θ*

_{σ}are given by

*a*)

*b*)

*t*⩾ 0,

*ρ*is the mass density, and $u\xa81$ and $u\xa82$ are the particle acceleration components in the (

*y*

_{1},

*y*

_{2}) coordinate system. Again, from the left-hand side of Eq. (55), we see that there exist two family of characteristic curves governed by the two ODEs given by Eq. (27). By introducing the curvilinear coordinates (

*s*

_{α},

*s*

_{β}) like in the quasi-static case, the equation of motion (55) can be recast into

*s*

_{α},

*s*

_{β}). Equation (56) is the same as Eq. (28) except the additional terms on the right-hand side. From the experimental perspective, at any given moment during the deformation, the displacement fields

*u*

_{1}and

*u*

_{2}can be measured, so is the principal stress orientation field

*θ*

_{σ}. If a series of such measurement is carried out over time, the particle velocity fields $u\u02d91$ and $u\u02d92$ can also be determined. Some experimental packages have already provided such options. With one more careful differentiation with respect to time, the particle acceleration fields $u\xa81$ and $u\xa82$ can be computed. Therefore, the numerical scheme discussed in this investigation can be used to solve the equation of motion (56), together with appropriate traction boundary conditions, to determine the stresses associated with the dynamic deformation process.

In present investigation, we limited ourselves to 2D deformations for that presently the state of the art technology in experimental mechanics only allows robust and accurate full-field deformation measurement in 2D over an entire domain. However, intensive research activities are underway to develop technologies for measuring true 3D deformations within a body. One such technique is the digital volume correlation (DVC), e.g., Refs. [19–21], a volumetric extension of the digital image correlation (DIC) technique, which we used in the present study. Therefore, we may inspect the prospect of extending the technique of the present study to 3D deformations.

**(**

*u***) is measured for $x\u2208R$, as a result, deformation gradient tensor**

*x***can be computed. The three principal stretches**

*F**λ*

_{i}(

*i*= 1, 2, 3), where

*λ*

_{1}⩾

*λ*

_{2}⩾

*λ*

_{3}, together with their directions can be determined as well. Such three directional vectors will form an orthogonal frame and can be specified by the three Euler angles, say (

*φ*,

*θ*,

*ψ*). The right stretch tensor

**can be written as**

*U***is the rotation tensor and (**

*R**φ*

_{λ},

*θ*

_{λ},

*ψ*

_{λ}) are the Euler angles associated with the right stretch tensor

**. The proper orthogonal tensor**

*U***=**

*Q*

*FU*^{−1}representing a local rotation and, therefore, can be written as

*φ*

_{q},

*θ*

_{q},

*ψ*

_{q}) are the Euler angles associated with such a rotation. Note that from the true 3D deformation measurement, principal stretches

*λ*

_{i}, Euler angles (

*φ*

_{λ},

*θ*

_{λ},

*ψ*

_{λ}) and (

*φ*

_{q},

*θ*

_{q},

*ψ*

_{q}) over the 3D body are obtainable. Now the Cauchy stress

**has the representation**

*σ**σ*

_{i}(

*i*= 1, 2, 3) are the three principal stresses and (

*φ*

_{σ},

*θ*

_{σ},

*ψ*

_{σ}) are the Euler angles associated with the principal stress directions. For

*isotropic*materials, we have

**are experimentally determinable. The 3D equation of equilibrium will have only three unknowns,**

*σ**σ*

_{i}(

*i*= 1, 2, 3), and three PDEs. Thus, the technique discussed in the present investigation can be applied to the 3D deformations. However, just like the DVC technique, which is facing the challenge of computational cost, the stress determination technique discussed here will face the same challenge in computations when applied to 3D deformations.

## Acknowledgment

Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration (NNSA) of U.S. Department of Energy (Contract No.89233218CNA000001).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The authors attest that all data for this study are included in the paper.