FEM Modeling of the Delamination Process in Fabric Composites

Abstract

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.

Keywords:
composite damage; RVE; DCB; VCCT; plain weave

1. 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.³ 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⁴. Skrzypek and Stadnicki⁵ presents a DCB test with a similar weave structure but the details of delamination propagation are not studied. Other studies like Joki et al.⁶ do not include weave geometry.

Modeling of the delamination process in the finite element method is done with two methods¹: 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.

2. 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-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² was chosen arbitrarily based on other studies on composite laminates^5,7,8.

Figure 1
Geometry of the DCB model. Geometry which represents weave fabric is present in the middle of the element. Numbers in squares depict boundary conditions applied to the model. The zoomed section represents one of the detailed zones. One of the 16 Representative Volume Elements (RVE) that build the upper weave pattern is highlighted. Units in mm.

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.

Figure 2
Boundary conditions applied to the finite element model. BC1 and BC2 (marked with corresponding numbers in squares) provide proper dilation of the composite layers. BC3 ensures that all of the nodes of the back area move dependently.

2.1. Weave geometry and material models

The geometry of the fabric can be obtained with simple periodic mathematical functions⁹. 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:

Figure 3
Geometry and FEM models of RVE. A – dimensions of the plain weave. B – weave pattern present in the RVE. Different threads orientations are marked with arrows. C – full FEM model of RVE consisting of four threads and a matrix.

$\begin{array}{c}{z}_y=-\frac{h_y}{2}cos\frac{\pi x}{a_x},\\ z_x=\frac{h_x}{2}cos\frac{\pi y}{a_y}.\end{array}$ (1)

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 ¹¹.

Table 1
Matrix material mechanical properties.

	E [MPa]	ν	G [MPa]
Polyester resin	3000	0.316	1139.8

Table 2
Mechanical properties of the thread materials.

	Exx[MPa]	Eyy[MPa]	Ezz[MPa]	νxy	νyz	νxz	Gxy[MPa]	Gyz[MPa]	Gxz[MPa]
UD230	121000	8600	8600	0.27	0.4	0.27	4700	3100	4700
UD395	209000	9450	9450	0.27	0.4	0.27	5500	3900	5500

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.

Figure 4
Six cases used in the homogenization process carried out on the RVE. Cases A1-3 represent the tension tests in three axes: x, y and z. Cases B1-3 represent the shear tests in three planes: xy, yz and xz. Shear test models include additional non-deformable elements to ensure pure strain state.

Table 3
Mechanical properties of the homogenized material.

UD230	Exx[MPa]	27321.81	Gxy[MPa]	2343.62	νxy	0.047	νxz	0.458
	Eyy[MPa]	27325.09	Gyz[MPa]	1442.69	νyz	0.459	νyx	0.047
	Ezz[MPa]	5217.75	Gxz[MPa]	1439.51	νzx	0.086	νzy	0.085
UD395	Exx[MPa]	44691.54	Gxy[MPa]	2654.84	νxy	0.038	νxz	0.477
	Eyy[MPa]	44451.56	Gyz[MPa]	1878.50	νyz	0.476	νyx	0.037
	Ezz[MPa]	2190.10	Gxz[MPa]	1875.01	νzx	0.061	νzy	0.060

2.2. 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.

Figure 5
Finite element model of the DCB test specimen.

Figure 6
FEM model of one of the weave patterns. Longitudinal threads are blue, and transversal threads are red.

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.

3. 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¹. 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.

Figure 7
Final deformation of the UD230 models. Displacement δ applied to the model is represented with an arrow. Units in mm. A – mesostructured frictionless model. B – homogenized model.

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.

Figure 8
Force vs. displacement in the DCB numerical test. δ stands for the displacement of the edge with the BC1 and F denotes the measured reaction force (Point 1 refers to the moment prior to the beginning of the propagation).

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.

Figure 9
Values of the energy release rate on the delamination front for point 1 of Figure 8. The dashed black line represents the critical value 166 mJ/mm².

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.

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 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¹².

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.

