A variational approach for predicting initiation of matrix cracking and induced delamination in symmetric composite laminates under in-plane loading

2.1 Perturbation stress field analysis

Consider an angle-ply [ϕm(2)/ψn(1)]s composite laminate that is cracked under tensional static loading. All of the cracks have equivalent distances between each other (Figure 1).

Figure 1:

Geometry of a symmetric laminate containing cracks in the middle layers in a global coordinate system.

A list of symbols and parameters is compiled in Table 1, as defined in the text and the following formulations.

The maximum number of layers in this analysis is four sublaminates; however, researchers are able to investigate for arbitrary lay-ups by handling their calculations following Ref. [11]. In addition, the thicknesses of these laminas can be considered arbitrary. As shown in many papers [11], [12], [24], these calculations are the same in estimating basic science; however, generally, their boundary conditions are different. The first step of calculations is assessing the stiffness degradation. For this purpose, it is needed to transfer the unit cell’s coordinate to the new one. In fact, if it is assumed to continue calculations based on the project just like Hashin did in 1985 [10], changes should be made in its primary strategy. In order to have 90° plies in the middle sublaminates, the rotation around Z axis should be [−ψn(1)+90]. Indeed, this method carries out the data to the previous calculations with cross-ply sublaminates. New lay-up for this composite can be displayed in the form of [S/90]s. The term S denotes [ϕm(2)−ψn(1)+90] in the general lay-up (Figure 2).

Figure 2:

Cracked laminate in the rotated coordinate.

As illustrated, under in-plane loading, which can be arbitrary, these cracks are constructed. After the rotation, it is expected that the type of forces and their amounts would be altered. Lowercase letters represent new coordinate characteristics. Next, it is shown that in some rotations depending on their angles and their quantities, occasionally they would be reversed into compressive loading. By these interpretations, we bring forward external loads for the primary coordinate as F_XX, F_XY, F_YY and for the secondary coordinate in F_xx, F_xy, F_yy. It is definitely seen that there are always some papers that have used some different types of loadings and lay-ups [34] (which work on stiffness reduction of off-axis laminates; neither induced delamination nor matrix cracking analysis for angle-ply laminates); however, in this research, it is tried to assess general loadings and stacking sequence.

Now, it is the time to explicate the composition of stress field in a cracked laminate. In a fragment, it can be expressed as an admissible stress field by a coincidence of the stresses in the uncracked material and some special unknown perturbation stresses that are due to the presence of cracks.

where i=1, 2 specifies the number of laminas, and u and p denote uncracked and perturbation, respectively. Proportionate with what is shown in Figure 1, h=t₁+t₂ and a dimensionless ratio λ=t2t1 are defined. By virtue of symmetry, it is tried to analyze only half of the laminate. In this model, as is observable in Figure 1, there are several boundary and continuity conditions that can help arrive at a solution for the perturbation stresses. The following conditions have been used in calculations:

1. There are no stresses on the out-plane surfaces (at z=∓h, σ˜mz(2)=0).

2. Using continuity conditions between interfaces, it can be seen that at z=±t₁, σ˜mz(1)=σ˜mz(2).

3. Free tractions on the crack surfaces denote that (x=±a, σ˜mx(1)=0).

4. The behavior of an uncracked laminate with a cracked laminate just like the work done by Vinogradov and Hashin [24] can be simulated; thus, the effects of local damages can be ignored.

Using the equilibrium equations and also balancing the membrane forces, some perturbation stresses that have impacts upon the stiffness degradation are achieved. Indeed, these amounts directly influence stiffness reduction. According to the steps that have been derived by Vinogradov and Hashin [24], the following can be written.

For lamina 1 (the middle lamina):

For lamina 2 (the outer ply):

where all of these stresses are the perturbation stresses in two sublaminates. The applied boundary conditions are denoted by

Note that these implicit results of boundary conditions will be used in calculating the final answers of ordinary differential equations of (θ, ϑ, η), which are derived in the next step.

2.2 Variational formulation

Now, the principle of minimum complementary energy is going to be used for obtaining the perturbation stresses. For this purpose, the strain energy release rate must be calculated for transverse cracked laminates. It is very important to know about how the strain energy has to be derived in a model. For a loaded unit cell, the complementary energy functional can be defined as

As it was mentioned in some papers like variational principles for generalized plane strain problems and their applications [32], the correct shape of the energy criterion should not be used in the form of strain energy; however, given that there is a linear elastic assumption, this equivalent is acceptable. Therefore, all equations that are considered have to be under this assumption. In the presence of cracks, for all elastic materials, Eq. (5) can be expressed in the following form:

