Stress Analysis and Stiffness Degradation of Open Cracks Composite Laminates Subjected to External Loads

TL;DR

This paper studies how open cracks in composite laminates affect stress distribution and stiffness under external loads.

Contribution

A new method is proposed to analyze stress distribution and stiffness degradation in cracked composite laminates using equilibrium differential equations and complementary energy minimization.

Findings

Fiber orientation and material properties significantly affect stress distribution in symmetric composite laminates.

Structural parameters influence the reduction of cross-layer shear modulus in laminates with transverse matrix cracks.

The results have practical implications for engineering applications involving composite materials.

Abstract

Composite laminated structures have extensive applications in the field of bionic engineering. Proficient comprehension of the mechanical properties of these structures is instrumental in the advancement of bionic composite materials. The objective of this study is to investigate the stress distribution and degradation of stiffness in composite laminates exhibiting open smooth surface cracks under varying external loads and structural parameters. Utilizing the general series function of the laminate’s axial stress, the general expression for the stress components of the damaged laminate is derived by integrating the equilibrium differential equation, boundary conditions, and stress continuity conditions. The influence of fiber orientation and material properties on the stress distribution within each layer of symmetric composite laminates was examined. Thereafter, the reduction in cross-layer shear modulus was assessed by employing the principle of complementary energy minimization. The impact of structural parameters on shear modulus reduction was explored. The findings indicate that structural and material parameters of symmetric laminates featuring transverse matrix cracks exert a notable influence on the stress distribution and degradation of stiffness within each layer, imparting practical significance to the research outcomes in engineering applications.

Figures

Full text

1. Introduction

