Nonuniform Stress Field Determination Based on Deformation Measurement

Liu, C.

doi:10.1115/1.4050535

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, σ₁ and σ₂, 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.

Even for the very simplified situations, like the uniaxial stress, where we indeed can talk about stresses, the stresses are calculated rather than directly measured in the usual sense. In a uniaxial stress experiment, say uniaxial tension, the gage section of the specimen has the shape of a uniform cylinder. A tensile force is applied to the cylinder along its axis. The experimenter must ensure that the deformation within the gage section is homogeneous and usually the experimenter feels confident about this requirement being met because the deformation can be directly measured and observed. Next, the experimenter also wants to make sure that the traction applied over the top and the bottom surfaces of the gage section can be approximated as σ = 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

σ_{i j} = σ δ_{i 3} δ_{3 j}, i, j = 1, 2, 3

Here, x₃ 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

Mathematically, a deformation is just a mapping of a material particle from its original position x, in the undeformed state, to the new position y 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

R_{ref}

⁠. The motion of the 2D object can be written as

y_{α} = {\hat{y}}_{α} (x_{1}, x_{2}), x_{α} \in R_{ref}, α = 1, 2

(1)

Here, y_α is the position vector at the current moment during the motion of a particle that is located at x_α in the reference region

R_{ref}

⁠. We can also write

y_{α} = x_{α} + u_{α} (x_{1}, x_{2}), x_{α} \in R_{ref}, α = 1, 2

(2)

where u_α(x₁, x₂) 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

{\hat{y}}_{α} (x_{1}, x_{2})

is twice differentiable continuous over region

R_{ref}

⁠, the mapping

{\hat{y}}_{α}

is one-on-one, and the Jacobian

J (x_{1}, x_{2}) = det {{\hat{y}}_{α},_{β} (x_{1}, x_{2})} > 0

for all

x_{α} \in R_{ref}

(α, β = 1, 2).

The deformation

y_{α} = {\hat{y}}_{α} (x_{1}, x_{2})

defined above can be viewed as locally homogeneous, so that it can be written in the form:

y_{α} = F_{α β} (x_{1}, x_{2}) x_{β} + b_{α} (x_{1}, x_{2}), α, β = 1, 2

(3)

where b_α is a vector and F_αβ is a tensor dependent on location (x₁, x₂). Since

J (x_{1}, x_{2}) = det {{\hat{y}}_{α},_{β} (x_{1}, x_{2})} = det {F_{α β}}

⁠, necessarily

det {F_{α β}} > 0

⁠. Meanwhile, the deformation gradient F_αβ can be related to the displacement u_α through

F_{α β} = δ_{α β} + \frac{\partial u_{α}}{\partial x_{β}}, α, β = 1, 2

(4)

Consider a “material fiber” with one end located at position (x₁, x₂) 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₁, x₂) 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

λ (x_{1}, x_{2}, θ) \equiv \frac{δ L}{δ ℓ} = (F_{α β} F_{η α} m_{β} m_{η})^{1 / 2}

(5)

where

[m_{α}] = (\cos θ, \sin θ)^{⊤}

⁠. So, the stretch λ is a function of location (x₁, x₂) and angle θ. The extension of a material segment in the direction of m_α, is simply λ − 1.

If we fix the location (x₁, x₂) 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

\frac{d λ}{d θ} |_{θ = θ_{λ}} = 0

(6)

The stretches along the direction with the directional angle θ_λ are the principal stretches. In planar deformation, the direction perpendicular to (x₁, x₂) plane is always a principal direction. The other two principal stretches λ₁ and λ₂, where λ₁ ⩾ λ₂, are within the (x₁, x₂)-plane, and we assign the θ_λ to be measured from the x₁ axis to the major stretch λ₁ counterclockwise.

From experimental perspective, all the quantities discussed so far can be measured in an experiment over the sample domain $R_{ref}$ ⁠. For example, using the optical technique of digital image correlation (DIC), one can measure the displacement field u_α over the region $R_{ref}$ ⁠. 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₁, x₂) and for any arbitrary angle θ from Eq. (5), and this will lead to the determination of the principal stretches λ₁ and λ₂, and the principal stretch orientation θ_λ at any given location (x₁, x₂).

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.

We start with a general homogeneous deformation, which can be written as

y = F x + b, x \in R_{ref}, y \in R

(7)

where

det F > 0

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

F = Q U = V Q

(8)

where the right-stretch tensor

U^{⊤} = U

and the left-stretch tensor

V^{⊤} = V

also

det U > 0

and

det V > 0

⁠, and Q is a proper orthogonal tensor, i.e.,

Q^{- 1} = Q^{⊤}

and

det Q = 1

⁠. Finally,

V = Q U Q^{⊤}

⁠.

In Sec. 2, we have stated that for a two-dimensional domain

R_{ref}

⁠, the displacement field u_α(x₁, x₂) (α = 1, 2) can be experimentally measured, so is the deformation gradient tensor F_αβ. Meanwhile, the principal stretches, λ₁ and λ₂, and the principal stretch orientation θ_λ, can be determined as well. Since the principal stretches, λ₁ and λ₂, are the two eigen values of the right-stretch tensor U, the two eigenvectors associated with λ₁ and λ₂ are

N_{1} = (\cos θ_{λ}, \sin θ_{λ})^{⊤}, N_{2} = (- \sin θ_{λ}, \cos θ_{λ})^{⊤}

(9)

In the Cartesian coordinate system (x₁, x₂), the right-stretch tensor U can be written as

U = R (θ_{λ}) Λ R^{⊤} (θ_{λ})

(10)

where

R (θ_{λ}) = (\begin{matrix} \cos θ_{λ} & - \sin θ_{λ} \\ \sin θ_{λ} & \cos θ_{λ} \end{matrix}), Λ = (\begin{matrix} λ_{1} & 0 \\ 0 & λ_{2} \end{matrix})

One can show that

R^{- 1} (θ_{λ}) = R^{⊤} (θ_{λ})

⁠, so that R(θ_λ) is a proper orthogonal tensor, as well as a rotational tensor associated with angle θ_λ. Meanwhile, the inverse of the right-stretch tensor U can be written through the principal stretches, λ₁ and λ₂, and the principal stretch orientation θ_λ as

U^{- 1} = R^{- ⊤} (θ_{λ}) Λ^{- 1} R^{- 1} (θ_{λ}) = R (θ_{λ}) Λ^{- 1} R^{⊤} (θ_{λ})

which can be expressed in the matrix form

U^{- 1} = (\begin{matrix} \frac{\cos^{2} θ_{λ}}{λ_{1}} + \frac{\sin^{2} θ_{λ}}{λ_{2}} & (\frac{1}{λ_{1}} - \frac{1}{λ_{2}}) \sin θ_{λ} \cos θ_{λ} \\ (\frac{1}{λ_{1}} - \frac{1}{λ_{2}}) \sin θ_{λ} \cos θ_{λ} & \frac{\sin^{2} θ_{λ}}{λ_{1}} + \frac{\cos^{2} θ_{λ}}{λ_{2}} \end{matrix})

(11)

As a result, the orthogonal tensor Q can be written as

Q = F U^{- 1}

(12)

For planar or 2D deformation, Q represents a rigid-body rotation, and we may express the orthogonal tensor Q as a rotation tensor R(θ_q), where the rotation angle θ_q is completely determined once the deformation gradient F_αβ, the principal stretches λ₁ and λ₂, and the principal stretch orientation θ_λ are experimentally measured.

The Lagrangian strain (or Green-Lagrangian strain), E, is defined by

E = \frac{1}{2} (U^{2} - I) = \frac{1}{2} (R (θ_{λ}) Λ^{2} R^{⊤} (θ_{λ}) - I)

or in matrix form

E = \frac{1}{2} (\begin{matrix} λ_{1}^{2} \cos^{2} θ_{λ} + λ_{2}^{2} \sin^{2} θ_{λ} - 1 & (λ_{1}^{2} - λ_{2}^{2}) \sin θ_{λ} \cos θ_{λ} \\ (λ_{1}^{2} - λ_{2}^{2}) \sin θ_{λ} \cos θ_{λ} & λ_{1}^{2} \sin^{2} θ_{λ} + λ_{2}^{2} \cos^{2} θ_{λ} - 1 \end{matrix})

(13)

One can show that the Lagrangian strain E has two eigen values,

(λ_{1}^{2} - 1) / 2

and

(λ_{2}^{2} - 1) / 2

⁠, with the two eigenvectors N₁ and N₂ given in Eq. . Therefore, the principal stretch orientation θ_λ also defines the direction of the major principal strain,

(λ_{1}^{2} - 1) / 2

⁠, of the Lagrangian strain E.

To summarize, from experiment, we can measure the displacement field u_α(x₁, x₂) (α = 1, 2), the deformation gradient tensor F, the principal stretches, λ₁ and λ₂, and the principal stretch orientation θ_λ, and consequently, we can also determine the components involved in the polar decomposition of the deformation gradient tensor F, U, and Q.

3.2 Notion of Tensor-Valued Isotropic Functions.

Let A and B be two second-order and symmetric tensors, and let φ be a tensor-valued function such that

B = φ (A)

(14)

