FEM Modeling of the Delamination Process in Fabric Composites

The article presents the process of modeling delamination of a woven composite material by the finite element method. The models contain a detailed mesostructure in the form of weave geometry. The delamination occurs in the DCB (Double Cantilever Beam) test with the critical energy release rate as a propagation criterion. The methods used, especially Virtual Crack Closure Techniques (VCCT), have allowed us to present the delamination front changes during the propagation. To further investigate the influence of the mesostructure, additional fully homogenous models were analyzed. Load-displacement graphs, typical for DCB tests, are presented. The obtained results show how the presence of detailed geometry of the composite influences the development of damage.


Introduction
One of the most common modes of composites failure is delamination 1,2 . In laminated composites it occurs when the adhesive bond between the layers fails. Delamination is not a destructive mode of failure, but it leads to significant reduction of laminate's stiffness and allows other modes to arise. Due to its interlaminar occurrence, it is difficult to detect the delamination process and it is even more challenging to determine its further development. To model this type of phenomenon, the finite element method (FEM) is commonly used 1,3 . Complex structure and behavior of composite materials are always limiting factors for all numerical methods, including FEM. To counteract these limitations, simplified geometry and isolation of failure modes are used. This approach might lead to a lower degree of accuracy 1,3 . Damage may occur in different places, but often leads to a similar failure mode.
There are three basic scales in the approach to modeling composite problems: macro, meso and microscale. Recent studies concerning failure modes in woven composites mainly describe multiscale modeling. Llorca et al. 3 describes how this approach allows engineers to model composite elements and account for many failure modes existing in different model scales. In multiscale modeling, damage modes of smaller scales affect other scales like material stiffness degradation. Nonetheless, the explicit effects of the geometry on behavior of the element are not simulated. In other studies concerning woven composites, research limits the analyzed geometry only to the weave itself 4 . Skrzypek and Stadnicki 5 presents a DCB test with a similar weave structure but the details of delamination propagation are not studied. Other studies like Joki et al. 6 do not include weave geometry.
Modeling of the delamination process in the finite element method is done with two methods 1 : Virtual Crack Closure Techniques (VCCT) and Cohesive Zone Model (CZM). To properly use the latter, material properties and failure behavior of the adhesive joint must be known. The VCCT method requires only one critical value -critical energy release rate G IC . In this approach bonds between the layers are modeled by pairing the corresponding nodes from different surfaces. They act as one up to the point when G I reaches the critical value. The separation of the layers is done by deleting the bond between the paired nodes.
In this study, the delamination process occurring in the composite fabric was simulated. A detailed weave pattern of the fabric is modeled near the delamination front. The specimens are subjected to the DCB test in order to compare results to the experiments. Due to the requirements of the CZM technique, the propagation of the delamination front is achieved by means of VCCT procedures. ANSYS engineering software was used to perform the FEM analysis.

Modeling of the DCB Test
The analyzed model simulates the four-layer fabric composite as shown in Figure 1. The entire model was divided into two zones due to high numerical cost of simulating the fabric weave in the entire model. The detailed sector includes the geometry of the plain weave. The second zone simulates composite layers as homogeneous material. To obtain homogeneity, the homogenization process was used. The detailed zone is present only in the middle layers and extends in the proximity of the initial delamination front. As it is widely accepted [1][2][3]5 , in the modeling of the delamination process the critical value of the energy release rate G IC for mode I fracture was chosen as a propagation criterion. Due to the qualitative character of the analysis, no specific adhesive joint between the layers was modeled. The critical value of G IC =166mJ/mm 2 was chosen arbitrarily based on other studies on composite laminates 5,7,8 .
*e-mail: msienkiewicz@meil.pw.edu.pl. Figure 2 depicts the boundary conditions (BC) introduced to the model. They are arranged to simulate the pulling present in DCB tests. Boundary conditions marked 1 and 2 provide dilation of the layers. BC1 pulls the upper edge of the composite in the z axis direction to the desired final displacement δ f (u z = δ f ) and BC2 restricts movement of the lower edge (u x =u y =u z =0) so that it acts like a hinge. Boundary condition 3 prevents lateral (u y =0) movement of the rear areas of the composite.