Bionics, also referred to as biomimicry or heuristics, entails the study of resolving artificial intelligence and technological challenges by emulating the characteristics and processes of biological systems [1,2,3]. This interdisciplinary field integrates knowledge and methodologies from various disciplines, including biology, physics, chemistry, mathematics, and computer science. The research in bionics is comprehensive, encompassing areas such as biomimetic robotics, biomimetic materials, and biomimetic algorithms. It represents the cutting edge in exploring the integration of natural intelligence with artificial constructs, propelling materials science toward new frontiers. In this domain, there is a growing interest in composite laminated structures due to their novel design concepts and practical utility. Composite sandwich structures, in particular, find diverse applications in bionic engineering, such as the development of bionic shrimp shell structures, bionic glass sandwich panels, and impact-resistant bionic mantis shrimp sandwich boards. However, the use of composite materials may result in internal cracks and damage within the laminated structures, which can severely compromise their mechanical properties. Consequently, numerous scholars have conducted extensive research on these layered structures. In general angle-laid laminates, the initial form of damage incurred under external loading is the occurrence of transverse matrix cracks in the 90° layer [4,5,6,7,8,9]. Such cracks are likely to decrease the stiffness of the laminates and pose a significant damage potential. Hajikazemi [10] proposed a more accurate method to calculate the effective thermal expansion coefficient of the cracked laminates and combined with the allowable stress field, the degradation of the thermoelastic constant of the laminates with transverse cracks could be predicted more accurately and then compared with the calculated results of various models. Finally, it is found that the calculated results of this method are more accurate than those of the shear hysteresis model. In the stress-based variational model established by Huang [11], the normal stress of two adjacent layers of the crack layer is introduced into the statically indeterminate boundary problem. Then, by using the principle of minimum complementary energy, the stress components of each layer can be derived by the variational method and solved by the governing equation. Then, according to the classical elastic theory of laminates, Huang [12] proposed that there are two types of transverse matrix cracks in angle-laid laminates, namely, open cracks and closed cracks. After variational analysis of these two kinds of cracks, it is found that there are great differences in stress distribution and effective mechanical properties between the open crack and the closed crack laminates. Huang [13] found that most current studies have ignored the influence of fiber bridging resistance between irregularly cracked surfaces on the effective mechanical properties of damaged laminates, and their research results are limited to situations where the cracked surface is smooth. In 2010, Cortes and Barbero [14] applied their method to any stacking order of symmetric laminates with diagonal cross-matrix cracks and matrix cracks in one or two directions, showing good prediction ability for crack initiation and evolution in different laminates. In 2014, Carraro [15] studied the stiffness degradation of multidirectional symmetric laminates under the action of off-axis cracks. A new analytical model is proposed, which can be used to predict the elastic properties of multidirectional symmetric laminates with arbitrary layering and one or more layers of cracks. The model is based on the so-called “optimal” shear hysteresis analysis of single-layer (symmetrically) cracked laminates and an analytical procedure that explains the presence of multilayer cracks while taking into account their interactions. The predicted results are in good agreement with the experimental data in the literature. The calculation results also show that considering the interaction of cracks between different layers is the basis for estimating stiffness changes reliably. Liang Hanbing [16] set the interlaminar shear stress in the form of a general function and obtained the stress component expression of the opening crack from the equilibrium equation and boundary conditions. On this basis, the governing equations of the opening crack were obtained by using the principle of minimum complementary energy, and then the stress distribution of the symmetric laminates of composite materials was finally obtained by analytical analysis and solution. The effect of fiber bridging on the shear modulus of composite laminates is analyzed, and the effect of bridging factors on the effective shear modulus is calculated. Huang [17] then obtained an exact solution of the stress in the cracked laminate, taking into account the series expansion form of the sine function of the shear stress between the 90° layer and the 0° layer. Meanwhile, the results of the interlayer stress show that the delamination occurs near the splitting position. Duan and Gan [18] proposed a semi-analytic finite element method to simulate composite laminates with arbitrary fiber orientation and anisotropic material properties in each layer. Hajikazemi and Sadr [19] established a variational model to analyze the stress field of symmetrically fractured laminates under general loads in the plane. Thus, the applicability of variational methods in considering multilayered symmetric laminates with arbitrary superimposed sequences is fundamentally extended. They compare the results of stress distribution in the literature for the case of orthogonal laminates. In this paper, Hajikazemmi [10] tries the energy variational principle for crack analysis of symmetric laminates with arbitrary lamination order under general plane loading. Some scholars, such as Zhang [20,21,22], assumed that the shear stress in laminates is linear through the entire thickness layer in the 90° layer, while the shear stress is equal to zero in the 0° layer. In view of the energy variational principle, Nairn [23], Varna and Berglund [24,25], and Berglund and Varna [26] further refined the model by assuming that in the 90° layer, the shear stress $[eqn]$ is a linear function of z and the stress $[eqn]$ is a quadratic function of z. The 0° layer may have a non-uniform stress distribution described by an exponential function with unknown shape parameters. The problem is reduced to a fourth-order differential equation with constant coefficients to find the minimum value of the additional energy. The equation and its solution contain an unknown shape parameter, which is calculated in the subsequent additional energy minimization. Berglund’s improvement introduced the statically indeterminate boundary problem. In 2014, Carraro [15] used an analytical method to provide a computational model for the elastic properties of any number of layered symmetric laminates, taking into account their interactions. For this reason, the two-dimensional optimal shear hysteresis analysis of laminates with single cracks was carried out for the first time. Gabbert et al. [7] studied unidirectional fiber composites with unsatisfactory interface conditions between reinforcement and filler. The microstructure is periodic and the phase is isotropic. The periodicity of the microstructure appears as a parallelogram. Using the concept of representative volume element (RVE), a finite element model of fiber distribution and incomplete contact conditions between phase interfaces is established. June M and Kayran A [27] conducted a finite element analysis of the geometric nonlinear progressive failure of perforated composite laminates under both in-plane and out-of-plane composite loads. Different failure criteria and material property degradation schemes based on layers and components were written, and a progressive failure analysis program was developed. Based on the finite element method, Vejen and Pyrz [28] studied the transverse crack propagation of long fiber composites. Three criteria of pure matrix growth, fiber/matrix interface growth, and fiber/matrix crack bending outside the interface are realized, and a crack propagation module of a software package is formed. The numerical crack path is compared with the experimental crack path. Finally, the crack propagation module is used to construct an example to obtain information about the influence of various fiber distributions on the crack propagation process. Ravi Joshi [29] applies evenly distributed transverse loads to the antisymmetric and symmetric composite laminates and conducts failure analysis for each layer in turn. Hashin [30] analyzed 90 of them using a variational method. The tensile and shear loads of the laminates with transverse matrix cracks are contained in the laminates. The allowable stress system satisfying the equilibrium condition and all boundary and interface conditions is constructed, and the best approximation is found by using the principle of minimum complementary energy. Xu Dong [31] not only considers three ideal crack conditions of smooth plane opening crack, smooth plane closing crack, and fully meshing opening crack, but also analyzes two conditions of non-ideal crack, that is, non-ideal opening crack and non-ideal closing crack. The study not only verifies the symmetry of the flexibility element and the equality of the change rate of the shear flexibility element under closed conditions but also shows that the non-ideal degree of crack has a great influence on the effective coupling coefficient and the effective shear flexibility. Giannadakis [32] applies a more detailed and accurate method to estimate stress using shape function dependencies with parameters minimized. Then, the results of coenergy are compared with the corresponding results in the literature, and finally the expression of shear modulus degradation is derived. An extensive review of the literature reveals that previous research has concentrated primarily on stress distribution within individual composites or solely on the 90° cracked layer, with a specific focus on the alterations in positive (tangential) stress distribution across the plane and between layers at varying crack densities (CD). Here, CD refers to the number of cracks per unit length. Limited studies have examined the influence of fiber orientation—defined as the angle between fiber direction and the laminate’s x-axis—on the interlaminar stress distribution’s variability, as well as the impact of diverse material characteristics on the internal stress distribution in composite laminates with exposed cracks. This paper introduces a stress distribution analysis model for fiber-reinforced composite laminates featuring smooth surface cracks, characterized by the absence of contact and shear force at the crack interface. An integral contribution of this study is the derivation of a general expression for the internal stress components in interlaminar plates, achieved through the employment of the variational method and energy principles. These expressions satisfy the criteria of elastomer equilibrium differential equations, stress boundary conditions, and stress continuity. Initially, the paper examines the impact of fiber orientation on the stress distribution within a single composite layer with smooth cracks subjected to two distinct external loads—tensile and shear—and elucidates the mechanisms behind in-plane tangential and interlaminar shear stress distribution. Subsequent sections provide an in-depth analysis of various composite materials, including glass/epoxy 1, glass/epoxy 2, and carbon fiber/epoxy, with a focus on stress distribution under different material properties. The paper concludes with an investigation into the influence of matrix crack propagation on the shear modulus of composite laminates, utilizing the principle of complementary energy minimization and accounting for different structural parameters, such as the number of 0° single layers, thickness, and the thickness of crack layers. These investigations comprehensively elucidate the effects of fiber orientation and material properties on the stress distribution and shear modulus following matrix crack diffusion. This research advances our understanding of how crack diffusion affects the mechanical properties of composite laminates, fostering the development of mechanics and stiffness degradation theories for these materials, and providing a robust theoretical foundation for the engineering applications of fiber-reinforced composites.