The function φ is called isotropic if for every proper orthogonal tensor Q [9],

Q B Q^{⊤} = φ (Q A Q^{⊤})

(15)

Since the proper orthogonal tensor Q represents a rotation, the tensor-valued isotropic function φ is invariant under arbitrary rotations. One consequence of an isotropic tensor-valued function φ that relates A and B is that tensors A and B are co-axial, or they share the same (orthonormal) eigenvectors. Also note that tensors A and B are co-axial if and only if 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.

The stress that is the work conjugate of the Lagrangian strain E is the second Piola–Kirchhoff stress S, where

S^{⊤} = S

⁠. Suppose that stress S has the principal components s₁ and s₂, where s₁ ⩾ s₂, and the major principal stress s₁ is in the direction described by the angle θ_s measured from x₁-axis to s₁ counterclockwise. The components of S in the Cartesian coordinate (x₁, x₂) can be related to the principal stresses s₁ and s₂ by

S = R (θ_{s}) \hat{S} R^{⊤} (θ_{s}), \hat{S} = (\begin{matrix} s_{1} & 0 \\ 0 & s_{2} \end{matrix})

(16)

where R(θ_s) is the rotation tensor associated with angle θ_s.

From the notion of isotropic tensor-valued function, discussed in Sec. 3.2, for isotropic materials, the principal stress orientation θ_s of the second Piola–Kirchhoff stress S, is the same as the principal stretch orientation θ_λ, or the principal orientation of the Lagrangian strain E. Since both S and E are expressed in the reference or undeformed configuration,

θ_{s} (x_{1}, x_{2}) = θ_{λ} (x_{1}, x_{2}), (x_{1}, x_{2}) \in R_{ref}

(17)

This is just the 2D specialization of the co-axiality of the Lagrangian strain E and the second Piola–Kirchhoff stress S for isotropic materials. As we have mentioned in Sec. 3.1 that the principal stretch orientation θ_λ(x₁, x₂) can be measured experimentally over the reference domain

R_{ref}

⁠, as a result, the principal stress orientation θ_s(x₁, x₂) of the second Piola–Kirchhoff stress S is also experimentally measurable over domain

R_{ref}

⁠.

3.4 Cauchy Stress and Principal Stress Orientation.

The stress associated with the deformed configuration is the Cauchy stress σ. The Cauchy stress σ is related to the second Piola–Kirchhoff stress S, through the deformation gradient tensor F, by

σ = \frac{1}{J} F S F^{⊤}, J = det F

(18)

Using the expression in Eq. , we have

\begin{aligned} σ & = \frac{1}{J} F R (θ_{s}) \hat{S} R^{⊤} (θ_{s}) F^{⊤} \\ = \frac{1}{J} R (θ_{q}) U R (θ_{λ}) \hat{S} R^{⊤} (θ_{λ}) U R^{⊤} (θ_{q}) \end{aligned}

(19)

Here, we have used the following facts, the co-axiality of the Lagrangian strain E and the second Piola–Kirchhoff stress S, i.e., θ_s = θ_λ, the polar decomposition F = QU, and Q = R(θ_q). Furthermore, using Eq. ,

\begin{aligned} σ & = R (θ_{q}) R (θ_{λ}) (\frac{1}{J} Λ \hat{S} Λ) R^{⊤} (θ_{λ}) R^{⊤} (θ_{q}) \\ = R (θ_{q} + θ_{λ}) (\frac{1}{J} Λ \hat{S} Λ) R^{⊤} (θ_{q} + θ_{λ}) \end{aligned}

(20)

The property of the rotation tensor in 2D, where R(θ_q)R(θ_λ) = R(θ_q + θ_λ), has been used in Eq. . Note that

Λ \hat{S} Λ / J

is a diagonal tensor, so Eq. provides one form of diagonalization of the Cauchy stress tensor σ for a planar or 2D deformation.

On the other hand, suppose that in the current or deformed state, the Cauchy stress σ has the principal components σ₁ and σ₂, where σ₁ ⩾ σ₂, with the principal stress orientation θ_σ, which is measured from the y₁-axis to σ₁ counterclockwise, we will have

σ = R (θ_{σ}) Σ R^{⊤} (θ_{σ}), Σ = (\begin{matrix} σ_{1} & 0 \\ 0 & σ_{2} \end{matrix})

(21)

By comparing Eqs. and , we conclude that

θ_{σ} = θ_{q} + θ_{λ}

or more precisely, the principal stress orientation θ_σ 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 Q, i.e.,

θ_{σ} (y_{1}, y_{2}) = θ_{q} (x_{1}, x_{2}) + θ_{λ} (x_{1}, x_{2}), x_{α} \in R_{ref}

(22)

where

y_{α} = {\hat{y}}_{α} (x_{1}, x_{2})

(α = 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. , 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.

In the absence of body force density, equation of equilibrium of a 2D deformation in a Cartesian coordinate system (y₁, y₂), together with the traction boundary condition, is given by

\begin{aligned} σ_{α β},_{β} & = 0, (y_{1}, y_{2}) \in R \\ σ_{α β} n_{β} & = t_{α}, (y_{1}, y_{2}) \in \partial R \end{aligned}

(23)

Here, we consider the Cauchy stress σ, where

σ^{⊤} = σ

⁠. Domain

R

represents the deforming region in the current configuration and its boundary

\partial R

has the unit outward normal vector n.

We consider the principal components of the Cauchy stress σ in the 2D region

R

⁠, σ₁ and σ₂, and with the associated principal stress orientation θ_σ, measured from y₁-axis to the major principal stress σ₁ counterclockwise. In the Cartesian coordinate (y₁, y₂), the Cauchy stress σ can be expressed in terms of its principal components σ₁ and σ₂, and the principal stress orientation θ_σ by Eq. , or to be more explicit,

\begin{aligned} σ_{11} & = σ_{1} \cos^{2} θ_{σ} + σ_{2} \sin^{2} θ_{σ} \\ σ_{12} & = σ_{21} = (σ_{1} - σ_{2}) \sin θ_{σ} \cos θ_{σ} \\ σ_{22} & = σ_{1} \sin^{2} θ_{σ} + σ_{2} \cos^{2} θ_{σ} \end{aligned}

Substitute the stresses into Eq. , the transformed version of the equation of equilibrium is given by

\begin{aligned} M_{1} (θ_{σ}) \frac{\partial σ}{\partial y_{1}} + M_{2} (θ_{σ}) \frac{\partial σ}{\partial y_{2}} + (σ_{1} - σ_{2}) M_{3} (θ_{σ}) \nabla θ_{σ} = 0, \\ (y_{1}, y_{2}) \in R \end{aligned}

(24)

with the traction boundary condition

Q (θ_{σ}, n) σ = t, (y_{1}, y_{2}) \in \partial R

(25)

where

σ = (σ_{1}, σ_{2})^{⊤}

and

t = (t_{1}, t_{2})^{⊤}

⁠, and

\begin{aligned} M_{1} & = \frac{1}{2} (\begin{matrix} 2 \cos^{2} θ_{σ} & 2 \sin^{2} θ_{σ} \\ \sin 2 θ_{σ} & - \sin 2 θ_{σ} \end{matrix}) \\ M_{2} & = \frac{1}{2} (\begin{matrix} \sin 2 θ_{σ} & - \sin 2 θ_{σ} \\ 2 \sin^{2} θ_{σ} & 2 \cos^{2} θ_{σ} \end{matrix}) \\ M_{3} & = (\begin{matrix} - \sin 2 θ_{σ} & - \cos 2 θ_{σ} \\ \cos 2 θ_{σ} & \sin 2 θ_{σ} \end{matrix}) \end{aligned}

and

\begin{aligned} Q = (\begin{matrix} \cos θ_{σ} m_{1} (θ_{σ}) \cdot n & - \sin θ_{σ} m_{2} (θ_{σ}) \cdot n \\ \sin θ_{σ} m_{1} (θ_{σ}) \cdot n & \cos θ_{σ} m_{2} (θ_{σ}) \cdot n \end{matrix}) \\ m_{1} (θ_{σ}) = (\begin{matrix} \cos θ_{σ} \\ \sin θ_{σ} \end{matrix}), m_{2} (θ_{σ}) = (\begin{matrix} - \sin θ_{σ} \\ \cos θ_{σ} \end{matrix}) \end{aligned}

Note that we can decouple the derivative terms of the principal stresses or diagonalize the first two matrices in Eq. simultaneously, so that the equation of equilibrium becomes

\begin{aligned} \frac{\partial σ_{1}}{\partial y_{1}} + \tan θ_{σ} \frac{\partial σ_{1}}{\partial y_{2}} = (\tan θ_{σ} \frac{\partial θ_{σ}}{\partial y_{1}} - \frac{\partial θ_{σ}}{\partial y_{2}}) (σ_{1} - σ_{2}) \\ \frac{\partial σ_{2}}{\partial y_{1}} - \cot θ_{σ} \frac{\partial σ_{2}}{\partial y_{2}} = (\cot θ_{σ} \frac{\partial θ_{σ}}{\partial y_{1}} + \frac{\partial θ_{σ}}{\partial y_{2}}) (σ_{1} - σ_{2}) \end{aligned}

(26)

where

(y_{1}, y_{2}) \in R

⁠. In Sec. 3.4, we have established that for isotropic materials, the principal stress orientation θ_σ can be related to the principal stretch orientation θ_λ through Eq. , 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₁, y₂) imposes a constraint to the PDEs in Eq. (26), where we have two unknown functions, the principal stresses σ₁ and σ₂, and two PDEs. With appropriate traction boundary conditions, the stress field, (σ₁, σ₂), 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.

