SIMULATION OF THE EFFECT LAMINATE SEQUENCE ON DELAMINATION MODE-I WITH ELASTIC COUPLINGS

A numerical study using the ANSYS 19.R3 environment is discussed in this research. This environment depends on the Virtual Crack Closure Technique (VCCT) method to test a double cantilever beam (DCB) according to the ASTM D5528 standard. Four kinds of laminate stacking sequences were considered. According to the results, the distribution of the strain energy release rates obtained along the delamination front in bending-extension and extension-twisting coupling had a good affinity with bending-extension coupling. At the same time, critical fracture toughness values were estimated to be around 87.9% of critical fracture toughness values bending-extension coupling. These results are proof of the bending-extension and extensiontwisting coupling success while testing the proximity to bending-extension coupling results of the DCB beam. These findings are compatible with the standard ASTM D5288. Therefore, the bending-extension and extension-twisting coupling provide a good indication of the delamination resistance during buckling tests of


INTRODUCTION
Composite structures are used in many industrial applications, such as starships, wind turbines, prosthetic limbs for athletes, ships, and cars. This is because of the excellent mechanical properties and performance of such structures. Laminated composite structures often suffer from structural defects during application. Therefore, delamination between laminated layers is one of the most frequent defects in these materials due to their weak interlaminar strengths. A detailed overview of the delamination problems in laminated composites is given by Bolotin [1].
The criterion used to express delamination growth is the strain energy release rate (SERR). It regards the energy dissipated during delamination growth per unit area. The energy required to obtain a crack to identify the onset of delamination within laminates has been researched in numerous studies [2,3].
The double cantilever beam (DCB) test is one of the most common experimental fracture toughness (GI) estimates offered by the American Society for Testing and Materials (ASTM) code ASTM D5528 [4]. This code describes the mechanism used to determine the fracture toughness between layers concerning mode I for composite materials reinforced with specific fibers. The distribution of SERR in mode I is irregular widthwise of the DCB specimen. The GI values can drop at the beam edges, as has been shown by many researchers [5,6,7].
Several researchers [8,9] have noticed different types of defects during the compression (buckling) of laminate at critical states. According to their results, the delamination phenomenon was a prominent defect in the composite materials, focusing on understanding this phenomenon during application problems and load-bearing.
Morais et al [10] were focused on the angles of fiber direction (0/θ, θ/-θ) interface during the DCB delamination tests. According to the results, delamination growth could occur at several fiber angles (e.g., 0/90°, 45°/-45° and 0/45°) [11]. Therefore, it was difficult to determine fracture toughness critical (GIc), and the best solution was 2 DIAGNOSTYKA, Vol. 23. No. 2 (2022) Musefer MA, Ameer ZA, Hamzah AF.: Simulation of the effect laminate sequence on … found by inserting many plies (0/90°) to make the layers more resistant to loads in any direction. Composite structures are expensive, and another important thing is the equipment needed to analyze the results. Thus, the procedures required to solve most problems will be reduced. The development of technology has led to finite element analysis (FEA), which provides a virtual solution at a low cost and little amount of time to delamination problems [12].
With the FEA, a way to obtain the requirements of ASTM D5528 may be offered. For example, delamination growth in composites, in which the crack growth should be slow and stable, may be carefully followed up on the results. FEA is a valuable technique for predicting what happens during the delamination growth of composite materials. Therefore, knowing the fracture toughness of materials may be confirmed during numerical simulation. There are many ways to determine SERR values and delamination growth, including the virtual crack closure technique (VCCT), which is the most widely used method [13,14]. In VCCT, the strain energy released during delamination growth is assumed to be equivalent to the effort necessary to close the fracture back to its original length [15]. The VCCT technique was found to be largely dependent on the size of the finite elements (mesh) in the region where the delamination growth begins and on the size of the load step [16,17].
The beginning was from a suggestion made by Rybicki and Kanninen [14] to apply the VCCT. After that, Shivakumar et al. [18] developed it to include a three-dimensional simulation and their findings were interesting. Therefore, its popularity in numerical analysis has increased. Xie and Biggers Jr. [19,20] used the VCCT method to calculate the fracture toughness and determine the delamination direction front in DCB. De Carvalho et al. [21] used the VCCT method to simulate the development of delamination in cross-ply laminates in the composite structure. Ricco et al. [22] suggested a new approach using the VCCT technique to understand the delamination of laminates subjected to pressure. Turon et al. [23,24] used the VCCT technique to obtain equations to determine the maximum size of elements to ensure the accuracy of the simulation results and their agreement with the experimental results. Therefore, small and organized element sizes and nodes are required by the VCCT technique to obtain precise results despite the increased computational cost of the simulation process [12].
This research article aimed to design a set of composite laminates at different angles. This research was necessary to determine the mechanical behavior of these materials during the DCB test using the advantages and possibilities provided by the virtual crack closure technique (VCCT), depending on the ANSYS 19.R3 environment.