2. Boundary Conditions and Stress Components

2.1. Element Volume Model

The early experimental results [4,5,6,7,8,9] demonstrate that as the crack density increases, transverse cracks within the 90° layer of the matrix-cracked laminates are uniformly distributed along the length of the plate, extending across the full width and throughout the entire thickness. Tensile stress $[eqn]$ is exerted on the composite laminates $[eqn]$ with laterally cracked matrix. As shown in Figure 1a,b, the thickness of the 90° layer and $[eqn]$ layer is $[eqn]$ . $[eqn]$ is the thickness of the single-layer plate. Subscript 1 represents the 90° layer, and subscript 2 represents $[eqn]$ layer. The distance between the two cracks is $[eqn]$ . Assuming that the cracks are equally spaced, the crack density is CD = $[eqn]$ . The total thickness of the laminate is $[eqn]$ . The whole layer is laid in accordance with $[eqn]$ , 90°, $[eqn]$ , and the structure of the laminate has upper and lower symmetry. Therefore, only the upper part of the volume element needs to be studied.

2.2. No Thickness Spring Model

Based on the actual defects present in the material, the non-ideal matrix crack is conceptualized as a defect interface. Taking into account the continuity of stress across the defect and the discontinuity of displacement, the theoretical model of spring without thickness is established.

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ represent imperfect interface parameters, and expressed in units of Pa/m. $[eqn]$ and $[eqn]$ denote the differences in normal and tangential displacements, respectively, between the two opposite surfaces of the defect. The values of the two independent imperfect interface parameters directly indicate the extent of ideal versus non-ideal defects. ① When $[eqn]$ , $[eqn]$ , the normal and tangential displacements are discontinuous, and the corresponding normal and tangential stresses are zero, indicative of a smooth plane opening crack. ② When $[eqn]$ , $[eqn]$ , the normal displacement is continuous, the normal displacement remains continuous, the normal stress is non-zero, the tangential displacement is discontinuous, the tangential stress is zero, and a smooth plane closed crack can be characterized. ③ When $[eqn]$ , $[eqn]$ , the normal displacement is discontinuous, the normal stress is zero, the tangential displacement is continuous, the tangential stress is non-zero, and a fully meshed opening crack can be characterized. (It should be noted that if the mesh represents a closed crack, the normal and tangential displacements are continuous, and the corresponding normal and tangential stresses are non-zero, which may indicate no damage). ④ When $[eqn]$ , $[eqn]$ , the normal and tangential displacements are continuous, and the corresponding normal and tangential stresses are non-zero, indicative of a non-destructive injury scenario.

2.3. Boundary Conditions

For the statically indeterminate boundary, two undetermined zonal loads are introduced in accordance with tensile loading. Subsequently, the directional stress and tangential stress within regions (1) and (2) at the RVE boundary, under various loading conditions, must satisfy the following criteria:

For tensile loading, $[eqn]$ and $[eqn]$ are introduced:

$$[eqn]$$