Figure 11
Distribution of the energy release rate on the delamination front prior to the propagation start. A – the UD230 model. B – the UD395 model. The yellow dashed line represents the delamination front.

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⁸. 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.

Figure 12
Force – displacement graphs obtained from DCB experiments⁸.

4. 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.

5. References

• ¹
  Camanho PP. Numerical modelling of failure in advanced composite materials. San Diego: Elsevier Science; 2015.
• ²
  Daneshjoo Z, Shokrieh M, Fakoor M, Alderliesten R. A new mixed mode I/II failure criterion for laminated composites considering fracture process zone. Theor Appl Fract Mech. 2018;98:48-58. http://dx.doi.org/10.1016/j.tafmec.2018.09.004
  » http://dx.doi.org/10.1016/j.tafmec.2018.09.004
• ³
  Llorca J, González C, Molina-Aldareguía JM, Segurado J, Seltzer R, Sket F, et al. Multiscale modeling of composite materials: a roadmap towards virtual testing. Adv Mater. 2011;23(44):5130-47. http://dx.doi.org/10.1002/adma.201101683
  » http://dx.doi.org/10.1002/adma.201101683
• ⁴
  Riva E, Nicoletto G. Modeling and prediction of the mechanical properties of woven laminates by the finite element method. In: Guagliano M, Aliabadi MH, editors. Fracture and damage of composites. Boston: WIT Press; 2005. p. 105-125. (WIT Transactions on State of the Art in Science and Engineering). http://dx.doi.org/10.2495/978-1-85312-669-7/05
  » http://dx.doi.org/10.2495/978-1-85312-669-7/05
• ⁵
  Skrzypek J, Stadnicki J. Virtualization of interlayer crack propagation in CRF composite. Mechanik. 2015;(2):146/35. http://dx.doi.org/10.17814/mechanik.2015.2.57
  » http://dx.doi.org/10.17814/mechanik.2015.2.57
• ⁶
  Joki R, Grytten F, Hayman B, Sørensen B. Determination of a cohesive law for delamination modelling: accounting for variation in crack opening and stress state across the test specimen width. Compos Sci Technol. 2016;128:49-57. http://dx.doi.org/10.1016/j.compscitech.2016.01.026
  » http://dx.doi.org/10.1016/j.compscitech.2016.01.026
• ⁷
  Matthews FL, Rawlings RD. Composite materials: engineering and science. Boca Raton: CRC Press; 1994.
• ⁸
  Sarrado C, Turon A, Costa J, Renart J. On the validity of linear elastic fracture mechanics methods to measure the fracture toughness of adhesive joints. Int J Solids Struct. 2016;81:110-6. http://dx.doi.org/10.1016/j.ijsolstr.2015.11.016
  » http://dx.doi.org/10.1016/j.ijsolstr.2015.11.016
• ⁹
  Kuhn JL, Charalambides PG. Modeling of plain weave fabric composite geometry. J Compos Mater. 1999;33(3):188-220. http://dx.doi.org/10.1177/002199839903300301
  » http://dx.doi.org/10.1177/002199839903300301
• ¹⁰
  Grufman C, Ellyin F. Numerical modelling of damage susceptibility of an inhomogeneous representative material volume element of polymer composites. Compos Sci Technol. 2008;68(3-4):650-7. http://dx.doi.org/10.1016/j.compscitech.2007.09.018
  » http://dx.doi.org/10.1016/j.compscitech.2007.09.018
• ¹¹
  Ansys. Ansys® Academic Research Mechanical, Release 19.2. Canonsburg: Ansys.
• ¹²
  Sun C, Zheng S. Delamination characteristics of double-cantilever beam and end-notched flexure composite specimens. Compos Sci Technol. 1996;56(4):451-9. http://dx.doi.org/10.1016/0266-3538(96)00001-2
  » http://dx.doi.org/10.1016/0266-3538(96)00001-2

Publication Dates

• Publication in this collection
  26 June 2020
• Date of issue
  2020

History

• Received
  13 Dec 2019
• Reviewed
  22 Apr 2020
• Accepted
  24 May 2020