CLASSICAL LAMINATION THEORY
The deformation of a laminate undergoing mechanical loading can be described base on the Classical Lamination Theory (CLT). The moment (M) and force (N) vectors in a laminate composite are defined as follows in CLT [25,26,27]: where: [ ] ̅̅̅ represents the converted stiffness matrix [26]; n represents the current number of plies and k represents the total number of layers; zk represents the length from the laminate mid-plane to the top of the kth ply [28]; K and Ꜫ (or γ) represent curvature and strain, respectively; the upper index ''0" represents the neutral plane. The constitutive equations are written as a simple matrix [29]: (4) is the kth ply center of position dimension in reference to the laminate mid-plane, and tkis the kth ply thickness, with i, j = 1, 2, 3.
The coupling matrix [B], has six different coupled forms: (BS, B0, BF, BL, BT, and BLT ) described in detail by York [29]. As a result, there are twenty-four different types of connected laminates, each with thousands of layups, depending on the total number of plies (n) in the sequence layers.
[B11], [B22], and [B12] affect the coupling between flexural bending and extension; the [B16] and [B26] affect the coupling between flexural twist and extension; and [B66] affects the coupling between flexural twist and shear, [25]: Therefore, the effect of ply sequencing, which is the subject of this paper, should be investigated, which can significantly affect the Mode-I SERR distribution in DCB.
This has been proven by several researchers, for example, de Morais [30] investigated the effect of the delamination front in the measurement of GI might be minimized by an appropriate specimen stacking sequence in DCB, using finite element analyses (FEA).
Thus, a set of class sequences has been selected for consideration which can be expressed using the term layers design and expressed by code. Through these previous studies [31,32,33] the non-symmetry in the strain energy release rate (SERR) over a specimen width in DCB has been quantified in two non-dimensional parameters: The first parameter, Dc is the bending stiffness ratio, which is defined as a function of the [D] matrix terms: The second parameter, Bt, is a measurement of bending-twisting coupling (BT) intensity, which is particularly essential in the DCB test: Where Dij denotes the elements of the specimen's bending stiffness matrices. While determining the stacking sequence, consideration should be given to obtaining with smaller [Dc] ratios, which will reduce the restrictions imposed on [D16, D26]. Therefore, the values required of [D16, D26] three times less than the rest of the elements of the stiffness matrix [D] To eliminate of the effect bending-twisting coupling (BT) on the distribution of the SERR [34], As shown on the values of Bt in the Table 2.
For the DCB test, the SERR (GI) distribution is generally symmetric with minima at the outward borders with respect to the center of the specimen. The "anticlastic curvature" induced by deformation of the specimen arms due to transverse and longitudinal bending coupling causes the difference in G along the delamination front [2]. Davidson [35] suggested that specimens should have a (Dc) of less than or equal to 0.25. As a result, it's necessary that both Dc and Bt values are minimal to provide a uniform delamination onset and SERR (GI) distribution along the crack front.