Weave geometry and material models
The geometry of the fabric can be obtained with simple periodic mathematical functions 9 . This method leads to an idealized weave pattern where all the threads of the reinforcement have the same dimensions. In this study, the plain weave type, presented in Figure 3A, was chosen.
The following equations were used to create weave geometry: The parameters used in the formulas correspond to the dimensions shown in Figure 3A. The following values were adopted: h x =h y = 0.1 mm, a x =a y = 1 mm, g = 0.1 mm and H = 0.3 mm.
The Representative Volume Element (RVE) was created in this way ( Figure 3C). It consists of four threads of the weave connected to a polyester resin. The RVE must represent a repeatable segment of the fabric 4,10 . The dimensions of the analyzed RVE are as follows: 1 × 1 × 0.3 mm (width × height × length). The composite obtained with Equation 1 with the above constant values has the volumetric ratio of reinforcement equal to 41.74%.
The composite fabric consists of two materials: reinforcing threads and resin matrix. Two cases for the thread material were analyzed, both arbitrarily chosen unidirectional composites: UD230 and UD395. In case of the matrix, the polyester resin was chosen for both cases. All the mechanical properties are presented in Table 1 and Table 2 11 .
As mentioned before, the FEM model consists of two zones. The detailed zone involves the thread and the matrix. The homogenous zone needs new material properties that will represent the mechanical characteristic of the weave. To achieve this, the homogenization, well established in the composite material modeling procedure, was used 1,3,7 . This method involves analyzing a small repetitive element of the geometry that can be treated as representative. In our case we performed the homogenization in relation to one RVE. To obtain homogenized mechanical properties, the RVE was analyzed in six numerical tests: three tension tests, one for each of the main axis x, y, z, and three shear tests,  one for each plane: xy, xz, yz. All test configurations are presented in Figure 4. Using the collected data (applied load and obtained displacements) one can calculate all needed material constants for the homogenized material. The results of this procedure are presented in Table 3.

FEM modeling
The geometry of the analyzed case and FEM mesh are presented in Figure 1 and Figure 5, respectively. The model is adjusted to the DCB test taking appropriate boundary conditions into consideration. Thus, the delamination front is present in the center of the model between the middle layers of the composite. As shown in Figure 6 the detailed zone in one layer is 4 by 4 RVEs wide, so there are 32 RVEs present in total. This zone is limited only to the two inner layers. There are two reasons for that: first, the detailed zone requires a very fine mesh, therefore it significantly increases the numerical cost of the analysis and, second, as the damage occurrence is included between the inner layers we limit our interest only to these layers. The mesh in the vicinity of the delamination front must be adjusted to VCCT. The technique requires that two contacting surfaces have an identical mesh in order to create bonding between relevant nodes pairs. Also, due to the method of   numerical calculations it is advised to refine the mesh in the zone of the delamination area. The resultant FEM mesh consists of 235583 elements and 283426 nodes. This size is sufficient for this type of analysis. The detailed zone consists of 210242 elements. This configuration made the model behave more realistically. The obtained distributions of G I , displacements and stresses are smooth therefore proving the mesh to be adequately dense. A crucial aspect in composite modeling is to properly model contact behavior between different materials. In the fabric weave two types of contact occur: thread to thread and thread to matrix. In case of the latter, the issue is simple. The matrix must join the separate materials with the adhesive bond. Therefore, the contact type between the matrix and the thread can be assumed to be bonded, which prevents nodes of contacting surfaces to move independently. Contacting surfaces cannot disconnect. The thread to thread contact, on the other hand, depends on the analyzed case. In some studies 1,3 , a thin layer of matrix is modeled between threads. Another possibility is that the threads are not bonded by any medium and can move separately. In this way another type of damage mode can be introduced to the model -thread debonding. Both cases were analyzed giving us the full spectrum of the real material. The debonded cases (marked further in the text after the type of the element -FRICTIONLESS) model a material containing additional damage. This material has lower strength in bending. On the other hand, a fully bonded material (marked as BONDED) represents the best possible outcome -material without any internal flaws.
Additional analysis of a completely homogenized model was performed for comparative purposes. All assumptions, techniques and boundary conditions were identical to the main case, the only difference being lack of the mesostructured zone.