The transformed equation of equilibrium is a set of first-order and linear PDEs. In its present form, we observe that it has two-family of characteristic curves, or lines, governed by two ordinary differential equations (ODEs),

\begin{aligned} σ_{1} : \frac{d y_{2}}{d y_{1}} = \tan θ_{σ} (y_{1}, y_{2}), α -lines \\ σ_{2} : \frac{d y_{2}}{d y_{1}} = - \cot θ_{σ} (y_{1}, y_{2}), β -lines \end{aligned}

(27)

where

(y_{1}, y_{2}) \in R

⁠. Here, we designate the characteristic curves associated with the major principal stress σ₁ as the α-lines, and the ones associated with the minor principal stress σ₂ as the β-lines. Once the field θ_σ(y₁, y₂) 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₂/dy₁, 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₁, y₂). Unlike the Cartesian coordinates (y₁, y₂), where there is a well defined origin y₁ = y₂ = 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₁, y₂), the curvilinear coordinate (s_α, s_β) should form a right-handed system.

To recast the equation of equilibrium into the curvilinear coordinate system (s_α, s_β), consider the derivative of an arbitrary variable f along an α-line and along a β-line,

\begin{aligned} \frac{\partial f}{\partial s_{α}} & = \cos θ_{σ} \frac{\partial f}{\partial y_{1}} + \sin θ_{σ} \frac{\partial f}{\partial y_{2}} \\ \frac{\partial f}{\partial s_{β}} & = - \sin θ_{σ} \frac{\partial f}{\partial y_{1}} + \cos θ_{σ} \frac{\partial f}{\partial y_{2}} \end{aligned}

Now

\begin{aligned} \frac{\partial σ_{1}}{\partial s_{α}} & = \cos θ_{σ} \frac{\partial σ_{1}}{\partial y_{1}} + \sin θ_{σ} \frac{\partial σ_{1}}{\partial y_{2}} = - (σ_{1} - σ_{2}) \frac{\partial θ_{σ}}{\partial s_{β}} \\ \frac{\partial σ_{2}}{\partial s_{β}} & = - \sin θ_{σ} \frac{\partial σ_{2}}{\partial y_{1}} + \cos θ_{σ} \frac{\partial σ_{2}}{\partial y_{2}} = - (σ_{1} - σ_{2}) \frac{\partial θ_{σ}}{\partial s_{α}} \end{aligned}

Denote that

\frac{\partial θ_{σ}}{\partial s_{α}} = θ_{σ},_{s_{α}} (s_{α}, s_{β}), \frac{\partial θ_{σ}}{\partial s_{β}} = θ_{σ},_{s_{β}} (s_{α}, s_{β})

where

θ_{σ},_{s_{α}} (s_{α}, s_{β})

⁠, and

θ_{σ},_{s_{β}} (s_{α}, s_{β})

are two known functions of the coordinate (s_α, s_β), the recast PDEs in Eq. become

\begin{aligned} \frac{\partial σ_{1}}{\partial s_{α}} + (σ_{1} - σ_{2}) θ_{σ},_{s_{β}} (s_{α}, s_{β}) & = 0 \\ \frac{\partial σ_{2}}{\partial s_{β}} + (σ_{1} - σ_{2}) θ_{σ},_{s_{α}} (s_{α}, s_{β}) & = 0 \end{aligned}

(28)