MODELING FOR DCB BASED ON VCCT
The VCCT approach, which is based on assumptions of "linear elastic fracture mechanics", was the main method for modelling delamination growth [36].
This method is the most extensively used for anticipating delamination growth, which indicates a composite structure's failure. Using threedimensional finite element models, equations are provided to calculate the strain energy release rates. In general, the closer the nodes are to one another (a), the more accurate the equation. The strain energy release rate values will be predicted using equations (9, 10, and 11). The strain energy release rate's Mode-I, Mode-II, and Mode-III components, GI, GII, and GIII, are computed as follows: ∆ is the closed area, ∆ is the length of the elements at the delamination front, and b is the width of the elements when ∆ = ∆ . The forces at the delamination front in column L, row i are represented by ( , , and ). The identical displacements behind the delamination at the top face node row ℓ are denoted as ( ℓ , ℓ , and ℓ ), and at the lower face node row ℓ* are indicated as ( ℓ * , ℓ * , and ℓ * ) as shown in fig 1.
Choosing the modeling of DCB configurations is done by agreeing to the appropriate fraction criterion to describe this behavior during testing and which corresponds to standard double cantilever beam ASTM D5288. Therefore, a suitable fracture criterion including both (Mode-I or Mode-II) and mixed modes should be accepted.
The formula that is used to express a prediction of delamination growth is Reeder Law [37] , which represents the fracture criterion as in equation (12). Geq-c is the equivalent critical, (GIC, GIIC, GIIIC) are critical SERR for Mode-I, Mode-II and Mode-III, respectively and GT = GI + GII + GIII.
Relying on the Reeder Law it is possible to distinguish between different values of mode II and mode III fracture toughness. Therefore, it satisfies the assumption. thus, as this study about to the DCB   Table 3. This convert the Reeder. G eq ≥ G eq−c (14) Note that the criterion was determined in numerical analyses using Finite Element models based on the critical displacement value at which crack onset occurs according to the mathematical formal follows [6] : To apply the above criterion to obtain the critical value of displacement ( critical), the simulation will be performed with the lowest value of displacement ( ), depending on the eq.15.

ANALYSIS
On carbon-epoxy laminates with varied ply sequences in angles, ANSYS (ACP) analysis of Mode I SERR distributions at delamination front as well as delamination length was done. The DCB model had a form of a beam that comprised of two comparable sub-laminates bonded together, and this specimen was divided into two branches at one end by employing an insert that acts as a delamination starter at a length of a˳= 50 mm [41], according to ASTM D 5528 as shown Figure 2. As indicated in Figure 3, the DCB Specimen dimensions were total length L= 150 mm, width b= 20 mm, and total thickness of laminates h= 4 mm, which was determined by the overall number of layers in the laminate (thickness of one layer 0.2 mm). Block branches were pulled in different directions vertically to induce delamination, which was accomplished by applying displacement ( ) to the upper branch.  True boundary conditions (TBC) that comply to standard Mode-I were utilized in this study, as shown in Figure 4: • the loading block lower branch was restrained from i.e, displacement (uix = uiy = uiz =0) and also rotate around the x-axis (θiy = θiz = 0), only the rotation θix=free. • the loading block upper branch was restrained similarly to the lower branch, except the displacement along the y-axis has a variable value (uiy = ). The behavior of DCB beams up to the growth of delamination propagation is described in this work, at different displacement values as to obtain on critical ( ). Table 3. Material properties obtained for carbon-epoxy laminate.
The accuracy of the VCCT calculation depends on the mesh suitable (Quad 4 type) to ensure the greatest accuracy use equal element sizes at initial the crack-tip node (fine) as shown in Figure 5. Table. 4 represented the number and node of elements in ANSYS for DCB. Finite element analysis was conducted on a workstation [Core (TM) i7-9700 CPU-3.00GHz, 32 GB RAM] by using ANSYS WORKBENCH 2019.R3 environment. Fig. 4. DCB beam configuration with fixed support and displacement of applied.

