By Dr. Barna Szabó
Engineering Software Research and Development, Inc.
St. Louis, Missouri USA

In mechanical and structural engineering, strength criteria are commonly expressed as allowable limits on scalar measures derived from the stress tensor field in a linear elasticity problem. These limits are determined from experimental observations and engineering safety requirements.

We must remember that continuum mechanics assumes materials have no microstructure and therefore treats them as capable of subdivision without limit. Real materials do not satisfy this assumption. Consequently, pointwise stresses computed from a continuum model should be interpreted as averages over a representative volume element (RVE).

As long as the stress field varies smoothly, such that the stress changes only slightly across an RVE, continuum theory may be understood as providing RVE-averaged stresses.  In the following, we consider situations in which the averaging process must be made explicit.

Setting the Goals of Computation

In regions where the stress field varies significantly over distances comparable to the characteristic material length scale, we must account for the nonlocal nature of the quantity governing failure initiation. This observation has motivated a substantial body of research on nonlocal approaches to strength assessment, including the Theory of Critical Distances, line- and volume-averaging methods, stress-field intensity approaches, and highly stressed volume methods. Reviews of these approaches can be found in [1,2].

The common theme underlying nonlocal approaches is that failure initiation can be correlated with a scalar-valued quantity derived from the elastic stress field. Although failure initiation is inherently nonlinear and irreversible, the process zone is surrounded by material that remains essentially elastic. The nonlinear aspects of the failure process are not modeled explicitly but are instead incorporated through calibration.

Therefore, all nonlocal approaches require the accurate computation of the quantity of interest (QoI). Without rigorous control of numerical approximation errors, the calibration process would tend to compensate for discretization errors, thereby compromising the predictive performance of the resulting model.

Strong and Generalized Formulations

In linear elasticity, we assume strains are much smaller than unity and that the stress tensor is proportional to the strain tensor. From the equilibrium equations and the strain-displacement relations, we obtain a second-order partial differential equation in terms of the displacement field. Therefore, this formulation, which belongs in the category of strong formulations, is meaningful only when the second derivatives of the displacement field are finite. 

The finite element method is based on a generalized formulation, namely the principle of virtual work, which requires only that the first derivatives of the displacement field be square-integrable. Consequently, the formulation encompasses a much larger class of admissible solutions than the strong formulation [3].  This observation is particularly relevant to stress-based failure assessment because failure criteria are often evaluated at locations where singular or nearly singular stress fields occur. For example, solutions of the form:

\mathbf{u} = r^\lambda \boldsymbol{\Phi}(\theta), \qquad \lambda \geq 0 \qquad \qquad (1)

are admissible. Here, u is the displacement field, r and θ are polar coordinates and Φ is a smooth vector function. Numerical estimation of λ­­­ and Φ is discussed in [5]. 

When λ < 1, which may occur at material interfaces, then the strain and hence stress are infinite at r = 0. Consequently, the generalized formulation admits solutions that are mathematically admissible yet incompatible with the fundamental assumption of linear elasticity that strains remain small.  The engineering challenge is to interpret such results correctly and determine the extent to which they provide useful information. 

A common mistake is to report stresses when the underlying exact solution has such singularities. The maximum stress computed from a finite element solution will be finite for any discretization, but will not tend to a finite value as the number of degrees of freedom increases.

Example

We consider the test specimen shown in Figure 1(a). It consists of a plate measuring 100 mm in length and 40 mm in width. A circular hole of diameter 10 mm is located at the center of the plate. The material is an IM7/8551-7A prepreg with a [0₁₀/90₁₀/0₁₀] layup. Each lamina is 0.13 mm thick. The material properties are listed in Table 1. The specimen is loaded by uniform tension of 100 MPa in the x-direction. The domain of primary interest is the neighborhood of the hole, shown in Figure 1(b).

Figure 1: (a) Finite element mesh (156 elements). (b) Hashin’s matrix tension failure index.