Therefore, the boundary conditions are expressed as follows [31]:

$$[eqn]$$

In symmetrically laminated plates, the external shear stress at the middle surface (Z = 0) of layer (1) is zero, viz:

$$[eqn]$$

On the surface of the laminate, at the top of layer (2) (Z = h), the surface shear stress and normal stress are 0, that is,

$$[eqn]$$

2.4. The Stress Components Within Each Layer

The load is uniformly distributed along the Y-axis. Hence, the stress component may be presumed to be a function of the variables x and z. This allows for the additional simplification that the out-of-plane shear stress is linearly dependent on the thickness of each layer, while the in-plane stress within each layer is independent of z. Consequently, under tensile loading, the interlaminar stress between region (1) and the adjacent region (2) is represented by Equation (7) as follows:

$$[eqn]$$

where in the superscript (12) denotes the interface between layers (1) and (2) and $[eqn]$ is a dimensionless variable. The combined formulas from (3) to (6), along with the interlayer stress conditions, must remain continuous. These can be individually derived for the stress components within each zone under tensile loading.

For the 90° layer,

$$[eqn]$$

For the $[eqn]$ layer,

$$[eqn]$$

Taking into account the continuity of stress, Equations (8) and (9), when combined with the boundary conditions, yield the following:

$$[eqn]$$

Within the formula $[eqn]$ , $[eqn]$ and $[eqn]$ represent three independent undetermined functions that are employed to depict the interlayer stress $[eqn]$ , and $[eqn]$ denotes undetermined load distribution coefficients. Each of these unknown functions and coefficients can be resolved using the governing equation.

3. Governing Equation and Supplementary Equation

3.1. Establishing Governing Equations and Supplementary Equations

For open cracks:

In adherence to the principle of minimum complementary energy, the governing equations for the functions $[eqn]$ , $[eqn]$ and $[eqn]$ can be established through variational analysis and are expressed as follows:

$$[eqn]$$

where U is the complementary energy given by

$$[eqn]$$

in which

$$[eqn]$$ $$[eqn]$$

and $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ are the normal and shear strain components, respectively.

Substituting Equations (8) and (9) into Equation (13) and applying Hooke’s law and Equation (11) yield the non-zero governing equation, given by $[eqn]$ , $[eqn]$ and $[eqn]$ :

$$[eqn]$$

The coefficients ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) and the constants $[eqn]$ are listed in Appendix A.

The governing equations for the functions $[eqn]$ and $[eqn]$ are obtained as

$$[eqn]$$

The function $[eqn]$ can be obtained as

$$[eqn]$$

The coefficients ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) and constants $[eqn]$ are listed in Appendix B.

3.2. Solution to Governing and Supplementary Equations

The equations can be solved using linear differential equation theory. Formula (15) is the sum of the characteristic constants of the general solution of the corresponding homogeneous equation and the functions $[eqn]$ and $[eqn]$ of the nonhomogeneous equations. The general solution of the homogeneous equations can be written by the root $[eqn]$ (j = 1, 2, …, 6) of the characteristic equation, such that

$$[eqn]$$

and the eigenvector $[eqn]$ can be obtained from the following equation:

$$[eqn]$$

For the distinct roots in our following numerical examples, some of the characteristic roots may be the same, the general solutions to the homogeneous equations cannot be expressed as exponential function $[eqn]$ s, and the forms of the solutions need to be slightly modified to accommodate them; however, these procedures are all standard treatments in ordinary differential equations. Because the algebraic equation of order is an even order equation of r, the six roots of the equation correspond to three pairs of positive and negative roots $[eqn]$ (j = 1, 2, 3). Therefore, the solution to a system of nonhomogeneous linear equations can be written as

$$[eqn]$$

For different roots, where $[eqn]$ and $[eqn]$ (j = 1, 2, 3) are unknown constants $[eqn]$ , and $[eqn]$ are undetermined constants. Using the boundary conditions (10), we have $[eqn]$ = $[eqn]$ . The solution (19) can easily be rewritten as

$$[eqn]$$

where $[eqn]$ .

For open cracks, in order to determine the constants $[eqn]$ (j = 1, 2, 3), particular solutions $[eqn]$ and $[eqn]$ can be obtained by combining the following homogeneous equations:

$$[eqn]$$

4. Stiffness Degradation Research

This section investigates the reduction in the shear modulus of orthogonal laminates with matrix cracks, utilizing a theoretical model based on the minimization of complementary energy. Variations in the shear modulus of the laminates are examined by manipulating the structure and material parameters. As per the classical lamination theory (CLT), the following is observed:

$$[eqn]$$