RESULTS AND DISCUSSIONS
The DCB numerical simulations resulted in interesting findings with respect to the couplings chosen for the models. The stiffness extensional matrix, coupling matrix, and bending matrix are shown in Table 2, for which the order of the sequence's laminate was divided between symmetric and antisymmetric relations. There was a balance for both cases, and its purpose was to make the values of the matrices A12 and A26 equal to zero. This was done to eliminate the effect of in-plane shear deformations under in-plane normal loading and to eliminate the effect of normal in-plane deformations under in-plane shear loading. Figure 6 depicts force-displacement curves that were obtained during the DCB pre-crack cycle. The curves are linear before the crack begins to propagate. This occurs at a used displacement value ( ) of about 2 mm, at which point the value that is applied to all the designs is represented by the primitive value of the crack growth. This means that the growth occurs in the insert region; therefore, each design has a certain value for the force that corresponds to the displacement. Figure 7 plots a strain energy release rate as a function of a sample curve width depending on the displacement value ( ) of 2 mm. The behavior of the rate of energy release against the width of the sample shows that the value of the energy is different. The reason for this is that the energy is at a specific place in the insert region depending on the sequence of layers. Therefore, the value may not be regarded as a crucial crack growth value. was applied to address this using the primitive displacement value of 2 mm and the energy release rate value GT for Mode I, which was obtained from a VCCT simulation shown in Figure  7. GT represented a maximum of GI for each design as well as a GIc value for the strain energy release rate in Table 3, which was obtained based on previous studies [13]. After a simulation was conducted based on the critical displacement value, a force displacement diagram was created (Figure 8), which showed the real behavior when the crack growth began for each design. The findings were supported by a delamination front at which a DCB beam had a deformed shape as shown in Figure 9.  Figure 10 shows the critical displacement values for each design. That is, it shows the critical behavior of each design's energy release (GI), which is one of the most important things to consider. It is key because crack growth at the edge of and along the DCB is dependent on this value. The behavior that was found during these simulations was explained by the rate of energy release (the Mode-I GI) for each design. The effect parameters Dc and Bt are shown in Figure 10. After studying, researching, and reviewing many studies, the designs for the layer sequences were selected [42,34,33]. Figure 10 shows the importance of the search results, for which details are listed successively. The GI distributions for all the design laminates are displayed in this figure. GI is symmetrically distributed with respect to the center line.
It was important to follow the ASTM criteria in order to achieve compatibility between the practical and simulated versions of the coupled laminates. The presence and absence of angle θ, which was 45°, clearly affected the layer sequence; this was seen in the front figure at the crack tip. Further, for all the cases, the SERR (GI) was symmetrical with a delamination propagation that began on both sides of the beam. The coupled and laminated bending extension (BE) beam model exhibited flat GI distributions more than the coupled and laminated bending twisting (BT), extension twisting (ET), and bending extension-extension twisting (BE-ET) models. Moreover, the BE model did not exhibit inflections on the beam borders. This result was comparable to a simple laminate and could be explained by the fact that the model produced only an in-plane effect along the z axis. However, the ET coupled design had only an angle of 45° for layer sequences. The laminated beam models exhibited a more complicated GI distribution for the layup type and were narrower (i.e., the ET model exhibited an inflection at the beam borders). The GI values showed upward inflections on both sides of the x axis, and this inflection appeared to be higher for the ET model than for the BT and BE-ET models In an attempt, from a design point of view, angle (45°) was combined with angle (0°,90°) in (BE) to an attempt to reduce the deviation on both sides of DCB in (ET) and this design was successful as it showed the values of GI averaged between two cases, i.e. (BE) and (ET).
From a design perspective, the 45° angle was combined with a 90° angle in the BE model to attempt to reduce the deviation on both sides of the DCB in the ET model. This design was successful as it showed that the GI values averaged between the two models. Further, inserting a percentage of layers also had an effect on the GI distribution and the values of the upward deflections on the two sides. The design of the BT model contained four layers at a 45° angle; thus, the model could be expressed by an angle-ply laminate. The design of the BE-ET model contained six layers also at 45° angles despite the increase in the number of angled layers (45°), but it could be seen that the GI distributions on both sides did not cause upward trends; rather, they approached zero. A reason for this opposite effect was a distribution of the 45° angle between the negative and the positive. That is, it was the distribution between the first layer (45°) and the final layer (−45°). Another reason was a combination of cross-ply and angle-ply that created an equilibrium with respect to matrix stiffness components (as shown in Table 2).
The second coupled and laminated BE design had a layer sequence that appeared to be more closely aligned with ASTM D5528 standards than the other designs; this was because the BE design in question had a 90° angle ply laminate. Also, design.4 (BE-ET) added a set of angles (45°, -45°) to (0°, 90°), It seems closer to agreeing with the criterion ASTM D 5528 and this is what was concluded in this study.
It is necessary to take into account the important parameters that are determined for numerical analysis to avoid undesirable results and to closely detail the critical point, such as the size of the time steps, and types of mesh used to control the crack propagation analysis. During the analysis of the results, was used time steps as shown in Table 5 for each design Which has a clear effect on the accuracy of the results, The ASTM D8825 results are always in good agreement with the peak sawtooth fine values. which is clear in Figure 11.
To obtain values that comply with a standard, it is necessary to focus on choosing the type of mesh in the analysis. So, it was used meshes with uniform crack region element length as shown in Figure 5. These parameters have been discussed in detail by R. Krueger [6].
To explain Figure 11 and 12, the forcedisplacement behavior of each design, it was noticed that in (BE) design it showed the highest stiffness strength by about (83.18 N) due to the fact that the arrangement of layers alternated between (0° and 90°) and this arrangement shows resistance between bending-extension represented in DCB for the ASTM D8835 standard. BE showed higher the length of the delamination is about (83.3 mm), and thus led to an increase in the GI the highest value compared with other designs. The design (ET) showed the least stiffness strength to delamination growth by about (46.67 N) because this design consists of alternating layers between (45°, -45°) which exhibit extension-twisting behavior, which does not meet the ASTM D8825 criterion. Thus, it showed the least slit length growth of about (69.77 mm) with the lowest value GI. Therefore, when combining the design BT and ET in terms of the arrangement of the angles of the layers and get on the designs BT and BE-ET. Therefore, will notice that these two new designs had stiffness strength to them (62.85 and 64.22 N, respectively) and the growth of the delamination length (78.1 and 75.8 mm respectively) and also GI an average between the design BT and the design ET. As shown in Figure 13, a typical curve of GIc versus delamination length, a, of the propagating delamination. The "R-curve" is to be noted with the delamination initiating from the pre-crack at a high value of the interlaminar fracture energy, GIc(init). The value of GIc then rapidly drops as the delamination increase, finally, the value of GIc become rises and forms a plateau. This behavior was observed in the Mode I for all model coupling and appears to be related to the development of a large damage zone at the crack tip as the crack propagates. Where apparent that the design BE showed the highest values GIc than the rest of the designs, while the design ET was the opposite, and the design BE-ET and BT showed an average value of GIc. These results for the values of GIc prove the credibility of all the results that were previously explained in this paper.   Understanding the influence of factor Dc on the distribution (SERR) at the DCB beam is critical. Values obtained from Table 2 that the lowest was for the design BE, which did not show any upward deviation on both sides of the beam for SERR, while the highest value for the design ET, which shows a skew on both sides beam clearly. The new thing in this paper is that the order of layers was obtained, the skew on the side beam of its distribution SERR was very low. This is what was referred to by Davidson [35], the value of Dc must be less or equal than 0.25. Design BE-ET has a value (Dc= 0.22) which is less than (0.25) as shown in Figure 14, Thus that all results are close to design BE. Therefore, parameter Dc is very important for predicting the behavior of a distribution SERR in a DCB beam in the delamination front.