For later convenience, we also write the principal stress orientation field as θ_σ(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, σ₁ and σ₂ once the functions $θ_{σ},_{s_{α}} (s_{α}, s_{β})$ and $θ_{σ},_{s_{β}} (s_{α}, s_{β})$ are given, will be based on Eq. (28).

5 Numerical Scheme

The general strategy of numerically solving the PDEs in Eq. is based on the following observation. Take the first PDE in Eq. as an example, suppose that the principal stress σ₂ is known as a function of the curvilinear coordinate (s_α, s_β) along a given α-line, say L_α, the PDE can be written as

\frac{\partial σ_{1}}{\partial s_{α}} + θ_{σ},_{s_{β}} (s_{α}, s_{β}) σ_{1} = θ_{σ},_{s_{β}} (s_{α}, s_{β}) σ_{2} (s_{α}, s_{β})

(29)

where (s_α, s_β) ∈ L_α. Note that the first-order linear PDE in the above equation only involves a single-derivative term, ∂σ₁/∂s_α, and thus, it is integrable along the characteristic α-line and one can formally write the expression for σ₁ as

\begin{aligned} σ_{1} (s_{α}, s_{β}) & = c_{α} (s_{β}) \exp (- \int^{s_{α}} θ_{σ},_{s_{β}} (ξ, s_{β}) d ξ) \\ + \exp (- \int^{s_{α}} θ_{σ},_{s_{β}} (ξ, s_{β}) d ξ) \\ \times \int^{s_{α}} {θ_{σ},_{s_{β}} (η, s_{β}) σ_{2} (η, s_{β}) \exp (- \int^{η} θ_{σ},_{s_{β}} (ξ, s_{β}) d ξ)} d η \end{aligned}

(30)

where 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_α.

Similarly, suppose that the principal stress σ₁ is known as a function of the curvilinear coordinate (s_α, s_β) in the second PDE in Eq. along a given β-line, L_β, we will have

\frac{\partial σ_{2}}{\partial s_{β}} - θ_{σ},_{s_{α}} (s_{α}, s_{β}) σ_{2} = - θ_{σ},_{s_{α}} (s_{α}, s_{β}) σ_{1} (s_{α}, s_{β})

(31)

where (s_α, s_β) ∈ L_β, and the formal expression for σ₂ can be written as

\begin{aligned} σ_{2} (s_{α}, s_{β}) & = c_{β} (s_{α}) \exp (\int^{s_{β}} θ_{σ},_{s_{α}} (s_{α}, ξ) d ξ) \\ - \exp (\int^{s_{β}} θ_{σ},_{s_{α}} (s_{α}, ξ) d ξ) \\ \times \int^{s_{β}} {θ_{σ},_{s_{α}} (s_{α}, η) σ_{1} (s_{α}, η) \exp (- \int^{η} θ_{σ},_{s_{α}} (s_{α}, ξ) d ξ)} d η \end{aligned}

(32)

where 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, σ₁ and σ₂, 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, σ₁ and σ₂, at discrete locations inside the deformed domain $R$ ⁠.

Suppose that for the given principal stress orientation field θ_σ(y₁, y₂) and according to Eq. , we already determined the two family of characteristic curves

L_{α}^{(i)}, L_{β}^{(j)}, i = 1, 2, \dots, M; j = 1, 2, \dots, N

Along a given α-line

L_{α}^{(i)}

⁠, the coordinate is described by (s_α, s_β) and in particular

0 ⩽ s_{α} ⩽ S_{α}^{(i)}

⁠, where

S_{α}^{(i)}

is the total length of the α-line

L_{α}^{(i)}

⁠. Similarly, coordinate s_β along a β-line

L_{β}^{(j)}

is such that

0 ⩽ s_{β} ⩽ S_{β}^{(j)}

⁠, where

S_{β}^{(j)}

is the total length of the β-line

L_{β}^{(j)}

⁠. If characteristic curves

L_{α}^{(i)}

and

L_{β}^{(j)}

intersect, the intersecting point is described by the coordinates

(s_{α}^{(i, j)}, s_{β}^{(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 and to obtain numerical values of the principal stresses, σ₁ and σ₂, at every intersecting point

(s_{α}^{(i, j)}, s_{β}^{(i, j)})

⁠. Note that in some domain configurations, e.g., multiply-connected domain, for a given pair (i, j), the characteristic lines

L_{α}^{(i)}

and

L_{β}^{(j)}

may not intersect, and as a result, we will have less than M × N intersecting points.

We start from the first characteristic line, say

L_{α}^{(1)}

in the α-line family, and assume that

L_{α}^{(1)}

is part of the domain boundary. We also assume that the traction boundary condition is such that the principal stress σ₂ is given over

L_{α}^{(1)}

⁠, denoted as

σ_{2} = σ_{2}^{(1)} (s_{α}, s_{β})

⁠. Note that the traction boundary condition over

L_{α}^{(1)}

cannot be described in terms of the principal stress σ₁. This is related to the restrictions of the so-called Cauchy problem and Cauchy data. Now, we have the equation for σ₁ over the one-dimensional line

L_{α}^{(1)}

(from the first PDE in Eq. ),

\begin{aligned} \frac{\partial σ_{1}}{\partial s_{α}} + θ_{σ},_{s_{β}} (s_{α}, s_{β}) σ_{1} = θ_{σ},_{s_{β}} (s_{α}, s_{β}) σ_{2}^{(1)} (s_{α}, s_{β}) \\ (s_{α}, s_{β}) \in L_{α}^{(1)} \end{aligned}

(33)

At the ends of the α-line

L_{α}^{(1)}

⁠, i.e., s_α = 0 or

s_{α} = S_{α}^{(1)}

⁠, the α-line

L_{α}^{(1)}

will intersect β-lines,

L_{β}^{(1)}

and

L_{β}^{(N)}

⁠. We assert that both

L_{β}^{(1)}

and

L_{β}^{(N)}

also form part of the boundary of domain

R

⁠, where the principal stress σ₁ is prescribed. In other words, the above equation has the boundary condition at both ends of characteristic line

L_{α}^{(1)}

for σ₁; thus, it can be numerically solved to have the value of σ₁ at discrete locations

(s_{α}^{(1, j)}, s_{β}^{(1, j)})

for j = 1, 2, …, N. The value of

σ_{2} = σ_{2}^{(1)} (s_{α}, s_{β})

at those locations is already known since it is part of the boundary condition.

Now, suppose that along the ith α-line,

L_{α}^{(i)}

⁠, the principal stresses, σ₁ and σ₂, have been obtained, and they are given at discrete locations

(s_{α}^{(i, j)}, s_{β}^{(i, j)})

as

(σ_{1}^{(i, j)}, σ_{2}^{(i, j)}), j = 1, 2, \dots, N

The second PDE in Eq. indicates that along the (i + 1)th α-line,

L_{α}^{(i + 1)}

⁠, the value for principal stress σ₂ can be estimated or updated to be

\begin{aligned} σ_{2}^{(i + 1, j)} & = σ_{2}^{(i, j)} - (σ_{1}^{(i, j)} - σ_{2}^{(i, j)}) θ_{σ},_{s_{α}} (s_{α}^{(i, j)}, s_{β}^{(i, j)}) δ s_{β}^{(i, j)}, \\ j = 1, 2, \dots, N \end{aligned}

(34)

at discrete locations

(s_{α}^{(i + 1, j)}, s_{β}^{(i + 1, j)})

⁠. In the above expression, every quantity on the right-hand side is known. In particular,

δ s_{β}^{(i, j)}

is the length of the line segment

L_{β}^{(j)}

between location

(s_{α}^{(i, j)}, s_{β}^{(i, j)})

on line

L_{α}^{(i)}

and location

(s_{α}^{(i + 1, j)}, s_{β}^{(i + 1, j)})

on line

L_{α}^{(i + 1)}

⁠. The line segment

δ s_{β}^{(i, j)}

is in general curved, only when lines

L_{α}^{(i)}

and

L_{α}^{(i + 1)}

are very close,

δ s_{β}^{(i, j)}

can be approximated by the direct distance between locations

(s_{α}^{(i, j)}, s_{β}^{(i, j)})

and

(s_{α}^{(i + 1, j)}, s_{β}^{(i + 1, j)})

⁠. As a result, the value of the principal stress σ₂ over line

L_{α}^{(i + 1)}

⁠, denoted as

σ_{2} = σ_{2}^{(i + 1)} (s_{α}, s_{β})

⁠, is known and the principal stress σ₁ satisfies the following equation over line

L_{α}^{(i + 1)}

⁠,

\begin{aligned} \frac{\partial σ_{1}}{\partial s_{α}} + θ_{σ},_{s_{β}} (s_{α}, s_{β}) σ_{1} & = θ_{σ},_{s_{β}} (s_{α}, s_{β}) σ_{2}^{(i + 1)} (s_{α}, s_{β}), \\ (s_{α}, s_{β}) \in L_{α}^{(i + 1)} \end{aligned}

(35)

and can be solved, together with the boundary condition at either ends of the line

L_{α}^{(i + 1)}

⁠, at discrete locations

(s_{α}^{(i + 1, j)}, s_{β}^{(i + 1, j)})

for j = 1, 2, …, N.

Such a process is thus repeated for all characteristic lines $L_{α}^{(i)}$ ⁠, i = 1, 2, …, M. The principal stresses, σ₁ and σ₂, are therefore obtained at all locations $(s_{α}^{(i, j)}, s_{β}^{(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_{β}^{(1)}$ ⁠, where the value of the principal stress σ₁ is known due to the fact that $L_{β}^{(1)}$ is part of the boundary of the domain $R$ and the traction boundary condition for σ₁ is given. Then, the principal stress σ₂ can be solved by the equation over the line $L_{β}^{(1)}$ ⁠. Using the similar updating scheme and repeating for all the β-lines $L_{β}^{(j)}$ (j = 1, 2, …, N), the principal stresses, σ₁ and σ₂, are therefore obtained at all locations $(s_{α}^{(i, j)}, s_{β}^{(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.

The configuration of the Brazilian disk compression is that a circular disk, made of homogeneous, isotropic, and linearly elastic material, and with the radius R, subjects to a pair of concentrated forces F₂ 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]:

\begin{aligned} σ_{11} & = \frac{F_{2}}{π R} {\frac{2 x_{1}^{2} (x_{2} + R) R}{(x_{1}^{2} + (x_{2} + R)^{2})^{2}} - \frac{2 x_{1}^{2} (x_{2} - R) R}{(x_{1}^{2} + (x_{2} - R)^{2})^{2}} - 1} \\ σ_{22} & = \frac{F_{2}}{π R} {\frac{2 (x_{2} + R)^{3} R}{(x_{1}^{2} + (x_{2} + R)^{2})^{2}} - \frac{2 (x_{2} - R)^{3} R}{(x_{1}^{2} + (x_{2} - R)^{2})^{2}} - 1} \\ σ_{12} & = \frac{2 F_{2}}{π R} {\frac{x_{1} (x_{2} + R)^{2} R}{(x_{1}^{2} + (x_{2} + R)^{2})^{2}} - \frac{x_{1} (x_{2} - R)^{2} R}{(x_{1}^{2} + (x_{2} - R)^{2})^{2}}} \end{aligned}

(36)

Also, for infinitesimal deformation, we do not need to distinguish the reference and the deformed (current) configurations.

Fig. 1

View large Download slide

(a) Configuration of the Brazilian disk compression test and (b) contour plot of the principal stress orientation θ_σ and associated characteristic curves

For the Brazilian disk compression, the principal stress orientation θ_σ can be expressed as

θ_{σ} (x_{1}, x_{2}) = \tan^{- 1} \frac{2 x_{1} x_{2}}{R^{2} + x_{1}^{2} - x_{2}^{2}}

(37)

The contour plot of θ_σ(x₁, x₂) 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. , which can be specialized for the Brazilian disk so that the α-lines are given by

\begin{aligned} x_{1}^{2} + {(x_{2} - \frac{R^{2} + a^{2}}{2 a})}^{2} = {(\frac{R^{2} - a^{2}}{2 a})}^{2}, \\ for - R ⩽ a ⩽ R and x_{1}^{2} + x_{2}^{2} ⩽ R^{2} \end{aligned}

(38)

and the β-lines are given by

\begin{aligned} {(x_{1} + \frac{R^{2} - b^{2}}{2 b})}^{2} + x_{2}^{2} = {(\frac{R^{2} + b^{2}}{2 b})}^{2}, \\ for - R ⩽ b ⩽ R and x_{1}^{2} + x_{2}^{2} ⩽ R^{2} \end{aligned}

(39)

Here, a and b are two parameters, where a is the vertical coordinate of the intersecting point of an α-line with the x₂-axis and b is the horizontal coordinate of the intersecting point of a β-line with the x₁-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 σ₁ = 0. As shown in both plots in Fig. 1(b), we should have the boundary condition σ₁ = 0 everywhere along the circular edge of the disk, except at the two poles and at the two poles, the minor principal stress σ₂ = σ = −F₂/R.

The gradient of the principal stress orientation θ_σ in the curvilinear coordinate (s_α, s_β), formed by the two family of characteristic curves, can be expressed as

\begin{aligned} \frac{\partial θ_{σ}}{\partial s_{α}} & = \frac{2 x_{2}}{\sqrt{x_{1}^{4} + (x_{2}^{2} - R^{2})^{2} + 2 x_{1}^{2} (x_{2}^{2} + R^{2})}} \\ \frac{\partial θ_{σ}}{\partial s_{β}} & = \frac{2 x_{1}}{\sqrt{x_{1}^{4} + (x_{2}^{2} - R^{2})^{2} + 2 x_{1}^{2} (x_{2}^{2} + R^{2})}} \end{aligned}

(40)

and the contour plots of the gradients are presented in Fig. 2. Now, all the required information for solving the equation of equilibrium are ready. Next, we apply the numerical scheme to obtain the field of the principal stresses σ₁ and σ₂.

Fig. 2

View large Download slide

Contour plot of the gradient of the principal stress orientation θ_σ in curvilinear coordinate (s_α, s_β)

To start the numerical solution, we consider the first α-line L_α(a₁), which is very close to the loading point near the bottom of the disk. As shown in Fig. 3(a), the minor principal stress σ₂ can be approximated by using the Boussinesq solution or

σ_{2}^{(1)} (s_{α}, s_{β} = r (θ)) = - \frac{2 F_{2}}{π} \cdot \frac{\sin θ}{r (θ)}

where r(θ) is the distance from the bottom pole to any location on L_α(a₁) and θ is the polar angle shown in Fig. 3(a). Once σ₂ is known over L_α(a₁), Eq. is used for solving σ₁ over L_α(a₁), with the boundary condition σ₁(s_α = 0, s_β) = 0. Equations and are used for solving the principal stresses, σ₁ and σ₂, over the subsequent α-lines L_α(a_i) (i = 1, 2, …, M). The update of the principal stress σ₂ from line L_α(a_i) to line L_α(a_i+1) is sketched in Fig. 3(b).

Fig. 3

(a) Boussinesq approximation of the minor principal stress σ2 along the first α-line Lα(a1), (b) update of the principal stress σ2 from line Lα(ai) to line Lα(ai+1), and (c) σ1 from line Lβ(bj) to line Lβ(bj+1)

View large Download slide

(a) Boussinesq approximation of the minor principal stress σ₂ along the first α-line L_α(a₁), (b) update of the principal stress σ₂ from line L_α(a_i) to line L_α(a_i+1), and (c) σ₁ from line L_β(b_j) to line L_β(b_j+1)

We can also switch the role of α- and β-lines and start from the first β-line L_β(b₁), where σ₁ ≡ 0 according to the traction boundary condition. The minor principal stress σ₂ over L_β(b₁) can be obtained by solving the integrable PDE with the end condition given by the Boussinesq solution. The update of the major principal stress σ₁ along α-lines L_α(a) will lead to the value of σ₁ 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, σ₁ and σ₂, from the closed-form solution. The numerical results of the two principal stresses, σ₁ and σ₂, 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 σ₂ 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, σ₁ and σ₂, converge to the closed-form solution, as shown in the bottom plot of Fig. 4(b), where max (δs_β) ∼ 0.004R. 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.

Fig. 4

(a) Principal stresses (σ1, σ2) of Brazilian disk compression from the closed-form solution. (b) Principal stresses (σ1, σ2) from the scheme shown in Fig. 3(b), from the top: 101 × 41 data points and max (δsβ) ∼ 0.02R, 201 × 41 data points and max (δsβ) ∼ 0.01R, and 501 × 41 data points and max (δsβ) ∼ 0.004R. (c) Principal stresses (σ1, σ2) from the scheme shown in Fig. 3(c), from the top: 21 × 21 data points and max (δsα) ∼ 0.1R, 41 × 51 data points and max (δsα) ∼ 0.04R, 41 × 101 data points and max (δsα) ∼ 0.02R.

View large Download slide

(a) Principal stresses (σ₁, σ₂) of Brazilian disk compression from the closed-form solution. (b) Principal stresses (σ₁, σ₂) from the scheme shown in Fig. 3(b), from the top: 101 × 41 data points and max (δs_β) ∼ 0.02R, 201 × 41 data points and max (δs_β) ∼ 0.01R, and 501 × 41 data points and max (δs_β) ∼ 0.004R. (c) Principal stresses (σ₁, σ₂) from the scheme shown in Fig. 3(c), from the top: 21 × 21 data points and max (δs_α) ∼ 0.1R, 41 × 51 data points and max (δs_α) ∼ 0.04R, 41 × 101 data points and max (δs_α) ∼ 0.02R.

The numerical results of the two principal stresses, σ₁ and σ₂, 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 σ₁. 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₁), using the Boussinesq solution. By choosing, the L_α(a₁) very close to the loading point will reduce the approximation error. The second factor is the finite difference scheme used for updating σ₂ from L_α(a_i) to L_α(a_i+1), or updating σ₁ 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₂ 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.

Fig. 5

View large Download slide

An infinite sheet with a circular hole and subject to remote tension

The closed-form expression of the stresses in the Cartesian coordinate system (x₁, x₂) is listed below [15]

\begin{aligned} \frac{σ_{11}}{σ} & = \frac{a^{2} (x_{1}^{2} - x_{2}^{2})}{2 (x_{1}^{2} + x_{2}^{2})^{2}} + \frac{a^{2} (x_{1}^{4} - 6 x_{1}^{2} x_{2}^{2} + x_{2}^{4})}{(x_{1}^{2} + x_{2}^{2})^{3}} - \frac{3 a^{4} (x_{1}^{4} - 6 x_{1}^{2} x_{2}^{2} + x_{2}^{4})}{2 (x_{1}^{2} + x_{2}^{2})^{4}} \end{aligned}

(41a)

\begin{aligned} \frac{σ_{22}}{σ} & = 1 + \frac{3 a^{2} (x_{1}^{2} - x_{2}^{2})}{2 (x_{1}^{2} + x_{2}^{2})^{2}} - \frac{a^{2} (x_{1}^{4} - 6 x_{1}^{2} x_{2}^{2} + x_{2}^{4})}{(x_{1}^{2} + x_{2}^{2})^{3}} + \frac{3 a^{4} (x_{1}^{4} - 6 x_{1}^{2} x_{2}^{2} + x_{2}^{4})}{2 (x_{1}^{2} + x_{2}^{2})^{4}} \end{aligned}

(41b)

\begin{aligned} \frac{σ_{12}}{σ} & = - \frac{a^{2} x_{1} x_{2}}{2 (x_{1}^{2} + x_{2}^{2})^{2}} + \frac{4 a^{2} x_{1} x_{2} (x_{1}^{2} - x_{2}^{2})}{(x_{1}^{2} + x_{2}^{2})^{3}} - \frac{6 a^{4} x_{1} x_{2} (x_{1}^{2} - x_{2}^{2})}{(x_{1}^{2} + x_{2}^{2})^{4}} \end{aligned}

(41c)

over the infinite domain

R = {(x_{1}, x_{2}) | x_{1}^{2} + x_{2}^{2} ⩾ a^{2}

⁠. Although we have the closed-form solution for the stress field, unlike the Brazilian disk compression case, other quantities cannot be expressed in simple algebraic forms and we have to rely on numerical calculations.

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 σ₁ and the β-lines (blue or running left and right) to principal stress σ₂. These characteristic curves have some characteristics, especially around the circular cavity, as shown in Fig. 6(b).

Fig. 6

View large Download slide

(a) Contour plot of the principal stress orientation θ_σ and the associated characteristic curves and (b) characteristic curves near the circular cavity edge

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.

Next, we need to determine the traction boundary conditions both far away from the circular cavity and along the edge of the circular hole. The far-field boundary condition can be approximated to be

\begin{aligned} σ_{1} & = σ, | x_{1} | ⩽ H and x_{2} = \pm H \\ σ_{2} & = 0, | x_{2} | ⩽ H and x_{1} = \pm H \end{aligned}

(42)

where the parameter 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

\begin{aligned} σ_{2} & = 0, on \sqrt{x_{1}^{2} + x_{2}^{2}} = a and | x_{2} | ⩽ \sqrt{3} a / 2 \\ σ_{1} & = 0, on \sqrt{x_{1}^{2} + x_{2}^{2}} = a and | x_{1} | ⩽ a / 2 \end{aligned}

(43)

and these boundary conditions are exact.

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.

Fig. 7

View large Download slide

Contour plot of the gradient of orientation field θ_σ in curvilinear coordinate system (s_α, s_β)

The numerical scheme of solving for the principal stresses, σ₁ and σ₂, 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.

Fig. 8

Numerical scheme for solving principal stresses (σ1, σ2)

View large Download slide

Numerical scheme for solving principal stresses (σ₁, σ₂)

We start from the far-field, or far below the circular cavity. Along the first β-line,

L_{β}^{(1)}

⁠, as shown in the upper plot in Fig. 8, we have the approximated boundary condition σ₁ = σ, so that the minor principal stress σ₂ satisfies the PDE,

\begin{aligned} \frac{\partial σ_{2}}{\partial s_{β}} - θ_{σ},_{s_{α}} (s_{α}, s_{β}) σ_{2} = - σ θ_{σ},_{s_{α}} (s_{α}, s_{β}) \\ (s_{α}, s_{β}) \in L_{β}^{(1)}, σ_{2} (s_{α}, s_{β} = 0) = 0 \end{aligned}

(44)

This PDE is integrable and the principal stress σ₂ is solved along the first β-line,

L_{β}^{(1)}

⁠. Now suppose that σ₁ and σ₂ along β-line

L_{β}^{(j)}

are known, we are seeking the solution along β-line

L_{β}^{(j + 1)}

⁠. The major principal stress σ₁ can be updated along the characteristic α-lines as

\begin{aligned} σ_{1}^{(i, j + 1)} & = σ_{1}^{(i, j)} - (σ_{1}^{(i, j)} - σ_{2}^{(i, j)}) θ_{σ},_{s_{β}} (s_{α}^{(i, j)}, s_{β}^{(i, j)}) δ s_{α}^{(i, j)}, \\ i = 1, 2, \dots, M \end{aligned}

(45)

after which the minor principal stress σ₂ satisfies the PDE with some initial conditions,

\begin{aligned} \frac{\partial σ_{2}}{\partial s_{β}} - θ_{σ},_{s_{α}} (s_{α}, s_{β}) σ_{2} & = - θ_{σ},_{s_{α}} (s_{α}, s_{β}) σ^{(j + 1)} (s_{α}, s_{β}), \\ (s_{α}, s_{β}) \in L_{β}^{(j + 1)}, σ_{2} (s_{α}, s_{β} = 0) = 0 \end{aligned}

(46)

and the principal stress σ₂ is solved over the β-line

L_{β}^{(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 σ₁ along α-lines, only solutions (σ₁, σ₂) along a portion of the previous β-line are used, while the solution (σ₁, σ₂) 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 σ₁ 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 σ₁ = 0. Therefore, the principal stress σ₁ 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.

Fig. 9

Contour plots of the closed-form solution for the two principal stresses (σ1, σ2)

View large Download slide

Contour plots of the closed-form solution for the two principal stresses (σ₁, σ₂)

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 = 4a. 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 σ₂ 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.

Fig. 10

The field of principal stresses (σ1, σ2) from the proposed scheme and H = 4a. From the top: with 552 data points and coarse mesh density surrounding the circular hole, with 2568 data points and median mesh density, and with 5022 data points and fine mesh density.

View large Download slide

The field of principal stresses (σ₁, σ₂) from the proposed scheme and H = 4a. From the top: with 552 data points and coarse mesh density surrounding the circular hole, with 2568 data points and median mesh density, and with 5022 data points and fine mesh density.

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 = 4a, to relatively farther away from the hole, H = 8a. 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 σ₂ is positive. As H becomes larger, these two regions become closer to that in the closed-form solution.

Fig. 11

The field of principal stresses (σ1, σ2) from the proposed scheme with median mesh density surrounding the circular hole. From the top: parameter H = 4a, 6a, and 8a, respectively.

View large Download slide

The field of principal stresses (σ₁, σ₂) from the proposed scheme with median mesh density surrounding the circular hole. From the top: parameter H = 4a, 6a, and 8a, respectively.

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.

Fig. 12

View large Download slide

(a) Specimen of PVC sheet with a circular hole and (b) overall response of a PVC sheet specimen

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⁻³s⁻¹. 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₁ represents the displacement in the horizontal direction and u₂ in the vertical direction. The PVC sheet specimen experiences quite large deformation with the axial extension $δ / L \sim 55 %$ 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.

Fig. 13

View large Download slide

Displacement field over the PVC sheet sample surface at selected moments indicated in Fig. 12(b)

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.

Fig. 14

View large Download slide

Deformation gradient field F_αβ (α, β = 1, 2) over the sample surface at moment (D)

From the discussions in Sec. 2 and according to Eqs. (5) and (6), the fields of principal stretches, λ₁ and λ₂, 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.

Fig. 15

Major principal stretch λ1, minor principal stretch λ2, and principal stretch orientation θλ at selected moments indicated in Fig. 12(b)

View large Download slide

Major principal stretch λ₁, minor principal stretch λ₂, and principal stretch orientation θ_λ at selected moments indicated in Fig. 12(b)

The principal stretches, λ₁ and λ₂, and the principal stretch orientation θ_λ will determine the right stretch tensor U according to Eq. (10), and from Eq. (12), the orthogonal tensor Q can be obtained. For 2D deformation of the PVC sheet, Q represents a rigid-body rotation and can be represented by the rotation tensor 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).

Fig. 16

View large Download slide

Principal stress orientation θ_σ field in the current state and at selected moments indicated in Fig. 12(b)

The evolution of the principal stress orientation field θ_σ, as shown in Fig. 16, forms the basis of the technique for determining the stresses associated with the deformations as presented in Fig. 15. Such stress field, which is nonuniform in general, will be determined and studied in Sec. 7.2.

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 σ₂ = 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 ${\bar{λ}}_{1} σ$ ⁠; here, ${\bar{λ}}_{1}$ is the average principal stretch λ₁ along the top and the bottom edges of the gage section, as shown in Fig. 17(a).

Fig. 17

View large Download slide

(a) Traction boundary condition along the external edges of the PVC sheet specimen and (b) traction boundary condition along the edge of the cavity by tracking the variation of the principal stress orientation θ_σ

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 σ₂ = 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 σ₁ = 0. However, the span of the polar angle at the top and the bottom poles of the cavity, where σ₁ = 0, is only about 7.2 deg, so that the majority of the cavity edge is controlled by the boundary condition of σ₂ = 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 σ₁ (α-lines) and the blue lines (running left and right) are the ones associated with the minor principal stress σ₂ (β-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 σ₂ = 0 along the lateral edges of the sheet specimen.

Fig. 18

View large Download slide

(a) Characteristic curves associated with the principal stress orientation field θ_σ shown in Fig. 17(a). (b) Gradient fields of the principal stress orientation θ_σ at moment (D) in the current configuration and projected into the curvilinear coordinate (s_α, s_β).

The last ingredient for solving Eq. is the gradient of the principal stress orientation θ_σ 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:

\nabla_{y} θ_{σ} = (G^{- 1} F) \nabla_{x} θ_{σ}, G = F F^{⊤}

(47)

where

\nabla_{x}

and

\nabla_{y}

denote the gradient operators in the reference and the current states, respectively, and F 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 x₁ and x₂. As a result, the gradient

\nabla_{x} θ_{σ}

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. 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 σ₁ and σ₂ in the PVC sheet. The computed principal stress fields (σ₁, σ₂) 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.

Fig. 19

Computed principal stress fields (σ1, σ2) in PVC sheet specimen at moment (D): (a) 2438 data points, (b) 5994 data points, and (c) 7990 data points

View large Download slide

Computed principal stress fields (σ₁, σ₂) in PVC sheet specimen at moment (D): (a) 2438 data points, (b) 5994 data points, and (c) 7990 data points

Finally, we would like to point out, again, that in determining the principal stress fields, σ₁ and σ₂, 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 ${\bar{λ}}_{1}$ over the top and the bottom edges of the gage section. The principal stress orientation field θ_σ and the average stretch ${\bar{λ}}_{1}$ are the results from the deformation measurement. The principal stress fields, σ₁ and σ₂, 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.

We choose the hyperelastic constitutive model, the generalized neo-Hookean (GNH) model [18], to predict the stresses in the PVC sheet. We further assume that the PVC sheet material is incompressible. Hyper-elasticity assumes that a deforming body possesses an elastic potential

W (x) = W (F (x), x), x \in R_{ref}

(48)

such that the first Piola–Kirchhoff (nominal or engineering) stress field associated with the deformation, τ(x) is given by

τ = \frac{\partial W}{\partial F}, on R_{ref}

(49)

Here, F is the deformation gradient tensor and in the present study, it is experimentally measurable. Further, the Cauchy stress σ is related to the nominal stress τ by

σ = τ F^{⊤}, on R

(50)

For homogeneous material, the elastic potential W will depend on position x only through the deformation gradient tensor F, i.e., W(x) = W(F(x)).

The elastic potential function of the GNH model has the form

W = \frac{μ}{2 b} {{[1 + \frac{b}{n} (I_{1} - 3)]}^{n} - 1}

(51)

where μ > 0, b > 0, and n > 0 are constants, and

I_{1} = λ_{1}^{2} + λ_{2}^{2} + λ_{3}^{2}, λ_{3} = \frac{1}{λ_{1} λ_{2}}

and λ_i (i = 1, 2, 3) are the local principal stretches, where λ₁ ⩾ λ₂ ⩾ λ₃.

For the two-dimensional deformation considered in the present study, the Cauchy stress components σ_αβ, according to the GNH model, are given by

\begin{aligned} σ_{α β} & = μ {1 + \frac{b}{n} (I_{1} - 3)}^{n - 1} (F_{α γ} F_{β γ} - \frac{1}{λ_{1}^{2} λ_{2}^{2}} δ_{α β}), \\ α, β = 1, 2 \end{aligned}

(52)

The two principal stresses σ₁ and σ₂ are given by

\begin{aligned} σ_{1} & = μ {1 + \frac{b}{n} (I_{1} - 3)}^{n - 1} (λ_{1}^{2} - \frac{1}{λ_{1}^{2} λ_{2}^{2}}), \\ σ_{2} & = μ {1 + \frac{b}{n} (I_{1} - 3)}^{n - 1} (λ_{2}^{2} - \frac{1}{λ_{1}^{2} λ_{2}^{2}}) \end{aligned}

(53)

The predictions of Eq. , when the local principal stretches λ₁ and λ₂ 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.

Fig. 20

View large Download slide

(a) Displacement fields in a uniaxial tension specimen made of PVC sheet and (b) variation of nominal axial stress Σ as function of the axial stretch λ from the uniaxial tension tests (solid lines) and fitting result of the GNH model (open diamonds).

For uniaxial tension, the GNH model relates the axial nominal stress Σ to the axial stretch λ through

Σ = μ {1 + \frac{b}{n} (λ^{2} + \frac{2}{λ} - 3)}^{n - 1} (λ - \frac{1}{λ^{2}})

(54)

Therefore, one can fit the expression in Eq. to the experimental data shown in Fig. 20(b) to determine the material constants μ, 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 λ₁ and λ₂ have been experimentally measured, as shown in Fig. 15. The predicted stress fields σ₁ and σ₂, based on the GNH model, are presented in Fig. 21 for the selected moments indicated in Fig. 12(b).

Fig. 21

View large Download slide

Stress field in the PVC sheet based on the GNH model prediction and at selected moments indicated in Fig. 12(b)

Figure 21 shows only the two principal stresses σ₁ and σ₂; however in general, three components are required to describe the stress field in a planar deformation. By examining Eq. , we note that

\begin{aligned} σ & = μ {1 + \frac{b}{n} (I_{1} - 3)}^{n - 1} (F F^{⊤} - \frac{1}{λ_{1}^{2} λ_{2}^{2}} I) \\ = μ {1 + \frac{b}{n} (I_{1} - 3)}^{n - 1} Q (U^{2} - \frac{1}{λ_{1}^{2} λ_{2}^{2}} I) Q^{⊤} \\ = μ {1 + \frac{b}{n} (I_{1} - 3)}^{n - 1} R (θ_{λ} + θ_{q}) (Λ^{2} - \frac{1}{λ_{1}^{2} λ_{2}^{2}} I) R^{⊤} (θ_{λ} + θ_{q}) \end{aligned}

where we have used Eqs. and , the notion Q = R(θ_q) for planar deformation, and the property of the rotation tensor R. Combining the above expression with Eq. , one will reach Eq. . Therefore, the principal stress orientation θ_σ = θ_λ + θ_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, σ₁ and σ₂, 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).

Fig. 22

Computed principal stress field (σ1, σ2) in the PVC sheet at selected moments indicated in Fig. 12(b)

View large Download slide

Computed principal stress field (σ₁, σ₂) in the PVC sheet at 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 $R_{ref}$ ⁠. The boundaries or edges of domain $R_{ref}$ 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.

Fig. 23

View large Download slide

Various domains in the determination of the stress field associated with a nonuniform deformation field: (a) Reference domain $R_{ref}$ over which the principal stretch orientation θ_λ is measured. (b) Current domain $R$ where the orientation of principal stress orientation θ_σ is measured, and the computational domain where the stresses σ₁ and σ₂ are determined.

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 σ₁, 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 σ₂ 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, σ₂ is small but not zero as the boundary condition requires. Along the boundary of the cavity, σ₂ 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₁ along the top and the bottom edge of the gage section is approximately zero and the vertical displacement u₂ 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.

Again, take the PVC sheet with a circular hole and subject to uniaxial tension as the example. The traction boundary condition along the lateral edges of the PVC sheet is exact, while the traction boundary condition over the top and the bottom edges of the computational domain is approximate. The top and the bottom edges are two characteristic β-lines,

L_{β}^{(1)}

and

L_{β}^{(N)}

⁠, over which major principal stress σ₁ needs to be described as additional boundary conditions. We formally write the boundary condition as

\begin{aligned} σ_{1} (s_{α}, s_{β}) & = f_{1} (s_{α}, s_{β}), (s_{α}, s_{β}) \in L_{β}^{(1)} \\ σ_{1} (s_{α}, s_{β}) & = f_{N} (s_{α}, s_{β}), (s_{α}, s_{β}) \in L_{β}^{(N)} \end{aligned}

where f₁(s_α, s_β) and f_N(s_α, s_β) are yet to be determined, but they should satisfy the constraint

\begin{aligned} \int_{0}^{S_{β}^{(1)}} f_{1} (s_{α}, s_{β}) \sin θ_{σ} (s_{α}, s_{β}) d s_{β} & = {\bar{λ}}_{1} {\bar{λ}}_{2} σ W \\ \int_{0}^{S_{β}^{(N)}} f_{N} (s_{α}, s_{β}) \sin θ_{σ} (s_{α}, s_{β}) d s_{β} & = {\bar{λ}}_{1} {\bar{λ}}_{2} σ W \end{aligned}

Here, W is the original width of the PVC sheet sample, and

S_{β}^{(1)}

and

S_{β}^{(N)}

are the total lengths of the β-lines,

L_{β}^{(1)}

and

L_{β}^{(N)}

⁠, respectively.

Furthermore, we assume that f₁(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

\begin{aligned} f_{1} (s_{α}, s_{β}) & = \sum_{i = 1}^{K} a_{i}^{(1)} ξ_{i} (s_{α}, s_{β}) \\ f_{N} (s_{α}, s_{β}) & = \sum_{i = 1}^{K} a_{i}^{(N)} ξ_{i} (s_{α}, s_{β}) \end{aligned}

For any given set of coefficients

a_{i}^{(1)}

and

a_{i}^{(N)}

(i = 1, 2, …, K), following the scheme discussed in the present investigation, we may determine a stress field

\begin{aligned} σ_{1} (y_{1}, y_{2}; a_{i}^{(1)}, a_{i}^{(N)}, i = 1, 2, \dots, K) \\ σ_{2} (y_{1}, y_{2}; a_{i}^{(1)}, a_{i}^{(N)}, i = 1, 2, \dots, K) \end{aligned}

Together with the experimentally measured deformation field, e.g., λ₁ and λ₂, we may define an energetic measure

Π (a_{i}^{(1)}, a_{i}^{(N)}, i = 1, 2, \dots, K)

⁠. Now we may invoke some variational principles to choose the set of coefficients

a_{i}^{(1)}

and

a_{i}^{(N)}

(i = 1, 2, …, K) such that the energetic measure Π is minimized. Once the set of coefficients

a_{i}^{(1)}

and

a_{i}^{(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_{β}^{(1)}$ and $L_{β}^{(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 σ₁ and σ₂, 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.

The second area of extension is the case of dynamic or impact loading. For dynamic deformations, where the inertia of material particles cannot be neglected, we need to consider the equation of motion. The equation of motion, in the absence of body force density and in terms of the two principal stresses σ₁ and σ₂, and the principal stress orientation θ_σ are given by

\begin{aligned} \frac{\partial σ_{1}}{\partial y_{1}} + \tan θ_{σ} \frac{\partial σ_{1}}{\partial y_{2}} & = (\tan θ_{σ} \frac{\partial θ_{σ}}{\partial y_{1}} - \frac{\partial θ_{σ}}{\partial y_{2}}) (σ_{1} - σ_{2}) \\ + ρ ({\ddot{u}}_{1} + \tan θ_{σ} {\ddot{u}}_{2}) \end{aligned}

(55a)

\begin{aligned} \frac{\partial σ_{2}}{\partial y_{1}} - \cot θ_{σ} \frac{\partial σ_{2}}{\partial y_{2}} & = (\cot θ_{σ} \frac{\partial θ_{σ}}{\partial y_{1}} + \frac{\partial θ_{σ}}{\partial y_{2}}) (σ_{1} - σ_{2}) \\ + ρ ({\ddot{u}}_{1} - \cot θ_{σ} {\ddot{u}}_{2}) \end{aligned}

(55b)

where

(y_{1}, y_{2}) \in R

and t ⩾ 0, ρ is the mass density, and

{\ddot{u}}_{1}

and

{\ddot{u}}_{2}

are the particle acceleration components in the (y₁, y₂) coordinate system. Again, from the left-hand side of Eq. , we see that there exist two family of characteristic curves governed by the two ODEs given by Eq. . By introducing the curvilinear coordinates (s_α, s_β) like in the quasi-static case, the equation of motion can be recast into

\begin{aligned} \frac{\partial σ_{1}}{\partial s_{α}} + (σ_{1} - σ_{2}) θ_{σ},_{s_{β}} (s_{α}, s_{β}) & = ρ {\ddot{u}}_{s_{α}} (s_{α}, s_{β}) \\ \frac{\partial σ_{2}}{\partial s_{β}} + (σ_{1} - σ_{2}) θ_{σ},_{s_{α}} (s_{α}, s_{β}) & = ρ {\ddot{u}}_{s_{β}} (s_{α}, s_{β}) \end{aligned}

(56)

where

{\ddot{u}}_{s_{α}} = {\ddot{u}}_{1} \cos θ_{σ} + {\ddot{u}}_{2} \sin θ_{σ}

and

{\ddot{u}}_{s_{β}} = - {\ddot{u}}_{1} \sin θ_{σ} + {\ddot{u}}_{2} \cos θ_{σ}

⁠, i.e.,

({\ddot{u}}_{s_{α}}, {\ddot{u}}_{s_{β}})

is the projection of the acceleration vector

({\ddot{u}}_{1}, {\ddot{u}}_{2})

onto the curvilinear coordinates (s_α, s_β). Equation is the same as Eq. except the additional terms on the right-hand side. From the experimental perspective, at any given moment during the deformation, the displacement fields u₁ and u₂ 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

{\dot{u}}_{1}

and

{\dot{u}}_{2}

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

{\ddot{u}}_{1}

and

{\ddot{u}}_{2}

can be computed. Therefore, the numerical scheme discussed in this investigation can be used to solve the equation of motion , 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.

Assuming that the 3D displacement field u(x) is measured for

x \in R

⁠, as a result, deformation gradient tensor F can be computed. The three principal stretches λ_i (i = 1, 2, 3), where λ₁ ⩾ λ₂ ⩾ λ₃, 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 U can be written as

U = R (ϕ_{λ}, θ_{λ}, ψ_{λ}) Λ R^{⊤} (ϕ_{λ}, θ_{λ}, ψ_{λ}) Λ = (\begin{matrix} λ_{1} & 0 & 0 \\ 0 & λ_{2} & 0 \\ 0 & 0 & λ_{3} \end{matrix})

where R is the rotation tensor and (φ_λ, θ_λ, ψ_λ) are the Euler angles associated with the right stretch tensor U. The proper orthogonal tensor Q = FU⁻¹ representing a local rotation and, therefore, can be written as

Q = R (ϕ_{q}, θ_{q}, ψ_{q})

and (φ_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

σ = R (ϕ_{σ}, θ_{σ}, ψ_{σ}) Σ R^{⊤} (ϕ_{σ}, θ_{σ}, ψ_{σ}), Σ = (\begin{matrix} σ_{1} & 0 & 0 \\ 0 & σ_{2} & 0 \\ 0 & 0 & σ_{3} \end{matrix})

where σ_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

R (ϕ_{σ}, θ_{σ}, ψ_{σ}) = R (ϕ_{q}, θ_{q}, ψ_{q}) R (ϕ_{λ}, θ_{λ}, ψ_{λ})

As a result, the principal stress directions of the Cauchy stress σ 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.

References

1.

Strang

,

G.

, and

Fix

,

G.

,

2008

,

An Analysis of The Finite Element Method

, 2nd ed.,

Wellesley-Cambridge Press

,

BerlinWellesley, MA

.

2.

Reddy

,

J.

,

2006

,

An Introduction to the Finite Element Method

, 3rd ed.,

McGraw-Hill

,

New York

.

3.

Zienkiewicz

,

O.

,

Taylor

,

R.

, and

Zhu

,

J.

,

2013

,

The Finite Element Method: Its Basis and Fundamentals

, 7th ed.,

Butterworth-Heinemann

,

Woburn, MA

.

4.

Grédiac

,

M.

,

Toussaint

,

E.

, and

Pierron

,

F.

,

2002

, “

Applying the Virtual Fields Method to the Identification of Elasto-Plastic Constitutive Parameters

,”

Int. J. Solids. Struct.

,

39

(

10

), pp.

2691

–

2705

. 10.1016/S0020-7683(02)00127-0

Google Scholar

Crossref

5.

Grédiac

,

M.

,

Toussaint

,

E.

, and

Pierron

,

F.

,

2002

, “

Special Virtual Fields for the Direct Determination of Material Parameters With the Virtual Fields Method. 2–Application to In-Plane Properties

,”

Int. J. Solids. Struct.

,

39

(

10

), pp.

2707

–

2730

. 10.1016/S0020-7683(02)00128-2

Google Scholar

Crossref

6.

Grédiac

,

M.

, and

Pierron

,

F.

,

2006

, “

Applying the Virtual Fields Method to the Identification of Elasto-Plastic Constitutive Parameters

,”

Int. J. Plast.

,

22

(

4

), pp.

602

–

627

. 10.1016/j.ijplas.2005.04.007

Google Scholar

Crossref

7.

Cameron

,

B.

, and

Tasan

,

C.

,

2021

, “

Full-field Stress Computation From Measured Deformation Fields: A Hyperbolic Formulation

,”

J. Mech. Phys. Solids.

,

147

, p.

104186

. 10.1016/j.jmps.2020.104186

Google Scholar

Crossref

8.

Chadwick

,

P.

,

1999

,

Continuum Mechanics: Concise Theory and Problems

, 2nd ed.,

Dover

,

New York, NY

.

9.

Itskov

,

M.

,

2019

,

Tensor Algebra and Tensor Analysis for Engineers with Applications to Continuum Mechanics

, 5th ed.,

Springer Nature Switzerland AG

,

Cham Switzerland

.

Google Scholar

Crossref

10.

Ogden

,

R.

,

1997

,

Non-Linear Elastic Deformations

,

Dover

,

New York, NY

.

11.

Holzapfel

,

G.

,

2000

,

Nonlinear Solid Mechanics: A Continuum Approach for Engineering

,

John Wiley & Sons Ltd.

,

West Sussex, England

.

12.

Hondros

,

G.

,

1959

, “

The Evaluation of Poisson’s Ratio and the Modulus of Materials of a Low Tensile Resistance by the Brazilian (Indirect Tensile) Test With Particular Reference to Concrete

,”

Australian J. Appl. Sci.

,

50

, pp.

243

–

268

.

13.

Liu

,

C.

, and

Thompson

,

D.

,

2012

, “

Deformation and Failure of a Heterogeneous High Explosive

,”

Philos. Mag. Lett.

,

92

(

8

), pp.

352

–

361

. 10.1080/09500839.2012.673021

Google Scholar

Crossref

14.

Liu

,

C.

,

2010

, “

Elastic Constants Determination and Deformation Observation Using Brazilian Disk Geometry

,”

Exp. Mech.

,

50

(

7

), pp.

1025

–

1039

. 10.1007/s11340-009-9281-2

Google Scholar

Crossref

15.

Timoshenko

,

S.

, and

Goodier

,

J.

,

1951

,

Theory of Elasticity

, 2nd ed.,

McGraw Hill

,

New York, NY

.

16.

Sutton

,

M.

,

McNeill

,

S.

,

Helm

,

J.

, and

Chao

,

Y.

,

2000

, “

Advances in Two-Dimensional and Three-Dimensional Computer Vision

,”

Top. Appl. Phys.

,

77

, pp.

323

–

372

.

Google Scholar

Crossref

17.

Sutton

,

M.

,

Orteu

,

J.

, and

Schreier

,

H.

,

2009

,

Image Correlation for Shapes, Motion, and Deformation Measurement: Basic Concepts, Theory and Applications

,

Springer

,

New York, NY

.

18.

Knowles

,

J.

,

1977

, “

The Finite Anti-Plane Shear Field Near the Tip of Crack for a Class of Incompressible Elastic Solids

,”

Int. J. Fract.

,

13

(

5

), pp.

611

–

639

. 10.1007/BF00017296

Google Scholar

Crossref

19.

Franck

,

C.

,

Hong

,

S.

,

Maskarinec

,

S.

,

Tirrell

,

D.

, and

Ravichandran

,

G.

,

2007

, “

Three-Dimensional Full-Field Measurements of Large Deformation in Soft Materials Using Confocal Microscopy and Digital Volume Correlation

,”

Exp. Mech.

,

47

(

3

), pp.

427

–

438

. 10.1007/s11340-007-9037-9

Google Scholar

Crossref

20.

Hu

,

Z.

,

Luo

,

H.

,

Bardenhagen

,

S.

,

Siviour

,

C.

,

Armstrong

,

R.

, and

Lu

,

H.

,

2015

, “

Internal Deformation Measurement of Polymer Bonded Sugar in Compression by Digital Volume Correlation of In-Situ Tomography

,”

Exp. Mech.

,

55

(

1

), pp.

289

–

300

. 10.1007/s11340-014-9856-4

Google Scholar

Crossref

21.

Croom

,

B.

,

Burden

,

D.

,

Jin

,

H.

,

Vonk

,

N.

,

Hoefnagels

,

J.

,

Smaniotto

,

B.

,

Hild

,

F.

,

Quintana

,

E.

,

Sun

,

Q.

,

Nie

,

X.

, and

Li

,

X.

,

2021

, “

Internal Deformation Measurement of Polymer Bonded Sugar in Compression by Digital Volume Correlation of In-Situ Tomography

,”

Exp. Mech.

,

61

(

2

), pp.

395

–

410

. 10.1007/s11340-020-00653-x

Google Scholar

Crossref

Nonuniform Stress Field Determination Based on Deformation Measurement

Abstract

1 Introduction

2 Finite Deformation

3 Principal Stretch Orientation and Principal Stress Orientation

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

3.2 Notion of Tensor-Valued Isotropic Functions.

3.3 Work Conjugate and the Second Piola–Kirchhoff Stress.

3.4 Cauchy Stress and Principal Stress Orientation.

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

4.1 Equation of Equilibrium.

4.2 Characteristic Curves and Recast Partial Differential Equations.

5 Numerical Scheme

6 Verification and Validation

6.1 Virtual Test Case I: Brazilian Disk Compression.

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

7 Uniaxial Tension of Polyvinyl Chloride Sheet With a Circular Hole

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

7.2 Stress Field in Polyvinyl Chloride Sheet.

7.3 Comparison With Constitutive Model Predictions.

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

8 Concluding Remarks

Acknowledgment

Conflict of Interest

Data Availability Statement

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Nonuniform Stress Field Determination Based on Deformation Measurement

Abstract

1 Introduction

2 Finite Deformation

3 Principal Stretch Orientation and Principal Stress Orientation

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

3.2 Notion of Tensor-Valued Isotropic Functions.

3.3 Work Conjugate and the Second Piola–Kirchhoff Stress.

3.4 Cauchy Stress and Principal Stress Orientation.

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

4.1 Equation of Equilibrium.

4.2 Characteristic Curves and Recast Partial Differential Equations.

5 Numerical Scheme

6 Verification and Validation

6.1 Virtual Test Case I: Brazilian Disk Compression.

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

7 Uniaxial Tension of Polyvinyl Chloride Sheet With a Circular Hole

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

7.2 Stress Field in Polyvinyl Chloride Sheet.

7.3 Comparison With Constitutive Model Predictions.

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

8 Concluding Remarks

Acknowledgment

Conflict of Interest

Data Availability Statement

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only