Equation (23) $[eqn]$ represents the shear modulus of $[eqn]$ layer, whereas $[eqn]$ denotes the shear modulus of the unidirectional (UD) composite. The laminate is assumed to be subjected to the average in-plane shear stress, denoted as $[eqn]$ . The shear strain of the undamaged laminate, denoted as $[eqn]$ , and the shear stresses within the laminate, denoted as $[eqn]$ , and $[eqn]$ are as follows:

$$[eqn]$$

4.1. Assumed Stress Distributions

Two unknown stress disturbance functions, denoted as $[eqn]$ and $[eqn]$ , are employed to characterize the shear stress distribution in the x-direction, and the following expression for shear stress is utilized:

$$[eqn]$$

The shear stress within the 90° layer is required to satisfy the condition of zero traction on the crack surface:

$$[eqn]$$

From the shear force balance

$$[eqn]$$

we obtain

$$[eqn]$$

leading to

$$[eqn]$$

substituting $[eqn]$ in the equilibrium equation

$$[eqn]$$

and integrating it with respect to x and applying stress free condition on the surface z = h we obtain:

$$[eqn]$$

The same equilibrium equation and boundary condition at z = 0 in 90-layer yields

$$[eqn]$$

It can be checked that the continuity of $[eqn]$ at the layer interface is satisfied. By introducing a new unknown function $[eqn]$ as

$$[eqn]$$

the following stress expressions can be written:

$$[eqn]$$

Finally, the unknown functions $[eqn]$ and $[eqn]$ have to satisfy boundary conditions:

$$[eqn]$$

4.2. Minimization of Complementary Energy

In the calculation of elastic properties, thermal stress is considered negligible [32], thereby allowing the thermal term to be omitted from the complementary energy expression. Consequently, the complementary energy is equated to the strain energy. The strain energy is associated with the repeating unit situated between two cracks in 90° layer:

$$[eqn]$$ $$[eqn]$$

Substituting Equation (34) into Equations (36) and (37) and conducting the Z-integral yields the subsequent expression for the complementary energy of the laminate in the event of a rupture in the 90-layer.

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$

The function Φ that minimizes (38) has to satisfy the Euler–Lagrange equation

$$[eqn]$$

which leads to the following equation for Φ(x) determination:

$$[eqn]$$

The solution of (42), which satisfies the boundary condition (35), is

$$[eqn]$$

The solution provided by Equation (43) yields a minimum value for the complementary energy corresponding to the selected value of $[eqn]$ , suggesting that an additional minimization of $[eqn]$ can be attained through a suitable selection. To facilitate this objective, the complementary energy expression (38) can be restructured as follows:

$$[eqn]$$

Replacing $[eqn]$ in the second integral via $[eqn]$ as follows from (42), we obtain

$$[eqn]$$

Substituting in (45), $[eqn]$ given by (39) and expressing $[eqn]$ from (43), we obtain

$$[eqn]$$

Selecting different functions $[eqn]$ , we will have different values of $[eqn]$ and µ, leading to different values $[eqn]$ . The lowest value will correspond to the most accurate stress distribution on average.

4.3. Shear Modulus Expression

The shear modulus of the damaged laminate can be determined using $[eqn]$ . The strain energy density of the RVE is calculated as follows:

$$[eqn]$$

The strain energy density can be obtained by dividing (46) with the RVE volume. Moreover, the shear stress $[eqn]$ in Equation (46) can be represented by the variable $[eqn]$ as per Equation (24), and the strain energy density for undamaged laminates is given by $[eqn]$ . Performing differentiation according to (47), we obtain

$$[eqn]$$

Since Hooke’s law for the damaged laminate has to have the form (24), the following shear modulus expression for the damaged laminate can be obtained:

$$[eqn]$$

For cross-ply laminates, $[eqn]$ and $[eqn]$ . Then, the expressions for the shear modulus can be re-written as

$$[eqn]$$

As the distance from the crack tip (layer interface) increases, it is anticipated that the shear stress disturbance will decay, necessitating a decrease in the absolute value of the function $[eqn]$ and its first derivative. Natural selection exhibiting such characteristics is described by an exponential function. The shape function incorporates a shape parameter n, which serves to differentiate the minimum value of the complementary energy specified by Equation (46). For the present study, the function that satisfies boundary conditions (35) is selected in the following form:

$$[eqn]$$

where n is the number of single-layer plates of 0° layer, e is the natural base number, and z is the thickness coordinate.

Using (18), the parameters $[eqn]$ , $[eqn]$ are calculated as follows:

$$[eqn]$$

5. Validation of Model Feasibility