where U^u and U′C are the strain energies for an uncracked body and due to perturbation stresses, respectively. As mentioned in Hashin’s papers [10], [24], [35], the value of the complementary energy using admissible stresses can be derived from the principle of minimum complementary energy or, in another words, it can be written as

where the terms C^∗ and C^u are the effective compliance matrix of cracked and uncracked bodies, respectively; V is the volume of the cracked laminate; σ¯ specifies the average in-plane stresses equal to applied stresses; and σ^p denotes the perturbation stresses that are constructed by the presence of cracks. Actually, it is needed to minimize the right-hand side of Eq. (7) in order to derive the closest amount of the strain energy to a practical value. For this purpose, it has been suggested to minimize an integration that is produced by multiplying the perturbation energy density by other quantities, which are described in equations.

Until here, it has been attempted to put forward a model that is restricted by two cracks in the volume of a fragment with 2a length and |z|≤h. Due to the factors in Figure 1 and conditions that have been explained before, for half of a unit cell, it can be written that

where W⁽ⁱ⁾ and εmnp(i) are the strain energy density and strain value of the ply (i). It is worth to mention that the strain energy density is a scalar function. The sum of integration of all sublaminates over the given volume provides the total strain energy of the body. To calculate the strain amount, or in another words (εmnp(i)=C¯(i)σmnp(i)), we have to compute the compliance of the matrix unidirectional composite. For this purpose, as it was suggested by Vinogradov and Hashin [24], it is better to rotate the coordinate as was mentioned in Section 2.1. Embedding the expressions of perturbation stresses [Eqs. (2) and (3)] and stress energy densities in the perturbation strain energy, an integral calculus with constant coefficients is derived that assists researchers in computing the value of this expression to totalize with the strain energy of the uncracked body. The general terms of this integration is expressed here:

For solving this integration from the Euler-Lagrange method, it must be converted to dimensionless form with these changes (ζ=xt1, ρ=at1). The above equation can be illustrated as

where

and the multipliers are the known amounts that are detailed in the Appendix. To compute the strain energy due to perturbation stresses, the unknown functions [θ(ζ), ϑ(ζ), η(ζ)] must be calculated from minimizing this integration. Because of three unknown functions, three Euler-Lagrange equations are needed, which are obtained from the below relations:

Substituting the result of Eq. (14) into Eqs. (12) and (13), a new one is derived that is only a function of two variables.

which helps us rewrite Eq. (11) as

As shown in the above equation, one of the unknown functions is omitted easily (η). The other independent unknown functions (θ, ϑ) will be determined in the coupled linear differential equations when Euler-Lagrange equations are implemented for two times. Thus, it can be written as

where Tia are the constant coefficients that can be obtained according to the formulations based on the previous coefficients before implementing the variational approach. Depending on the boundary conditions, this coupled linear differential equation can be converted to the model that is suggested by Refs. [10], [23], [24].

2.3 Compliance of cracked unit cell

In order to solve coupled linear differential equations, there are several solutions that are suggested by commercial mathematical software. The authors prefer to utilize a more rational way, which is an analytical solution. Indeed, from the standpoint of mathematics, it is not logical that numerical methods are used instead of analytical ways when both solutions are available.

For these grounds, let us consider a matrix that is deduced from the coupled linear differential Eqs. (17) and (18):

As it is clear to have a specified solution for Eqs. (17) and (18), the determinate of the above relation [Eq. (19)] should be zero, and thus it can be written that

where

Due to the simplification of the equations, it is assumed that k=r². Depending on the boundary conditions, it can be observed that there might be two complex roots and one real positive root.

This is not the right place to describe more points about the process of solution. With these interpretations, it is better to refer to differential equation books like that of Simmons [36]. It has indicators to solving these kinds of problems.

The general model of a system with these roots can be noted in the form of

It is essential to note that this problem has symmetry of about (x=0). Here, multipliers of θ(ζ) and ϑ(ζ) have relations with each other. Given Eqs. (21) and (22) in calculations lead authors to express that

The insertion of the expression in Eq. (4) helps us obtain the particular solutions of Eqs. (24) and (25). At first, it is needed to compute the constants C_(i) and then there are precise answers for unknown functions (θ, ϑ). This is a uniform solution with a less detailed description for calculating the stiffness degradation in comparison with Vinogradov and Hashin’s work [24]. Briefly using these functions, an equation can be obtained for compliance of the cracked body, which helps authors to evaluate a true value for strain energy in a cracked body. The expression of compliance for a cracked body can be obtained as

According to the obtained equation, it is noted that in the variational standpoint, the minimum amount of the perturbation strain energy in the unit cell has the principal impact on the effective compliance matrix (C^∗). Now, in this model, which is able to estimate the stiffness of a cracked body, it can be said that these calculations would not depend on the angles of sublaminates and it can be arbitrary.

