In this investigation, we demonstrate a technique that, under certain circumstances, will determine stresses associated with a nonuniform deformation field without knowing the detailed constitutive behavior of the deforming material. This technique is based on (1) a detailed deformation measurement of a domain (currently in 2D) and (2) the observation that for isotropic materials, the strain and the stress, which form the so-called work-conjugate pair, are co-axial, or their eigenvectors share the same directions. The particular measures for strain and stress chosen in this study are the Lagrangian (or Green-Lagrangian) strain and the second Piola–Kirchhoff stress. The deformation measurement provides the field of the principal stretch orientation θλ and since the Lagrangian strain and the second Piola–Kirchhoff stress are co-axial, the principal stress orientation θs of the second Piola–Kirchhoff stress is then determined. The Cauchy stress is related to the second Piola–Kirchhoff stress through the deformation gradient tensor, which can be measured experimentally. We then show that the principal stress orientation θσ of the Cauchy stress is the sum of the principal stretch orientation θλ and the local rigid-body rotation θq, which is determinable by the deformation gradient through polar decomposition. Such a relationship is valid for finite deformations. With the principal stress orientation θσ known, the equation of equilibrium, now in terms of the two principal stresses, σ1 and σ2, and θσ, can be solved numerically with appropriate traction boundary conditions. The stresses determined using this technique obviously satisfy the equation of equilibrium, in contrast to those obtained from a constitutive model with input from deformation measurement. The technique and the associated numerical scheme are verified and validated through two virtual test cases representative of the simply-connected and multiply-connected domains, where exact solutions are available. The technique is then applied to an experimental case of nonuniform deformation of a polyvinyl chloride (PVC) sheet with a circular hole and subject to uniaxial tension. In this case, the associated stress field is also determined through a constitutive model of hyperelasticity, the generalized neo-Hookean (GNH) model, calibrated for the particular PVC sheet. Limitations and restrictions of the technique and the associated numerical scheme, as well as possible extensions will be discussed.
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.
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 , 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. , 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
From experimental perspective, all the quantities discussed so far can be measured in an experiment over the sample domain . For example, using the optical technique of digital image correlation (DIC), one can measure the displacement field uα over the region . 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 (x1, x2) and for any arbitrary angle θ from Eq. (5), and this will lead to the determination of the principal stretches λ1 and λ2, and the principal stretch orientation θλ at any given location (x1, x2).
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.
To summarize, from experiment, we can measure the displacement field uα(x1, x2) (α = 1, 2), the deformation gradient tensor F, the principal stretches, λ1 and λ2, and the principal stretch orientation θλ, and consequently, we can also determine the components involved in the polar decomposition of the deformation gradient tensor F, U, and Q.
3.2 Notion of Tensor-Valued Isotropic Functions.
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. . 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.
3.4 Cauchy Stress and Principal Stress Orientation.
4 Equation of Equilibrium, Characteristic Lines, and Recast Partial Differential Equation
4.1 Equation of Equilibrium.
Such experimentally measured field θσ(y1, y2) imposes a constraint to the PDEs in Eq. (26), where we have two unknown functions, the principal stresses σ1 and σ2, and two PDEs. With appropriate traction boundary conditions, the stress field, (σ1, σ2), which in general is nonuniform, can be determined through numerical solutions to the PDEs in Eq. (26) over the given domain .
4.2 Characteristic Curves and Recast Partial Differential Equations.
Note that at each location in domain , the α-line and the β-line passing through that location are perpendicular to each other since the product of their slopes, dy2/dy1, 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 (y1, y2). Unlike the Cartesian coordinates (y1, y2), where there is a well defined origin y1 = y2 = 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 (y1, y2), the curvilinear coordinate (sα, sβ) should form a right-handed system.
Equation (28) represents the simplest form of the equation of equilibrium in a 2D deformation. Our subsequent exploration of solving the equation of equilibrium, together with some traction boundary conditions, for the principal stresses, σ1 and σ2 once the functions and are given, will be based on Eq. (28).
5 Numerical Scheme
Although Eqs. (30) and (32) indicate both the existence and the form of solutions for the two principal stresses, σ1 and σ2, under the given assumptions, for practical applications, we will not directly use Eq. (30) or (32). We will use a numerical form of Eqs. (29) and (31) to obtain the principal stresses, σ1 and σ2, at discrete locations inside the deformed domain .
Such a process is thus repeated for all characteristic lines , i = 1, 2, …, M. The principal stresses, σ1 and σ2, are therefore obtained at all locations , 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 , where the value of the principal stress σ1 is known due to the fact that is part of the boundary of the domain and the traction boundary condition for σ1 is given. Then, the principal stress σ2 can be solved by the equation over the line . Using the similar updating scheme and repeating for all the β-lines (j = 1, 2, …, N), the principal stresses, σ1 and σ2, are therefore obtained at all locations , 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 Brazilian disk compression test has a characteristic feature that all three stress components vanish along the disk boundary except at the two ends of the vertical diameter. Along the traction-free boundary, since the shear stress is already zero, one of the principal stresses must be zero as well. Such principal stress can be identified by tracking the variation of the principal stress orientation θσ along the boundary. There are only two choices for the principal stress orientation θσ along any traction-free boundary, either it is pointing to the normal direction of the boundary or it is pointing to the tangential direction of the boundary. If θσ is in the normal direction of the boundary, then σ1 = 0. As shown in both plots in Fig. 1(b), we should have the boundary condition σ1 = 0 everywhere along the circular edge of the disk, except at the two poles and at the two poles, the minor principal stress σ2 = σ = −F2/R.
We can also switch the role of α- and β-lines and start from the first β-line Lβ(b1), where σ1 ≡ 0 according to the traction boundary condition. The minor principal stress σ2 over Lβ(b1) can be obtained by solving the integrable PDE with the end condition given by the Boussinesq solution. The update of the major principal stress σ1 along α-lines Lα(a) will lead to the value of σ1 along the next β-line Lβ(b) similar to that shown in Fig. 3(b). This procedure is then repeated for all β-lines. Such a scheme is illustrated in Fig. 3(c).
Figure 4(a) presents the principal stress fields, σ1 and σ2, from the closed-form solution. The numerical results of the two principal stresses, σ1 and σ2, according to the proposed scheme, indicated in Fig. 3(b), are presented in Fig. 4(b). The top plot in Fig. 4(b) shows the result where 101 α-lines and 41 β-lines are used, and the maximum length of the line segment δsβ, used in the σ2 updating scheme, is about 2% of the disk radius R. The comparisons presented in Fig. 4(b) highlight the important role of the maximum length of the line segment δsβ, or step size, in the proposed numerical scheme. Only when the step size δsβ is small enough that the numerical solutions of the two principal stresses, σ1 and σ2, converge to the closed-form solution, as shown in the bottom plot of Fig. 4(b), where max (δsβ) ∼ 0.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.
The numerical results of the two principal stresses, σ1 and σ2, according to the alternative scheme shown in Fig. 3(c), are presented in Fig. 4(c). In the top plot of Fig. 4(c), we only use very few data points, but we see that the numerical result is already quite good, though slightly off symmetry in σ1. In the middle and the bottom plots of Fig. 4(c), we moderately increase the number of data points, and the numerical results become much better. Also, we see that the smallest step size is only about 2% of the disk radius R, but in the previous updating scheme, we have to use much smaller step size to reach the same accuracy.
The virtual testing case of Brazilian disk compression is unique in that every quantity for checking the proposed scheme can be written in closed form and these include each individual characteristic curve. As a result, it can isolate the factors that might affect the accuracy of the numerical scheme easily. We have shown that using the curvilinear coordinate formed by the two-family of characteristic curves makes the PDE and boundary condition much simpler and cleaner, since the boundary is one of the characteristic curves. The determination of the characteristic curves only requires the experimentally obtainable principal stress orientation θσ and based on the curvilinear coordinate system, we can easily convert solving the coupled PDEs to solving a series of integrable PDEs along one family of the characteristic curves.
There are two factors that control the accuracy and efficiency of the numerical scheme. The first factor is the approximation the point-force boundary condition by a traction along the α-line Lα(a1), using the Boussinesq solution. By choosing, the Lα(a1) very close to the loading point will reduce the approximation error. The second factor is the finite difference scheme used for updating σ2 from Lα(ai) to Lα(ai+1), or updating σ1 from Lβ(bj) to Lβ(bj+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 x2 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.
The field of the principal stress orientation θσ is shown in Fig. 6(a) for the virtual testing case. For the majority of the domain , θσ is around 90 deg, only near the edge of the circular cavity, θσ varies away from 90 deg. Since we already knew/measured the principal stress orientation θσ, the characteristic curves can be obtained by solving the two ODEs in Eq. (27), and they are presented also in Fig. 6(a). Again, we assign the α-lines (red or running up and down) to be associated with the principal stress σ1 and the β-lines (blue or running left and right) to principal stress σ2. These characteristic curves have some characteristics, especially around the circular cavity, as shown in Fig. 6(b).
The β-line far below the hole is approximately running from left to right and the farther away from the hole, the flatter the β-line is. As the β-line approaches the circular hole, which is traction-free along its boundary, the middle section of the β-line is “attracted” toward the circular hole and bent. At some point, the β-line touches the edge of the circular hole and wraps around a section of the cavity, so that the lower section of the circular cavity becomes part of the characteristic β-line. This section ends at two points, where the polar angles are θ = −60 deg and −120 deg, respectively, at which the β-line has sharp kinks and then turns away from the hole. Further up, the β-line becomes two branches on either side of the cavity, and they can be thought of being semi-infinite, starting from the edge of the circular cavity and extending to infinity. Also at points that the β-lines touch the hole edge, these β-lines are all perpendicular to the hole boundary. At locations where the polar angles θ = 60 deg and 120 deg, the two-branched β-lines merges into one and again, wraps around portion of the cavity boundary. Afterward, the β-line departs the hole and moves upward to infinity.
Similar observations can be made of the characteristic α-lines, from the right side of the hole to the left side of the hole. The α-line also touches the hole edge and wraps around the hole boundary between the polar angles θ = ±60 deg, and then becomes two branches, and finally merges into one again. Note that the characteristic curves around the circular cavity boundary are closely related to how the traction boundary conditions are described, i.e., the traction boundary condition is prescribed for the quantity whose characteristic curves are perpendicular to the boundary, not the one whose characteristic curves are part of the boundary.
The contour plots of the gradient of the principal stress orientation θσ in the curvilinear coordinate (sα, sβ), formed by the two family of characteristic curves, are shown in Fig. 7.
The numerical scheme of solving for the principal stresses, σ1 and σ2, is similar to that used in the Brazilian disk compression test. But, there are some differences due to the fact that the deforming domain is multiply connected. The similarities and differences of the numerical scheme are illustrated in Fig. 8 and discussed here.
However, the virtual test of an infinite sheet with a circular hole and subject to remote tension deals with multiply-connected domain. At some point, portion of the hole boundary becomes part of the β-line, where the next β-line has two branches, as illustrated in the bottom plot in Fig. 8. As a result, when we update the principal stress σ1 along α-lines, only solutions (σ1, σ2) along a portion of the previous β-line are used, while the solution (σ1, σ2) over part of the circular boundary is not used. The solution update scheme is carried out on the left and the right sides of the circular cavity, as seen in the bottom plot in Fig. 8, until the moment when the next β-line becomes a single line extending to infinity on both sides. The principal stress σ1 on the left and the right sides of the β-line is updated according to our scheme. However, according to the traction boundary condition, along the boundary of the circular hole, which is also part of the characteristic β-line, we have σ1 = 0. Therefore, the principal stress σ1 along the entire β-line is known. Now, the scheme illustrated in the middle plot in Fig. 8 can be applied for the rest of the domain .
We now compare the numerical results with the closed-form solution for the virtual test problem, where the closed-form solution for the two principal stresses is shown in Fig. 9.
We first study the effect of data density or mesh size in the region near the circular hole. We fix the far-field boundary at the distance of H = 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 σ2 from the closed-form solution, shown in Fig. 9, near the left and the right sides of the circular hole with those from the numerical solutions.
The more important factor that affects the accuracy of the numerical solution is the parameter H, i.e., at which location we approximate the far-field boundary condition. In Fig. 11, we fix the data density in the region around the circular hole and change the parameter H from relatively close to the hole, H = 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 σ2 is positive. As H becomes larger, these two regions become closer to that in the closed-form solution.
In many experimental measurements, the experimental sample design involves multiply-connected domain, for example, specimen with an internal crack and specimen with damage characterized by cavity of certain shape. The virtual test of an infinite sheet with a circular hole and subject to remote tension represents such class of specimen designs and test configurations. By tracking the variation of the experimentally measurable principal stress orientation θσ, the traction boundary condition along the edge of the cavity can be identified and such boundary conditions are exact. The far-field boundary condition can be approximated through remotely applied total force. Through careful sample design and the placement of far-field boundary, error due to approximated far-field boundary condition can be reduced and the accuracy of the numerical results improved, as shown in Fig. 11. However, the characteristic curves associated with the multiply-connected domain become more complicated surrounding the cavity and denser data grid is required.
We have transformed the solving of a set of PDEs into solving a series of integrable PDEs along one family of the characteristic curves and then updating the solution along the other family of the characteristic curves. When solving an integrable PDE along a given characteristic curve, it only involves “initial” condition, or boundary condition from one side of the curve. As a result, advantage can be gained by starting the numerical scheme from both sides of the domain, top and bottom or left and right, and the final results then meet somewhere in the middle. Such an approach will minimize the error when the finite difference scheme, for updating the solution, moves around the cavity.
7 Uniaxial Tension of Polyvinyl Chloride Sheet With a Circular Hole
In this section, an experimental investigation of the nonuniform deformation and the associated nonuniform stress field is carried out on a PVC sheet with a circular hole and subject to uniaxial tension. We use the optical technique of DIC [16,17] to measure the nonuniform displacement field over the sample surface and then apply the technique and the numerical scheme, discussed in Secs. 3–5, to determine the nonuniform stress field in the sample. This stress field is then compared with the prediction using a hyperelasticity constitutive model, the generalized neo-Hookean (GNH) model, which has been calibrated through uniaxial tension test for the PVC sheet.
7.1 Uniaxial Tension Test, Overall Response, and Deformation Measurement.
The material used in this experiment, PVC is a thermoplastic polymer. The sketch of the PVC sheet specimen is shown in Fig. 12(a) and the nominal dimensions of the sample are the gage length L = 20.84 mm, the width W = 37.50 mm, the diameter of the circular hole d = 5.31 mm, and the thickness of the sheet B = 0.254 mm.
We used an Instron test machine (model 1125) to apply the tensile load to the specimen and the load cell capacity is 500 N. The loading speed was controlled by the crosshead velocity of the test machine and it was set at 1.2 mm/min, which gives the equivalent strain rate of about 1.0 × 10−3s−1. The overall response of one of the PVC sheet specimens is presented in Fig. 12(b), where the applied tensile load P is normalized by the cross-sectional area of the sheet A = W B and the gage section elongation δ is normalized by the original gage length L.
Optical DIC is used for mapping the displacement field over the PVC sheet sample surface. Such displacement measurement is shown in Fig. 13 for the selected moments indicated in Fig. 12(b). In the contour plots of Fig. 13, u1 represents the displacement in the horizontal direction and u2 in the vertical direction. The PVC sheet specimen experiences quite large deformation with the axial extension at moment (D), while at the same time, the cavity at the center of the sheet specimen, which starts as a circle, is elongated over 220% in the vertical direction and contracts about 8% in the lateral direction. Meanwhile, both the top and the bottom edges of the gage section remain approximately straight during the tensile process.
With the displacement field measured, deformation gradient can be computed through numerical differentiation. The deformation gradient field Fαβ (α, β = 1, 2) at moment (D) is shown in Fig. 14.
From the discussions in Sec. 2 and according to Eqs. (5) and (6), the fields of principal stretches, λ1 and λ2, as well as the principal stretch orientation θλ can be determined. These fields are presented in Fig. 15 for the selected moments indicated in Fig. 12(b). For the isotropic PVC sheet, θλ also represents the principal stress orientation of the second Piola–Kirchhoff stress S as discussed in Sec. 3.3.
The principal stretches, λ1 and λ2, and the principal stretch orientation θλ will determine the right stretch tensor U according to Eq. (10), and from Eq. (12), the orthogonal tensor 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).
7.2 Stress Field in Polyvinyl Chloride Sheet.
To determine the stress field by solving the equation of equilibrium, Eq. (28), the traction boundary conditions along the external and the cavity edges of the PVC sheet specimen are required. Along the external side edges of the sheet, we have traction-free boundary condition. By noting the principal stress orientation θσ along the side edges, we have the condition σ2 = 0. Note that the nominal applied load at the top and the bottom edges of the gage section is σ, but in the current configuration, this applied load is represented by ; here, is the average principal stretch λ1 along the top and the bottom edges of the gage section, as shown in Fig. 17(a).
Along the edge of the cavity, we also have traction-free boundary condition, but we need to identify which principal stress vanishes over which portion of the boundary. Figure 17(b) presents the variation of the principal stress orientation θσ as a function of the polar angle θ over the edge of the cavity. The polar angle θ is measured from the center of the cavity. In Fig. 17(b), the tangential direction of the cavity edge is shown as dashed line, and we see that for the most part, the principal stress orientation θσ tracks the tangential direction, which indicates that the minor principal stress σ2 = 0 over the portion of the cavity edge that θσ aligns with the tangential direction. Only between the two open circles for either θ < 0 or θ > 0, i.e., near the top and the bottom poles of the cavity, that the principal stress orientation θσ deviates from the tangential direction and aligns with the normal direction of the edge. Therefore, near the top and the bottom poles of the cavity, the boundary condition becomes σ1 = 0. However, the span of the polar angle at the top and the bottom poles of the cavity, where σ1 = 0, is only about 7.2 deg, so that the majority of the cavity edge is controlled by the boundary condition of σ2 = 0.
The other ingredient for solving Eq. (28) is the curvilinear coordinate system (sα, sβ) formed by the two-family of characteristic curves obtained by solving the two ODEs in Eq. (27). For the principal stress orientation θσ field shown in Fig. 17(a) and at moment (D), the characteristic curves are given in Fig. 18(a). In Fig. 18(a), the red lines (running up and down) are the characteristic curves associated with the major principal stress σ1 (α-lines) and the blue lines (running left and right) are the ones associated with the minor principal stress σ2 (β-lines). One thing to note from Fig. 18(a) is that the characteristic curves (the α-lines) on the left and the right sides of the PVC sheet coincides with the specimen edge, where the PVC sheet specimen edge is shown as dashed line. This further confirms that the boundary condition σ2 = 0 along the lateral edges of the sheet specimen.
With all the ingredients ready, we follow the steps illustrated in Fig. 8 and use Eqs. (44)–(46) to determine the two principal stresses σ1 and σ2 in the PVC sheet. The computed principal stress fields (σ1, σ2) in PVC sheet specimen at moment (D) are presented in Fig. 19 for three different mesh sizes. We observe no significant differences when the mesh size decreases or the number of data point increases in the computation.
Finally, we would like to point out, again, that in determining the principal stress fields, σ1 and σ2, as shown in Fig. 19, the only information required is the principal stress orientation field θσ shown in Fig. 17(a), together with the average stretch over the top and the bottom edges of the gage section. The principal stress orientation field θσ and the average stretch are the results from the deformation measurement. The principal stress fields, σ1 and σ2, are then obtained by solving the equation of equilibrium over the domain shown in Fig. 18(a), with traction boundary conditions. Except the assumption of isotropy, no other knowledge of constitutive behavior of the deforming material gets involved.
7.3 Comparison With Constitutive Model Predictions.
We have mentioned in Sec. 1 that the usual approach for obtaining stresses associated with a nonuniform deformation field is to use a constitutive relation. Such a constitutive model needs to be calibrated for a particular material through independent experiment involving only simple loading and deformation configurations, or the stress state of the deforming sample in the calibration experiment is simple and can be calculated from applied force. In this section, we will follow such an approach to obtain the stresses associated with the deformation field of the PVC sheet shown in Figs. 13 and 15.
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.
With the calibrated parameters μ, b, and n for the PVC sheet, we can use Eq. (53) to predict the principal stresses in the case of PVC sheet with a circular hole and subject to uniaxial tension, since the principal stretch fields λ1 and λ2 have been experimentally measured, as shown in Fig. 15. The predicted stress fields σ1 and σ2, based on the GNH model, are presented in Fig. 21 for the selected moments indicated in Fig. 12(b).
Corresponding to those stress fields predicted from the GNH model, shown in Fig. 21, we also determine the two principal stresses, σ1 and σ2, using the technique proposed in the present study and based on the measurement shown in Fig. 16. These stress fields are shown in Fig. 22 for the selected moments indicated in Fig. 12(b).
We first note that the domain over which the stress fields are determined in Fig. 22 is different from that shown in Fig. 21. The domain, over which the stresses are given in Fig. 21, is the current or deformed domain. Figure 23(a) shows the sample domain in the undeformed or reference state, indicated as . The boundaries or edges of domain 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 , 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 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 . This observation highlights the necessity of careful design of both the specimen and the loading fixture when doing the experiment, so that the computational domain can be as close to the current domain as possible, thus minimizes the degree of approximation of the traction boundary conditions.
The computed stress field in Fig. 22 and the predicted stress field in Fig. 21 have similar structures and slightly different overall magnitude, especially for the component of the major principal stress σ1, where the predicted stress from the GNH model is higher than that of computed stress. Such a difference may be the result of the approximation of the traction boundary condition we just discussed above. The matching of the minor principal stress σ2 between the GNH model prediction and the computed stress is much better. On the other hand, we also see that the predicted stress from the GNH model does not satisfy the traction-free boundary conditions along the lateral edges and along the edge of the cavity. Along the lateral edges of the PVC sheet, σ2 is small but not zero as the boundary condition requires. Along the boundary of the cavity, σ2 obviously does not vanish. It is likely that the stress field predicted from the GNH model does not satisfy the equation of equilibrium as well.
Finally, we want to point out an interesting observation that in the computed stress field, stress concentration does not occur at the cavity edges, but slightly inward at location near the cavity boundary. In contrast, stress concentration occurs right at the cavity edge according to the stress field predicted from the GNH model, which is consistent with intuition and many closed-form solutions. The fact that stress concentration does not occur at the cavity edge in Fig. 22 might be due to the artifact of either the numerical scheme or the experiment itself. We note that the PVC sheet, or any polymer sheet, has very low bending stiffness and may be easily bent during a planar deformation causing out-of-plane motion, especially along edges of the specimen. Such out-of-plane motion, though may be exceedingly small, can influence the in-plane deformation measurement, and since the stress determination scheme proposed in this study relies on the detailed deformation field measurement, any error in deformation measurement will manifest itself in the resulting stress field.
7.4 Further Comments on the Traction Boundary Condition and Its Approximation.
In the experimental case of a PVC sheet with a circular hole and subject to uniaxial tension, discussed in this section, the imposed boundary conditions over the top and the bottom of the gage section are actually displacement boundary conditions. As demonstrated in Fig. 13, the lateral displacement u1 along the top and the bottom edge of the gage section is approximately zero and the vertical displacement u2 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.
Furthermore, we assume that f1(sα, sβ) and fN(sα, sβ) can be represented as the sum of a series of known functions, ξi(sα, sβ) (i = 1, 2, …, K), like
One can see that the approach sketched here for obtaining more accurate traction boundary conditions along the characteristic lines, and , is much more numerically involved. However, we have converted the stress determination problem into an optimization problem of minimizing the energetic measure Π, where the stresses are just the by-product of the process, as a result, more powerful numerical techniques, like machine learning and neural networking, can be applied.
This approach may go even further. We may directly consider the displacement boundary conditions as shown in Fig. 13. Rather than limited to the characteristic curve boundary of the computational domain, we will try to determine the tractions along the top and the bottom edge of the current domain shown in Fig. 23, where the displacement is measured experimentally. At any point along this displacement-prescribed boundary, we need to specify the value of a combination of both σ1 and σ2, as indicated in Eq. (25), which can be expressed as a sum of the expansion of known functions with yet-to-be-determined coefficients. Through the optimization or minimization of the energetic measure Π, both the tractions along the displacement-prescribed boundary and the stresses within the current domain can be obtained.
8 Concluding Remarks
In this study, we explored and demonstrated a technique that will determine the stresses associated with a nonuniform deformation field without knowing the detailed constitutive behavior of the deforming material. Here, we do not attempt to deal with the complete boundary value problem in continuum mechanics, but rather we focus on some special circumstances, where if the detailed deformation field over a domain has been measured and traction boundary condition can be identified or approximated, the equation of equilibrium can be solved numerically for the two principal stresses. Since the technique relies on the full-field experimental measurement of the deformation field, in particular, the field of the principal stretch orientation over the whole domain, we are limited to planar or two-dimensional (2D) field only. Meanwhile, the technique also relies on the observation that the Lagrangian strain and the second Piola–Kirchhoff stress for isotropic materials share the same eigenvectors; therefore, the technique proposed in the present study applies to isotropic solids only.
The experimentally measured deformation field over the whole domain plays the key role in the proposed technique. The measured deformation field not only provides the principal stretch orientation θλ at every location but also provides the rigid-body rotation θq of the material particle at that location. The principal stress orientation θσ is the sum of these two angles, according to Eq. (22). With the principal stress orientation θσ field specified, the equation of equilibrium becomes numerically solvable with only the traction boundary conditions. Thus, the numerical computation of the stresses closely intertwines with the experimental measurement.
The close linkage of experimental measurement and numerical computation also poses high demand regarding the accuracy and signal/noise ratios of the experimental technique used for deformation field measurements. Any experimental measurement results contain uncertainty and noise. How sensitive of the proposed numerical scheme on the measurement uncertainty and noise has not been addressed in the present study. It is believed that the optical technique DIC used in the present investigation provides very high accuracy in measuring the displacement field, especially for finite or large deformations. The close match of the lateral edges of the PVC sheet, identified from the DIC result, and the characteristic curves, which are the results of solving the ODEs based on the measured θσ field, is a clear indication of such accuracy of the DIC technique.
The numerical scheme proposed and used in this investigation for solving the equation of equilibrium and computing the stresses depends on the characteristic curves. Once the principal stress orientation θσ field is specified, the two family of the characteristic curves are completely determined. The role of characteristic curves in solving first-order PDEs is quite important. For a single PDE involving a single unknown function, if the PDE can be recast onto the characteristic line, it becomes an ODE and is readily solvable. The numerical scheme proposed in the present study uses similar idea. The traction boundary conditions are described over the characteristic curves that form the boundary of the computational domain, as discussed in Sec. 7.3. Some of these traction boundary conditions are exact and some are approximate, and, in addition to a better sample and loading fixture design, we have given a rough sketch in Sec. 7.4 for improving the accuracy of the tractions over some of the characteristic curves. Such improvement will be explored in a future study.
With the deformation field measured from experiment and complemented with the stresses determined using the technique and the numerical scheme discussed in this study, we have a useful tool for examining different constitutive models on different materials.
The technique discussed in this paper could be extended in different directions and we briefly discuss a few of them here. The first is the situation where domain 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 , so we should be able to determine the stresses within . 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 . 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.
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.
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.