Given that we are examining smooth open cracks, $[eqn]$ and $[eqn]$ are both equal to zero. Based on Equations (1) and (2), it can be inferred that $[eqn]$ and $[eqn]$ . To validate the accuracy of the variational model, analyses were conducted on laminates made from glass/epoxy resin 1 (or glass/epoxy resin 2). At a crack density of CD = 15 (cr/cm), the distribution of normal and shear stresses within the 90° layer was calculated using MATLAB 2018a software, as depicted in Figure 2a,b. In this paper, the normal stress of the 90° smooth crack layer is compared with the calculated data of Mingkang [8] and Liang Hanbing [16]. The basic theory of composite material mechanics and the energy principle are used in this paper to build the model, which is slightly different from the above references; however, in some aspects, the concrete verification process is used to consider the interface parameter Kt of the non-ideal crack in Mingkang [8]. However, this topic is to consider the smooth crack condition, so that Kt = 10 will degenerate into a smooth surface crack to discuss together, while the non-thickness spring model in this paper when considering the smooth crack condition, $[eqn]$ and $[eqn]$ are zero, can also degenerate into a smooth crack, and the simulation results of Liang Hanbing [16] are basically consistent. The results show that the model has higher precision, the modeling process is more concise, and the efficiency of the calculation results is improved.

6. Stress Distribution Analysis

6.1. Influence of Fiber Laying Angle on Stress Distribution of Laminates with Matrix Cracks

To examine the impact of fiber layering angle on the mechanical properties of symmetric laminates with opening crack damage when subjected to tensile loading, a set of glass/epoxy resin 2 laminates $[eqn]$ were numerically analyzed, and the material constants were derived from Mingkang [8]: $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and a thickness of $[eqn]$ for each lamina. The crack density was set at 0.6 (cr/mm).

Figure 3 illustrates that $[eqn]$ , the normal stress in the x-direction of the 90° layer, is zero at the unit boundary (crack surface), as dictated by the boundary conditions. $[eqn]$ represents the normal stress perpendicular to the fiber direction in the θ-degree layer, decreasing as the θ angle increases due to the reduced alignment of the fiber direction with the x-axis. The distribution of $[eqn]$ is inversely proportional to that of σx1, with $[eqn]$ being significantly greater in magnitude, signifying that the θ-degree layer is the primary load-bearing layer. Conversely, the load-bearing capacity of the 90-degree layer is substantially diminished following cracking. As depicted in Figure 4, $[eqn]$ represents the normal stress in the y-direction of the 90° cracked layer. Its distribution shape is opposite to that of $[eqn]$ , and its magnitude is lower than that of $[eqn]$ , which further indicates that the crack diminishes the bearing capacity. $[eqn]$ denotes the normal stress parallel to the fiber direction, and its magnitude increases with an increasing θ angle. This behavior is in contrast to the decrease in $[eqn]$ with increasing θ angle. In Figure 5, $[eqn]$ and $[eqn]$ represent the shear stresses within the xy plane of the 90° cracked layer and the θ° layer, respectively, with their distributions and shapes being mutually inverse. Due to the presence of the crack surface, the boundary values are zero, with $[eqn]$ exceeding $[eqn]$ in magnitude. The distribution of both normal and shear stresses within the plane is independent of the material thickness, and this distribution remains consistent across any thickness interface.

6.2. Influence of Material Properties on Stress Distribution of Laminates with Matrix Cracks

The present study is designed to investigate the influence of various material parameters on the behavior of damaged laminates, utilizing glass epoxy resin and carbon fiber as the selected materials. The material parameters are derived from the articles of Hashin [30]. Table 1 presents the material constants for two types of epoxy resin and carbon fiber:

The variation in surface stress of laminate $[eqn]$ with material parameters is discussed in the case of open matrix crack.

Figure 6 illustrates that the internal normal stress $[eqn]$ exhibits a descending order from glass/epoxy resin 1 to glass/epoxy resin 2, and finally to carbon fiber/epoxy resin. Moreover, at x = 0, the stress values are maximized for both materials.

Figure 7 demonstrates the in-plane normal stress $[eqn]$ in descending order, with the carbon fiber/epoxy resin exhibiting the highest stress, followed by glass/epoxy resin 2, and then glass/epoxy resin 1. For the carbon fiber/epoxy resin, the stress value increases gradually within the range of −0.8 $[eqn]$ to 0 $[eqn]$ (inclusive of 0.8 $[eqn]$ to 0 $[eqn]$ ), reaching a maximum of approximately 0.4. For glass/epoxy resin 2, the stress value remains stable at approximately 0.1 within the interval of −0.6 $[eqn]$ to 0 $[eqn]$ (inclusive of 0.6 $[eqn]$ to 0 $[eqn]$ ). For glass/epoxy resin 1, the stress value is stable at about 0.05 within the range of −0.4 $[eqn]$ to 0.41 $[eqn]$ .