The finite element mesh shown in Figure 1(a) consists of 156 hexahedral elements. This mesh design is typical for problems where boundary layer effects may have to be considered. The mapping functions were fifth-degree polynomials. The computations were performed using StressCheck. The goal is to demonstrate an application of a failure criterion based on the linear elastic stress field. We will use Hashin’s Matrix Tension Failure Index [4] for this purpose:

F_{mt} = \left(\frac{\sigma_{22} + \sigma_{33}}{Y_t}\right)^2 + \frac{\tau_{12}^2 + \tau_{13}^2}{S_\ell^2} + \frac{\tau_{23}^2 - \sigma_{22}\sigma_{33}}{S_t^2}, \qquad \sigma_{22} + \sigma_{33} \geq 0 \qquad \qquad (2)

where index 1 indicates the direction parallel to the fiber index, 3 is normal to the mid-surface of the lamina, and index 2 is normal to the 1-3 plane. The symbol Y_t represents the transverse tensile strength, S_l and S_t are the longitudinal (in-plane) and transverse shear strengths, respectively. Here we will use Y_t = 75.890 MPa, S_l = 90 MPa and S_t = 80 MPa. Matrix tension failure is predicted when F_mt ≥ 1. 

Table 1: Elastic material properties

Figure 2 shows the computed pointwise (local) values of F_mt in the region of primary interest for polynomial degrees ranging from 3 to 8 as a function of the number of degrees of freedom N. We used the finite element mesh shown in Fig. 1. Since the values increase monotonically with the polynomial degree, the computations were repeated on a refined mesh to assess convergence. Refinement was achieved by subdividing each hexahedron into eight hexahedra, resulting in 1248 elements. The results corresponding to the refined-mesh ­­­­­indicate convergence toward a limiting value of approximately 1.03.

Figure 2: Convergence curves for the local and nonlocal maxima of Hashin’s matrix tension failure index.

Referring to Fig. 1(b), we observe that high values of F_mt occur within a small volume bounded by a portion of the hole’s inner surface, extending through the thickness of the 0° lamina and over a short radial distance. Formally, we define this domain in cylindrical coordinates as follows:

\Omega_c = \left\{ r,\theta,z \;\middle|\; r_h < r < r_h + \Delta r,\; \theta_1 < \theta < \theta_2,\; z_{k-1} < z < z_{k-1} + t_{\mathrm{ply}} \right\} \qquad \qquad (3)

where r_h is the radius of the hole, Δr is a distance, defining the radial extent of the domain, k is the ply number and t_ply is the ply thickness. The values of θ₁ and θ₂ are determined from

\theta_i = \theta_0 + (-1)^i \frac{s}{2r_h}, \qquad i = 1,2 \qquad \qquad \qquad (4)

where θ₀ is the polar angle identifying the location of the maximum value of F_mt and s is an experimentally determined arc length. For the present problem, a reasonable initial estimate of s is s = 10t_ply, that is, the subtended angle is 10t_ply/r_h. To demonstrate the process, we chose Δr = 0.25 mm and s = 1.3 mm for this example.

Figure 2 shows the nonlocal (volume-averaged) values for the 156-element and 1248-element meshes. While the sequence based on the 156-element mesh shows no evidence of convergence to a limiting value, the sequence based on the 1248-element mesh converges clearly. Note that the computed value for the 156-element mesh at p = 8 is very close to the limiting value of the 1248-element sequence. However, this is not the case for the local values, indicating that nonlocal values converge faster than their local counterparts. Solution verification requires demonstrating that convergence has occurred.

In this example, the local maximum of F_mt converges to a finite value. However, even in cases where the local maximum does have a finite limit, the nonlocal values will converge to a well-defined limit.

The parameters Δr and s must be inferred from experimental data. They are crucial for establishing a predictive model based on linear elasticity. Proper inference requires that the exact QoIs are finite and that numerical errors are negligible compared to experimental errors.

Here, we computed nonlocal values of F_mt by volume averaging. An alternative approach is to estimate the nonlocal values of F_mt from volume-averaged stress components.