5.1 Cross-ply lay-ups

Initially, it is better to introduce the same materials that had been used in Nairn and Hu’s results [31].

EA=128 GPa, GA=4 GPa, νA=0.3.ET=7.2 GPa, GT=2.4 GPA, νT=0.5, t1=t2=0.14 mm.

According to the utilized conditions by Ref. [33], the applied stresses are σxx00=100 MPa,  σxy00=0, σyy00=0. If the presented Eqs. (35) and (42) are thoroughly investigated, the obtained results have incredible agreements to Ref. [33] for a (0₂/90₄)s carbon fiber/epoxy as shown in Figure 5, even though this process is much easier than theirs because of the ignorance of simplified additional equations.

Figure 5:

Comparison of the obtained energy release rate with Li’s results for [0₂/90₄]s carbon/epoxy laminate under σ₁=100 MPa [33].

Figure 5 illustrates the variation of energy release rates due to two forms of most possible damages with the increase of crack densities. It is shown that between 0 and 0.51 crack densities, the odds-on type of damage is formation of matrix cracking; however, from this critical point to the lowest distance between cracks, delamination must be observed. Generally, this critical point helps designers to know something more about the next possible damage. As shown in Figure 5, the size of the unit cell plays an essential role in the determination of the kind of damages. Moreover, this analogy proves the point that the new method of achieving energy release rate in the variational approach has the same precision that theirs have.

In the next attempts, the authors consider diverse cross-ply lay-ups of carbon fiber/epoxy, which are distinguished by their thicknesses. Likewise, the critical points are obtained easily, so that it is straightforward to estimate the effects of the number of cracked sublaminates in dislocation of critical points relative to crack densities. For this purpose, various types of fiber/epoxy cross-ply laminates are considered to report the influence of adding cracked laminas in Figures 6–8.

Figure 6:

Crack density vs. energy release rate for [0/90]s carbon/epoxy laminate under σ₁=100 MPa.

Figure 7:

Crack density vs. energy release rate for [0/90₂]s carbon/epoxy laminate under σ₁=100 MPa.

Figure 8:

Crack density vs. energy release rate for [0/90₃]s carbon/epoxy laminate under σ₁=100 MPa.

In all of these calculations, the considered stresses are  σxx00=100 MPa,  σxy00=0, σyy00=0; however, the authors try to draw some graphs that have non-zero σxy00 and σyy00. Nevertheless, Figures 9–11 are provided as samples with other remote loadings. In these paradigms, the stresses are  σxx00=100 MPa, σxy00=1 MPa, σyy00=10 MPa. As it is detectable in these cross-ply laminates, the increasing thicknesses of cracked sublaminates lead to regressive conditions, or, in another words, induced delamination is more prospective by declining the critical point relative to crack density.

Figure 9:

Crack density vs. energy release rate for [0/90]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 10:

Crack density vs. energy release rate for [0/90₂]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 11:

Crack density vs. energy release rate for [0/90₃]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

5.2 Angle-ply lay-ups

Diverse angle-ply laminates are chosen to compute their energy release rates when their stresses in local coordinates are the same as the precise amount that cross-ply laminates have. Indeed, the authors try to examine how differentiation of angles affects the energy release rate for induced delaminations and the next matrix cracking in a unit cell. According to these assessments, it is decided to provide Table 2, which introduces some essential factors about the stress and its angles for local coordinate.

Relatively straightforward calculations specify the point that in order to achieve angle-ply energy release rates, the expressions have to be investigated according to their local coordinates as Hosseini-Toudeshky et al. [8] did in their micromechanical approach. The results are drawn in diverse line charts, which are shown here.

5.2.1 [45/−45n]s Angle-ply laminates

Figures 12–14 show changes in the energy release rates of matrix cracking and induced delamination for [45/−45]s, [40/−45₂]s and [40/−45₃]s. There is a modest increase in the energy release rate for matrix cracking and then it plunges sharply, although this steep decrease is not the same as the decline that is seen in the induced delamination. The critical points for Figures 12–14 are 0.81, 0.605 and 0.47, respectively. It should be underlined that the decreasing trend of critical points tends to be relative to the increasing number of thicknesses in the cracked laminas. In addition, according to the data shown, the amounts of energy release rate by increasing the number of cracked laminas are escalated.

Figure 12:

Crack density vs. energy release rate for [45/−45]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 13:

Crack density vs. energy release rate for [40/−45₂]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 14:

Crack density vs. energy release rate for [45/−45₃]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

5.2.2 [40/−40n]s Angle-ply laminates