Figure 8 illustrates the in-plane shear stress $[eqn]$ in a descending order, with the carbon fiber/epoxy resin demonstrating the highest stress, followed by glass/epoxy resin 1, and then glass/epoxy resin 2. Within the interval of −0.55 $[eqn]$ to 0.55 $[eqn]$ , the maximum shear stress value for the carbon fiber/epoxy resin remains stable at approximately 0.292. For glass/epoxy resin 1, the maximum shear stress value is stable at approximately 0.1605, while for glass/epoxy resin 2, it is stable at approximately 0.1578. The subsequent discussion examines the variation in interlayer shear stress in relation to material parameters.

Figure 9 depicts the interlayer shear stress $[eqn]$ in a descending order, with glass/epoxy resin 1 exhibiting the highest stress, followed by glass/epoxy resin 2, and carbon fiber/epoxy resin. Additionally, the stress value is zero at x = 0, and it exhibits an increasing trend followed by a decreasing trend within the range of 0 to 1.

Figure 10 reveals the interlayer shear stress $[eqn]$ in a descending order, with the carbon fiber/epoxy resin exhibiting the highest stress, followed by glass/epoxy resin 2, and then glass/epoxy resin 1. Within the interval of −0.6 $[eqn]$ to 0.6 $[eqn]$ , the stress value remains constant at zero. Between x = 0.6 $[eqn]$ and 0.8 $[eqn]$ , the stress value exhibits a gradual increase. Between x = −0.6 $[eqn]$ and −0.8 $[eqn]$ , the stress value undergoes a gradual decrease.

Figure 11 shows the order of normal stress $[eqn]$ between layers from large to small: carbon fiber/epoxy resin, glass/epoxy resin 1, and glass/epoxy resin 2. When x is between −0.5 $[eqn]$ and 0.5 $[eqn]$ , the stress value reaches its peak, between 0 and 0.1.

7. Shear Modulus Evolution

To advance the investigation into the degradation of stiffness in composite laminates resulting from crack damage, orthogonal laminates $[eqn]$ composed of glass/epoxy resin 1 and glass/epoxy resin 2 were selected, respectively, with their material properties presented in Table 1. The variable n represents the number of single layers of 0° layers.

Figure 12 illustrates that for both glass/epoxy resin 1 and glass/epoxy resin 2, the shear modulus exhibits a rapid decline as the crack density (CD) increases. Moreover, this decreasing trend is more pronounced with an increasing number of layers in the 0° orientation. These findings suggest a direct correlation between the structural configuration of the laminates and their resistance to crack propagation. In Figure 13, The thickness of the outermost 0° layer is maintained constant at 1 mm, whereas the thickness of the 90° layer varies, yielding distinct values. Specifically, thicknesses ranging from 1 to 5 mm were examined, with an incremental step for the 90° layer thickness increase set at 1 mm. In this analysis, the values for both in-plane and out-of-plane shear moduli are maintained constant, with $[eqn]$ = 4.58 GPa and $[eqn]$ = 3.4 GPa. An inverse relationship is observed between layer thickness and the reduction of $[eqn]$ , with greater thickness correlating with increased reduction. The findings indicate that the thickness of the transverse layer influences the slope of the shear modulus $[eqn]$ curve. Additionally, it is concluded that crack density saturation occurs at an earlier stage with increasing layer thickness. An identical method was employed to quantify the impact of the thickness of the 0° layer. The 90° layer’s thickness was kept constant at 1 mm, while the 0° layer’s thickness was assigned values ranging from 1 to 5 mm, increasing in 1 mm increments. The results are presented in Figure 14. Figure 14 shows the significance of the longitudinal plies on the final shear properties of the cross-ply. This becomes evident as the reduction in shear modulus due to cracking at saturation reaches nearly 16% for 5 mm thick layers of 0, whereas it approaches 42% for 1 mm thick layers (i.e., a single layer of 0).

8. Conclusions

This manuscript is structured into two principal sections. The initial section delves into the stress distribution model of smooth cracked symmetric composite laminates, examining variations in fiber orientation and material properties within the framework of the variational method and the principle of minimum potential energy. Adhering to the principle of minimum complementary energy, the paper derives the governing equations through variable differentiation and establishes the solution constants and undetermined coefficients in accordance with the boundary conditions. This culminates in the determination of stress components across each laminate layer. To validate the accuracy of the developed simulation model, a comparative analysis was conducted using simulation data provided by Ming Kang [7] and Liang Hanbing [16]. In subsequent chapters, the paper employs the complementary energy minimization theory to forecast the degradation of the shear modulus in composite laminates under varying structural configurations. The findings of this research have significantly advanced the understanding of the mechanical properties of composite laminates with matrix cracks. They offer theoretical insights into how different structural parameters, external loads, and material properties influence mechanical properties, thereby contributing to the theoretical foundation and practical application of composite mechanics.