Results and Discussion
The main analyses were carried out until the delamination front approached the end of the weave structure. As it was stated before, the delamination begins when the value of the energy release rate G I reaches the critical value. Calculation of G I is held on the delamination front and is conducted with the J-integral 1 . This process adds a significant numerical cost to the whole analysis and requires appropriate load application. It was achieved by dividing it into two steps. The first step begins from zero and ends when the delamination starts. This step is divided into several substeps. The next load step begins afterwards and continues to the point when delamination reaches the end of the mesostructure. The division of this  step is much higher; to accurately present the propagation, over 100 substeps were required. The final displacement assigned to BC1 was different for each of the materials: 3.3 mm for UD230 and 4.5 mm for UD395.
Final total displacements of the UD230 model are presented in Figure 7. In case of the UD395 composite, the maximum total displacement is equal to 4.5016 mm.
Similarly to the DCB experiments, the load-displacement data were recorded. They are displayed in Figure 8 together with the results from the fully homogenized model.
The difference between the models can be visualized by showing the distribution of the energy release rate on the delamination front. Figure 9 shows the G I value distribution obtained in both cases of the modeled geometry. This distribution   is obtained in the last moment before crack propagation (point 1 in Figure 8). In all cases the G I value for several nodes surpasses the critical value marked with the dashed black line. It means that in those points the delamination front will move. The distribution differs significantly in each case. The fully homogenized model has a thumb-like shape of the G I value distribution. Additionally, many nodes surpassed the critical value and, therefore, the delamination front will "jump" to the next location. On the other hand, in case of the mesostructure models the distribution is strongly inhomogeneous. In addition, only a few nodes have the value of the energy release rate higher than the critical one. It means that the propagation occurs only in a few places.
The most convenient method to show how the delamination front changes during the separation of layers is to present equivalent stress on the matrix surface. Stress concentrations will occur on the line of the front. Figure 10 depicts those front shapes obtained. Numbers corresponding to the points of loading are indicated in Figure 8. The delamination front is the place of stress concentration.
The presented results are of qualitative character due to the geometrical dimensions of the models and the arbitrarily chosen materials. Nevertheless, extensive impact on the behavior of the specimen caused by the level of accuracy of the model (detailed or homogeneous) can be observed. In case of the fully homogenous model, the delamination front has the thumb shape, taking into account that the weave geometry causes a highly irregular shape of the front, as can be seen in Figure 10. This feature is associated with the distribution of the energy release rate G I across the specimen. It is worth noting that the distribution obtained in the fully homogenized model is similar to the experimental results for the case of the isotropic material of two aluminum plates subjected to the DCB test 12 . Figure 10. Shapes of the delamination front during loading. Numbers correspond to the specific points marked in Figure 8. A -different shapes observed. B -scheme of the initial front shape. C -example images of equivalent stresses observed at the delamination front.
The distribution of the energy release rate on the delamination front is strongly related to the weave geometry ( Figure 11). Peak values can be observed in places where longitudinal threads approach the delamination front. On the other hand, when the transversal thread approaches the front, lower values can be observed. Those differences are related to the directional stiffness of the different threads. During loading, the longitudinal threads are stretched in the direction of their greater stiffness, while for the transversal threads in their less stiff direction. Comparing values from Table 2, it can be noticed that Young's moduli for threads in the longitudinal direction -E xx are significantly greater (14 and 22 times higher for UD230 and UD395 respectively) than E yy . Other peaks in the G I values appear in the areas where larger quantities of the matrix are present in the vicinity of the delamination front. In case of the thread material UD230, the peak values are similar. For the stiffer material UD395 the peak value caused by the proximity of the longitudinal thread is higher than those caused by the resin rich areas.
The force vs displacements graphs (presented in Figure 8) recorded for the fully homogenized models show the behavior typical for VCCT. After the beginning of the propagation, the force decreases in a serrated manner. Including the mesostructure significantly changes the behavior of the element. The delamination starts sooner than in the homogenous models, but it does not lead to an immediate drop in the force value. The graphs even out and tend to approach a steady value. The reason behind that behavior lies in the way of the delamination front propagation. In the fully homogenized model the propagation occurs in "jumps". The energy release rate surpasses the critical value in many nodes. That causes this "jump" when the delamination front moves to the next layer of elements. On the other hand, in the mesostructure models the propagation is continuous. This is a result of the small number of nodes that are surpassing the G IC value every step. A smoother trend can be observed in experimental results for composite materials 8 . Figure 12 presents studies concerning the DCB test and fabric composites. Quantitative comparison to these experimental results is impossible because of the difference in mechanical properties and dimensions of our models. As the delamination front approaches the end of the weave zone, the values in the graphs increase due to the proximity of the homogeneous region.

Conclusions
The presented study showed the importance of modeling the detailed structure of the composite materials. Including the weave pattern significantly altered the behavior of the  models. The presented approach gave an insight into the influence of the geometry of the composite structure on failure development. In fully homogenized models we can observe the characteristic serration region after the initial delamination propagation. In contrast, the detailed zone leads to smoother behavior. This phenomenon is strongly related to the shape and method of the crack's front propagation during loading. Presence of the different materials altered the values of the energy release rate G I calculated on the delamination front. Since this quantity determines the onset and further development of propagation of the crack in VCCT techniques, the difference in behavior of the specimens after the beginning of the delamination can be noticed. The distribution of G I across the delamination front in the homogenized models is much more uniform than in those models with the structure of the composite fabric. The delamination front in case of the homogenized models moves by propagating to the next layer of elements in large groups. This creates a flat and uniform line. In cases when the woven structure was present the process of propagating the delamination front was continuous -some parts of the delamination front were still, while others were propagating. The result of this is a non-uniform delamination front. Further advancements of the presented method should include an increase of the mesostructure size.