Figures 15–17 depict how inclinations for the energy release rates have changed for matrix cracking and induced delamination as well as Figures 12–14 illustrated. The most obvious changes, which are notable in comparison with the previous figures, are its measures for energy release rates. Indeed, whatever the stacking sequences of local coordinates ([S/90]s) is closer to [180/90]s, the energy release rates see a slower change and vice versa. It can be clearly seen that the energy release rate for the next matrix cracking in [40/−40]s starts with 20.77 J/m²; in contrast, this amount for [45/−45]s begins with 3.81 J/m².

Figure 15:

Crack density vs. energy release rate for [40/−40]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 16:

Crack density vs. energy release rate for [40/−40₂]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 17:

Crack density vs. energy release rate for [40/−40₃]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Another distinguishing feature is their critical points. The variations of energy release rate versus crack density for matrix cracking and delamination intersect at 1.74 (1/mm), 0.725 (1/mm) and 0.413 (1/mm) in Figures 16–18. The trend of decreasing energy release rates is virtually the same as the other stacking sequences. As can be observed, some of these graphs are totally outnumbered by increasing the number of sublaminates.

Figure 18:

Crack density vs. energy release rate for [30/−30]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

5.2.3 [30/−30_n]s Angle-ply laminates

In this section, the widest difference can be seen in the amounts of energy release rates that Figures 18–20 demonstrate with the previous figures. As illustrated, the critical points are 1.96 (1/mm), 0.925 (1/mm) and 0.524 (1/mm), respectively, for Figures 18–20.

Figure 19:

Crack density vs. energy release rate for [30/−30₂]s carbon/epoxy laminate under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

Figure 20:

Crack density vs. energy release rate for [30/−30₃]s carbon/epoxy laminate for under σ₁=100 MPa, σ₁₂=1 MPa and σ₂=10 MPa.

This opportunity has been granted to express the standpoint that local stresses and also angles, when the middle layers are 90°, are determinative. In other words, for initiation of delamination and intralaminar matrix cracking, the key point is that how the global stresses would change and how these amounts affect the energy release rates according to their crack densities.

As the Euler-Lagrange equations are derived, it can be said that variational approaches can be acceptable if the orientation of fibers would be restricted in two layers and the intralaminar cracks would occur in the middle layer.

In Figure 21, the changes in critical densities for carbon/epoxy [0/90n]s, [45/−45n]s, [40/−40n]s, [30/−30n]s and [15/−15n]s laminates are illustrated with respect to the number of middle cracked plies. The results demonstrate a rational trend according to the altering angles. In other words, whatever the number of cracked laminates are increased, the risk of delamination increases in proportion to the particular lay-up that is being investigated. According to previous studies [8], it can be proved, when the coordinates and the amounts of stresses are equal in local calculations, that critical crack density points have broadly similar amounts.

Figure 21:

Number of middle cracked plies vs. critical crack densities for distinct carbon/epoxy laminates.

Table 3 shows the critical energy release rates, critical crack densities, and energy release rates due to matrix cracking and induced delamination for a minute unit cell (c=0.01). G_c denotes the critical energy release rate, which represents the precise energy release rate for exceeding induced delamination in relation to matrix cracking as possible damages. Table 3 depicts that whenever stacking sequences in local coordinates get far from [0/90n]s, the energy release rates increase dramatically. Furthermore, the amounts of critical crack densities are growing rapidly. In this table, it can be seen that almost all critical crack density values are decreasing by increasing the damaged area. In other words, decreasing from 0.81 to 0.47 in [0/90n]s shows that the possibility for matrix cracking is shrank by tripling the cracked sublaminates. This trend seems set to continue for other stacking sequences.

Table 3:

The critical energy release rates, critical crack densities and energy release rates due to induced delamination and matrix cracking for a unit cell with c=0.01 (1/mm) under σ₁=100 MPa, σ₁₂=1 MPa, σ₂=10 MPa.

No.	Lay-up	Critical crack density (1/mm)	Gc (J/m2)	Gm/c=0.01 (1/mm)	Gd/c=0.01 (1/mm)
1	[0/90]s	0.81	3.79	3.813	3.785
2	[0/902]s	0.61	11.06	13.636	11.032
3	[0/903]s	0.47	23.48	32.028	23.477
4	[45/−45]s	0.81	3.79	3.813	3.785
5	[40/−452]s	0.61	11.06	13.636	11.032
6	[45/−453]s	0.47	23.48	32.028	23.477
7	[40/−40]s	1.73	13.088	20.77	13.59
8	[40/−402]s	0.73	41.19	49.382	41.133
9	[40/−403]s	0.411	84.19	90.66	84.076
10	[30/−30]s	1.975	67.222	140.508	74.702
11	[30/−302]s	0.923	222.014	305.058	224.559
12	[30/−303]s	0.529	451.260	512.55	450.651