(1) In the instance of smooth open cracks, the alteration in the fiber orientation within θ layer significantly affects the internal positive (tangential) stress distribution. Due to the variation in the fiber orientation angle, the load component in the x (or y) direction corresponding to the fiber orientation also changes. Such variations result in distinct stress distributions across the layers. The specific manifestations are as follows. For in-plane stress $[eqn]$ , the stress value in the cracked layer increases concomitantly with the increase in the fiber orientation angle. When the spatial horizontal coordinate x is between −0.5 $[eqn]$ and 0.5 $[eqn]$ , the in-plane normal stress $[eqn]$ reaches its peak value. However, the stress value in the θ layer exhibits an opposite trend, reaching its minimum value within the same coordinate range. Similarly, for the in-plane normal stress $[eqn]$ , the stress value in the cracked layer decreases as the fiber orientation angle increases. When the spatial horizontal coordinate x is between −0.5 $[eqn]$ and 0.5 $[eqn]$ , the in-plane normal stress $[eqn]$ attains its peak value. Conversely, the θ layer exhibits an opposing trend, with stress values reaching their minimum within the same coordinate range. Concerning in-plane shear stress $[eqn]$ , the stress value in the cracked layer decreases as the fiber orientation angle increases. When the spatial horizontal coordinate x is at −0.3 $[eqn]$ and 0.3 $[eqn]$ , the in-plane shear stress $[eqn]$ reaches its peak value. Conversely, the stress value in θ layer follows an inverse pattern, achieving its lowest value at these coordinates. It can be observed that the bearing capacity of the cracked layer is greatly reduced due to the presence of matrix cracks, with the external load being mainly carried by the θ layer.

(2) For the same opening crack, variations in material parameters significantly influence both the in-plane and interlayer positive (tangential) stress distributions, with the impact varying depending on the stress type and the specific material parameter. For in-plane stress, the in-plane normal stress $[eqn]$ attains its peak when the horizontal coordinate x falls within the range of −0.2 $[eqn]$ to 0.2 $[eqn]$ . For in-plane shear stress $[eqn]$ , the stress reaches its peak when the spatial horizontal coordinate x is between −0.6 $[eqn]$ and 0.6 $[eqn]$ . The normal stress $[eqn]$ in the in-plane does not follow a definitive pattern. However, stress variation patterns differ based on varying fiber orientation angles. For interlaminar stress, when the horizontal coordinate x lies within the range of 0 $[eqn]$ to 1 $[eqn]$ (or 0 $[eqn]$ to −1 $[eqn]$ ), the interlaminar shear stress $[eqn]$ initially peaks before declining to zero. Within the horizontal coordinate range of 0.8 $[eqn]$ to 1 $[eqn]$ (or −0.8 $[eqn]$ to −1 $[eqn]$ ), the interlayer shear stress $[eqn]$ gradually increases, culminating at its peak value. The interlaminar shear stress $[eqn]$ attains its peak within the horizontal coordinate range of −0.6 $[eqn]$ to 0.6 $[eqn]$ . In the horizontal coordinate range from 0.5 $[eqn]$ to 1 $[eqn]$ (or 0.5 $[eqn]$ to −1 $[eqn]$ ), the interlayer shear stress progressively diminishes.

(3) Matrix cracking not only influences the stress distribution within composite laminates but also significantly impacts the evolution of their stiffness. For instance, variations in structural parameters, including the number of single-layer plates and the thickness of the 0° layer, are introduced. Moreover, this decreasing trend is more pronounced with an increasing number of layers in the 0° orientation. These findings suggest a direct correlation between the structural configuration of the laminates and their resistance to crack propagation. With the outermost 0° layer thickness maintained at 1 mm, the impact of transverse layer thickness on the shear modulus is investigated by varying the 90° cracked layer thickness from 1 to 5 mm in 1 mm increments. The findings indicate that the transverse layer thickness significantly influences the slope of the shear modulus curve. An increase in layer thickness is associated with a decrease in the slope of the shear modulus curve. This is primarily due to the increased proportion of the cracked layer within the laminate, which corresponds to a deeper matrix crack and a more rapid natural decline in the shear modulus. An identical methodology was employed to quantify the impact of 0° layer thickness. The 90° cracked layer thickness is maintained at 1 mm, with the 0° layer thickness varying from 1 to 5 mm in 1 mm increments. The slope of the shear modulus curve’s decline for the 0° layer is inverted relative to that of the cracked layer. At saturation, a 5 mm thick layer 0 exhibits a shear modulus reduction of approximately 16% due to cracking, whereas a 1 mm thick layer (i.e., a single layer 0) experiences a reduction of nearly 42%.