CONCLUSIONS
A numerical simulation study on the applicability of the DCB test for composite laminates with stacking ply sequence specimens was conducted. Elastic couplings with varied fiber orientation angles are shown in the stacking sequences (e.g., bendingtwisting, bending-extension, extension-twisting, and bending-extension with extension-twisting). VCCT was used to simulate the composite laminate design in ANSYS 19.R3 environment. The DCB beam models were loaded and supported following standard specifications.
The major objective of the current study was the estimate the influence of the stacking ply sequences on SERR distributions at the delamination front, critical fracture toughness (GIc), delamination onset, delamination length. Therefore, to try to get on indications for the buckling test which consider is one of the problems in composite structures.
It was concluded, that distribution of SERR is symmetric about the z-axis due to the effect of balance in layer sequences. Also, SERR Musefer MA, Ameer ZA, Hamzah AF.: Simulation of the effect laminate sequence on … distributions obtained along delamination front in bending-extension coupling (BE) and bendingextension with extension-Twisting coupling (BE-ET) were not skewed upward on either side of the beam, while there was a skewed upward on either side of the beam in extension-twisting coupling (ET) and bending-extension coupling (BE). On the other, that the critical fracture toughness (GIc) values of the ET, BT, and BE-ET equal 68.5%, 79.6%, and 89.9% respectively from the critical fracture toughness (GIc) value for BE. The value of the fracture toughness GIc along the DCB beam increases with the increase in the delamination length for each design.
The method used to calculate the critical value of the displacement at which crack formation occurs for each design is one of the most essential aspects of this study. The DCB model's ply stacking sequence had a significant effect on the results. Also, it could be concluded that the VCCT method gives a better predict delamination problems in a composite structure, especially in DCB beams.

Declaration of competing interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.