Key Takeaways

• The development of decision-grade predictive models requires both experimental calibration and rigorous control of numerical errors in the quantities of interest.
• Failure initiation criteria, including the Hashin, Tsai-Wu, Tsai-Hill, Puck, and LaRC formulations, are expressed in terms of local stress or strain quantities evaluated at material points. Applying these criteria to real materials requires the consideration of two cases. In the first case, the stress varies smoothly, exhibiting negligible variation over the size of an RVE. Consequently, the pointwise stress is virtually identical to the average stress over the RVE. In the second case, the stress gradient is steep, necessitating explicit stress averaging.  This case is of particular interest because such regions are the usual sites of failure initiation.
• Explicit stress averaging requires the introduction of length parameters which must be inferred from experimental observations.  Proper inference requires that the numerical errors in the QoI be negligibly small in comparison with experimental errors.  Therefore, solution verification is an essential technical requirement. The hp-version of the finite element method is ideally suited for this purpose. 
• The development and ranking of competing predictors of failure initiation must follow a well-defined, properly executed experimental program, as described in [6]. Such experimental investigations should be accompanied by the development of mathematical models capable of generalizing the resulting data beyond the specific conditions under which they were obtained.

References

[1] Moghtaderi, S.H., Jedi, A. and Ariffin, A.K. (2023). A review on nonlocal theories in fatigue assessment of solids. Materials, 16(2), p. 831.

[2] Delprete, C., Maggio, L.G.D. and Sesana, R. (2021). Theory of critical distances: A discussion on concepts and applications. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 235(21), pp. 5695-5708.

[3] Szabó, B. and Babuška, I. (2021).  Finite Element Analysis. Method, Verification and Validation. 2^nd Edition. John Wiley & Sons Inc.  New York.

[4] Hashin, Z. (1980). Failure criteria for unidirectional fiber composites. Journal of Applied Mechanics, 47(2), pp. 329–334.

[5] Yosibash, Z. (1997). Computing edge singularities in elastic anisotropic three-dimensional domains. International Journal of Fracture, 86(3), pp. 221-245.

[6] Szabó, B. and Actis, R. The demarcation problem in the applied sciences. Computers & Mathematics with Applications, 162, pp. 206-214, 2024.

Related Blog Posts:

• A Memo from the 5th Century BC
• Obstacles to Progress
• Why Finite Element Modeling is Not Numerical Simulation?
• XAI Will Force Clear Thinking About the Nature of Mathematical Models
• The Story of the P-version in a Nutshell
• Why Worry About Singularities?
• Questions About Singularities
• A Low-Hanging Fruit: Smart Engineering Simulation Applications
• The Demarcation Problem in the Engineering Sciences
• Model Development in the Engineering Sciences
• Certification by Analysis (CbA) – Are We There Yet?
• Not All Models Are Wrong
• Digital Twins
• Digital Transformation
• Simulation Governance
• Variational Crimes
• The Kuhn Cycle in the Engineering Sciences
• Finite Element Libraries: Mixing the “What” with the “How”
• A Critique of the World Wide Failure Exercise
• Meshless Methods
• Isogeometric Analysis (IGA)
• Chaos in the Brickyard Revisited
• Why Is Solution Verification Necessary?
• Variational Crimes and Refloating the Costa Concordia
• Lessons From a Failed Model Development Project
• Where Do You Get the Courage to Sign the Blueprint?
• Great Expectations: Agentic AI in Mechanical Engineering
• The Differences Between Calibration and Tuning
• Honored in the Breach
• Remembering Ivo Babuška
• Turtle Shells and Legacy Finite Element Codes: Evolutionary Constraints in the Age of Explainable AI
• Beyond the Black Box: Explainable AI Requires Explainable Simulation
• Trustworthiness in Simulation: Credibility or Decision-grade Reliability?
• From Finite Element Modeling to Simulation Science: Unlocking the Power of Numerical Simulation
• Why is the hp-version of the Finite Element Method Ideally Suited for Laminated Composites?