Detailed Finite Element Models for the Simulation of Laser Shock Wave Response of 3D Woven Composites

1. Introduction

Novel aircraft engine fan blades utilize a multi-material structure composed of a 3D woven CFRP core partially covered by an adhesively bonded thin metallic layer. The end-of-life management (reuse or recycle) of such a configuration involves many challenges, with one of them being the disassembly of the components [1,2,3,4]. Numerous methodologies have been explored to address the debonding challenge, encompassing various physical and chemical phenomena [5]. However, many of these methods have been proved less suitable in scenarios involving composite materials, especially when the preservation of the integrity of the materials is essential. Recent developments on the laser-shock technique have revealed the potential of the method to be used for disassembly of adhesively bonded CFRP and metallic parts [6]. One advantage of this method is it’s the precise load application at the bondline by calibrating the laser parameters [6].

The laser-shock technique involves the interaction between a high-power beam laser and a target material. This interaction initiates a plasma expansion process, induced by the recoil momentum of the ablated material, resulting in a pulsed pressure. Notably, when plasma expansion occurs within a confinement regime, where the dielectric material is transparent to the laser, both the pressure level and duration are increased. Thus, an elastic precursor shock followed by an elastic-plastic compression shock is generated and propagates in the material. After the plasma expansion, the surface unloads, and a plastic-decompression shock alongside an elastic-plastic decompression shock starts to propagate, creating a release wave. When the release wave interacts with the elastic precursor shock wave, it creates high localized tensile stresses. Demonstrating versatility, this technique has already been successfully used to induce controlled delamination in composite materials [7], paint stripping on metallic and composite substrates [6,8,9].

In [6], the authors have demonstrated experimentally for the first time the application of the laser shock technique for the disassembly of adhesively bonded CFRP/metallic coupons. To get a deeper insight of the physical mechanisms and to understand the role of the parameters that must be calibrated in order for the achieve debonding without damaging the CFRP part, a numerical simulation model of the process has been developed in [10]. Although that model has served its role efficiently a careful look at the numerical and experimental findings brings the need for a more detailed understanding of the shock wave propagation in the 3D woven CFRP material and for a validation of the predicted induced damage state.

In recent years, there have been reported several simulation attempts of the laser-shock process [8,9,10,11,12,13], with most of them focusing on the use of homogenized models to investigate or optimize the stress field in specific areas. While for metallic materials, there have been a few works studying in detail the laser shock wave propagation, such as the work of [14], for composites, and especially for 3D woven composites, there is a notable gap.

On the other hand, the extension of use of 3D woven composites in engineering applications has raised the need for developing reliable simulation models. Two useful reviews of such models can be found in [15]. While for simple geometries and loading conditions homogenized models [16] are sufficiently efficient, and at the same time very fast and easy to be developed, for complicated geometries and dynamic loading conditions [17,18], detailed models are needed. If the mechanical behavior of structural parts or structures is to be simulated, multi-scale models are used. Such models usually employ homogenization of a representative volume element which represents the exact geometry of the 3D woven composite [19,20,21,22,23,24,25,26]. Additionally, Meso-scale modeling has been applied to study structural components of bi-axial weft-knitted composites in [27]. Only one paper has been published on the modeling of shock loading of 3D woven composites. In that paper, Pankow et al. [28] have developed FE models to simulate the deformation response and damage development of 3D glass-reinforced composites subjected to shock loading. Two models were presented; the first comprised of a layered homogenous model and in the second model each layer is represented to capture the interactions between fiber tows and matrix, considering the discrete non-homogeneous material distribution in each layer.

In this paper, we have developed a detailed model to study the mechanical response of a 3D woven composite subjected to laser shock loading and to assess the efficiency of a homogenized model used for the same purpose.

2. Materials and Methods

2.1. Detailed FE Model

2.1.1. Design of the 3D Woven Geometry

The 3D woven geometry was designed in TexGen [29] using the Python code. The specimen's woven pattern comprises two noticeably distinct sections: the bulk and the surface. Combining these elements into a unified woven model necessitates the utilization of the scripting features of the TexGen software [29]. Initially, a layer-to-layer 3D weave is established to develop the bulk structure of the specimen. This procedure is direct and extensively detailed in the scripting manual, available on the TexGen website. To incorporate the top and bottom layer in the designed model, each yarn is designed individually. The implementation subroutine has been created parametrically in relation to the yarn dimensions and the inputs required are the positions of the offset for each layer, as it is described by Figure 1, and the identification numbers of the yarn. The geometrical parameters such as the warp dimensions and the gap distances are the global parameters that are used for the whole script. The nodes are progressively denser, while approaching the curvature of the offset and their position underneath the wrap fiber bundle, for the top binders, or on top of the wrap bundle, for the bottom binders follows the super-ellipse shape of the wrap bundle’s cross section. The node’s coordinates for the super-ellipse, $C_x$ and $C_y$, are calculated by [30]

$$C_x\left(t\right)=\frac{w}{2}\mathit{cos}2\pi t,0\le t\le 1$$

(1)

$$C_y\left(t\right)=\left\{\begin{array}{c}\frac{h}{2}{\left(\mathit{sin}\left(2\pi t\right)\right)}^n, 0\le t \le 0.5\\ -\frac{h}{2}{\left(\mathit{sin}\left(2\pi t\right)\right)}^n, 0.5\le t \le 1\end{array}\right.$$

(2)

where, h is the height of the super-ellipse, n is the order of the ellipse and t is a local variable ranging from 0-1 indicating the start and end of the super-ellipse. Combining the two geometries using the LayeredTextile function, the final geometry of the 3D woven composite, shown in Figure 2, was created.

2.1.2. Meshing of the Detailed Model

The model's discretization employs voxel elements which have been proven to be accurate even in the presence of high curvature in certain binder yarns [31]. This choice was based on two primary considerations. Firstly, the utilization of voxel elements permits the application of the advanced material models which are detailed in Section 4, and secondly, the voxel elements offer a more uniform distribution of stresses through the material's thickness compared to tetrahedral elements. The meshing process was conducted within the TexGen software and subsequently exported as an Abaqus file, which was then imported into the LS-Dyna. Specifically, the voxel elements employed in the model utilize ELFORM 1, which is a constant stress solid element.

The aforementioned process exhibits inherent limitations that warrant careful consideration. Specifically, the LS-Prepost software encounters challenges in converting the exported Abaqus input file to an LS-Dyna key file for denser meshes. Considering the scale of the geometry that was described in Section 2, it is not difficult to exceed the software capabilities. The conversion process’s upper limit is not described in any of the manuals, but from rigorous testing utilizing different high-performance computer systems it was concluded that the software was not designed to manage large file conversions. To address the issue, a reduced geometry was employed to conduct a mesh density study. The dimensions of the reduced model are 10 mm by 10 mm, which encloses the center of the original model, while preserving the full thickness.

2.1.3. Mesh Density Study for the Detailed Model

The selection of the mesh density in the detailed FE model is a trade-off process among the accurate representation of the woven structure of the material, the accuracy of the computed results and the capabilities and restrictions of the conversion software. To evaluate the effect of mesh density, square and rectangular elements of different sizes have been considered. The rectangular elements, although have an increased aspect ratio, allow for a more accurate modeling of the bundle’s curvature. Table 1 provides the details of the elements. Figure 3 shows the square and rectangular FE meshes of the reduced geometry of the 3D woven composite.

The evaluation of the different meshes was done in terms of back-face velocity measurements. The corresponding back-face velocity diagrams are compared in Figure 4. As shown, there is a minor effect of the mesh density. The shock wave velocity, derived from the shock front's arrival time, remains constant, indicating that the through-the-thickness stiffness remains unaffected by changes in the mesh size. The rectangular meshes R1 and R2 exhibit a marginal increase in the first reflected pulse, while the initial peak remains unaltered. In mesh SQ3, in some boundaries of the bundles, isolated elements have been created as can be seen in Figure 5, which is a known issue in coarse voxel meshes [32]. Considering these factors and the absence of substantial deviations noted in the mesh density study, the selected mesh for the simulations is the R1, which results to a final model comprised of 6.6 million elements. Figure 6 depicts the final FE mesh of the wrap and binder bundles.

Another crucial parameter in the accuracy of the FE model is the detailed description of the fiber bundles’ geometry. This depends on the proper description of material directions, particularly for the binder yarns that present high curvature and, in some regions, appear to be almost vertical. The description of the material direction in the model requires adding a curve that follows the geometry of the binders, for each binder separately, and rotating the local coordinate system of each according to the angle defined by the position of the element on the curve. Figure 7 illustrates a highly complex binder yarn and its respective material direction at each element.

2.2. Homogenized FE Model

The homogenized FE model was created by means of the LS-PrePost also using ELFORM 1 constant stress solid elements. The element shape and size were selected accordingly so as to enable a direct comparison with the detailed model. The FE mesh of the homogenized model is shown in Figure 8.

2.3. Material Models

2.2.1. Fiber Bundles

The fiber bundles were modeled as a high fiber content unidirectional composite. In order to model the composite bundles in LS-Dyna, MAT_162 was used. This material model is characterized by strain rate dependent constants and a progressive damage formulation. Both unidirectional and woven composite configurations can be applied within the framework of this constitutive model. To describe the failure mechanisms of composite materials, a variety of failure criteria can be incorporated into the MAT_162 formulation. Specifically, the model can simulate fiber failure, matrix damage, and delamination phenomena under various loading scenarios, including opening, closure, and sliding modes of failure. The degradation of material properties is accounted for by using a damage parameter in the stiffness matrix, which controls the progressive reduction in material strength and stiffness. A more detailed explanation of the material model can be found in [33].

The homogenized mechanical properties of the fiber bundles were determined using the the Chamis micromechanical equations [34]. The elastic modulus of the composite in the longitudinal direction was calculated by

$$E_{11}=V_f{E}_{f11}+V_m{E}_m$$

(3)

where $E_{f11}$ and $E_m$are the elastic moduli of the fiber in the longitudinal direction and the matrix respectively, and $V_f$ and $V_m$ indicate the fiber and matrix content. The elastic moduli in the transverse and normal directions are considered equal and were calculated by

$$E_{22}=E_{33}=\frac{E_m}{1-\sqrt{V_f\left(1-\frac{E_m}{E_{f22}}\right)}}$$

(4)

where $E_{f22}$ is the fiber’s transverse elastic modulus. The shear moduli were derived by

$$G_{12}=G_{13}=\frac{G_m}{1-\sqrt{V_f\left(1-\frac{G_m}{G_{f12}}\right)}}$$

(5)

$$G_{23}=\frac{G_m}{1-\sqrt{V_f\left(1-\frac{G_m}{G_{f23}}\right)}}$$

(6)

where $G_m$ indicates the shear modulus of the matrix, and $G_{f12}$ and $G_{f23}$ the shear moduli of the fiber. Additionally, the micromechanical equations for the Poisson’s ratios are

$$v_{12}=v_{13}=V_f{v}_{f12}+V_m{v}_m$$

(7)

$$v_{23}=V_f{v}_{f23}+V_m\left(2v_m-\frac{v_{12}{E}_{22}}{E_{11}}\right)$$

(8)

where, $v_{f12}$ and $v_{f23}$ are the Poisson’s ratios for the fibers describing the planes of 1-2 and 1-3 respectively. $v_m$ is the Poisson’s ratio of the matrix and $E_{11}$, $E_{22}$ are the homogenized elastic moduli for the directions parallel and perpendicular to the fibers respectively. Finally, the equations that calculate the homogenized strength properties of the composite are the following

$$F_{1t}=V_f{F}_{ft}+V_m{F}_{mt}$$

(9)

$$F_{2t}=F_{3t}=\frac{1}{k_{\sigma }}{F}_{mt}$$

(10)

$$F_{1c}\cong \frac{G_m}{1-V_f}$$

(11)

$$F_{2c}=F_{3c}=\frac{F_{mc}}{K_c}$$

(12)

$$F_{23}=\frac{F_{ms}}{K_{\tau }}$$

(13)

$$F_{12}=F_{13}=F_{ms}$$

(14)

where, the subscript $"t"$ indicates the tensile and $"c"$ the compressive properties respectively. $F_{ft}$ is the tensile strength of the fibers and $F_{mt}$, $F_{mc}$ are the tensile and compressive strengths of the matrix respectively, $F_{ms}$ indicates the shear strength of the matrix. $K_{\sigma }$, $K_c$, $K_{\tau }$ are stress concentration factors that can be derived from literature [35]. The mechanical properties of the IM7 fibers are listed in Table 2.

Regarding the 2896 epoxy matrix, the mechanical properties were primarily derived from the data sheet of the manufacturer and are listed in Table 3. Notably, the manufacturer did not specify the shear strength. Consequently, an assumption was made, aligning it with the comparable epoxy matrix LY556, acquired from bibliographic sources. Finally, the calculated mechanical properties of the fiber bundles, considering a high fiber content of 90% ($v_f=0.9$), are listed in Table 4.

2.2.2. Matrix

The matrix was modeled using the elastoplastic material model 81 of LS-Dyna. This model applies an elasto-visco-plastic material behavior, allowing the users to define stress-strain curves that may exhibit a strain rate dependency. Also, this model accommodates the influence of damage preceding rupture, primarily based on the measure of the effective plastic strain. The stress-strain curve that is imported into the material model describes a bilinear elastic-plastic behavior. The yield stress of the 2896 epoxy matrix is not available in the literature; thus, an assumption was made based on the available data. As presented in Table 3, the known properties of the matrix are the stress and strain at fracture, as well as the elastic modulus. In order to combine them into an elastic-plastic curve an assumed yield stress needs to be calculated using the available elastic modulus. The strain for the calculation was chosen to be 30% of the fracture strain in order to obtain a yield stress that is lower than the tensile strength.

The laser-shock process is characterized by very high strain rates. Therefore, modification of material properties is required. In this context, maintaining a consistent adjustment of the material properties with strain rate across material models is essential. To achieve this, the same strain rate dependency formulations that are utilized for the properties of the composite bundles using MAT_162 are also employed for the modification of the matrix properties. Specifically, the elastic modulus and strength of the matrix material were determined according to [33] using

$$\frac{{\sigma }_d}{{\sigma }_{qs}}=1+C_{rate}^{\sigma }\mathit{ln}\left(\frac{\dot{\epsilon }}{{\epsilon }_0}\right)$$

(15)

$$\frac{E_d}{E_{qs}}=1+C_{rate}^{El}\mathit{ln}\left(\frac{\dot{\epsilon }}{{\epsilon }_0}\right)$$

(16)

where a dynamic property, that is denoted by the subscript d, is calculated at a specified strain rate $\dot{\epsilon }$ using the quasi-static property, denoted by the subscript qs, that has been measured at the reference strain rate $\dot{{\epsilon }_0}$. Evidently, there is a lack of high strain rate experiments for the required epoxy matrix. In order to calculate the strain rate coefficients, strain rate dependent properties of a similar, well-studied, epoxy resin were utilized. The data for the LY 556 material that were used to correlate the increase in the elastic modulus and strength with the strain rate are listed in Table 5. The least squares methodology was applied for the fitting of the data and the determination of the two Crates. The calculated values are 0.05 and 0.07 for the $C_{rate}^{\sigma }$ and $C_{rate}^{El}$ respectively. The comparison of the experimental and computed values of the elastic modulus and strength with the strain rate is presented in Figure 9.

The previously mentioned methodology has given the Crates that were used to calculate the modified stress-strain curve of the matrix at a strain rate of 10⁶/s, which is the developed strain rate in the laser-shock process. The modified stress-strain curve is presented in Figure 10.

2.2.3. Homogenized Model

The homogenized model also employed MAT_162, incorporating the composite's homogenized properties in the formulation of a standard woven composite that is provided by the material model. These properties, that are well-tuned for laser shock applications, have been effectively validated against experimental back-face velocity measurements. Furthermore, accurate properties of the 3D woven composite at high strain rates were obtained through a comprehensive experimental campaign using the Split-Hopkinson apparatus in [10].

2.4. Laser Shock Modeling

In the laser shock experimental process, a protective, sacrificial layer is used on the top surface of the CFRP specimen. In this context, a thin aluminum layer was modeled on top of the 3D woven material, both in the detailed and the homogeneous models. The shock loading was modeled as a pressure on a circular spot with a diameter of 4 mm as shown in Figure 11. The pressure loading is following the normalized temporal profile shown in Figure 12, and the maximum pressure was calculated by [41]

$$P_{max}=0.01\sqrt{\frac{a}{2a+3}}\sqrt{Z{I}_0}$$

(17)

where, I₀ (GW/cm²) is the laser intensity, a is the part of the energy being used for the ionization, Z (g cm^-2/s^-1) is the relative acoustic impedance calculated by

$$Z=2\frac{Z_1{Z}_2}{Z_1+Z_2}$$

(18)

where, $Z_1$ and $Z_2$ are the acoustic impedance of the material and the confinement respectively.

As due to the inhomogeneity of the 3D woven composite in the through-the-thickness direction a different shock wave response is expected at different normal paths of the wave, two spot locations were modeled: one causing a shock wave that intersects the fiber bundles and one causing a shock wave that passes mainly through the matrix areas.

3. Results

3.1. Shock Wave Propagation

A step-by-step comparative analysis of the shock wave propagation was conducted between the detailed and homogenized models, as depicted in Figure 13. When observing the normal stress σ_z stress at a cross-section during mid-loading, it is apparent that the behavior of the shock front and the subsequent release wave, resulting from the surface unloading, align closely in both simulations. Notably, in the detailed model the shock front is confined within a smaller area, likely due to the energy dissipation from numerous small reflections during material alterations. This effect tends to diminish the σ_z stress at the edges of the shock front, where its values are generally lower. Furthermore, at t = 4 x 10^-6 s, where the two release waves, one traveling in the shock front's direction and the other occurring upon the shockwave's reflection at the material's back-face, converge, both models effectively depict the stress field. The primary difference between the models lies in the computed values of smaller stresses at the edges and tail of the shock wave.

In addition, the shock wave propagation was evaluated for both models by means of the back-face velocity computed at a singular node placed at the center of the loading region. The corresponding diagrams are compared in Figure 14. In the detailed model, spot 1 corresponds to the fiber bundles path and spot 2 to the matrix path. A similar wave arrival time is observed for both models. This is an indication of the efficient simulation of the elastic behavior by the models and more specifically, of the efficient modeling of the macroscopic through thickness properties of the 3D woven composite material. The two velocity curves predicted by the detailed model exhibit lower peak velocities. This can be ascribed to the increased energy dissipation of the additional reflections. In Figure 15, the shock wave propagation in the detailed model is compared for the two spots. For spot 2, additional reflections caused by the intersection of the wave with the binder bundles are visible. This explains the lower peak in the back-face velocity diagram for this spot.

3.2. Damage Accumulation

A big asset of the detailed model is the capability to predict even a minor damage, including matrix cracking between the fiber bundles and several damage modes within the bundles. Validation of the localized damage is very difficult as in most cases it is not observable in experiments. However, the insights offered by the detailed model can be proved highly advantageous in efforts aimed at preserving the composite's integrity during a disassembly process.

In the present laser shock simulations, the detailed model predicted only minor matrix cracking at the resin-rich area under the loaded surface (Figure 16), two matrix cracking modes within the bundles at the area below the loaded surface and at the opposite edge of the specimen (Figure 17 and Figure 18), and, finally, fiber damage in the bundles just below the loaded surface (Figure 19).

The homogenized model has predicted only matrix cracking, parallel and perpendicular to the fibers at the loaded area as shown in Figure 20. Contrary to the detailed model, no information is given by the homogenized model on whether the predicted damage is in the resin-rich areas or within the bundles. Moreover, the homogenized model did not predict any fiber damage, which is a significant omission.

4. Conclusions

To study in detail the laser shock loading of 3D woven composites, a detailed FE model was developed and compared with a homogenized model. Observations on the shock wave propagation were made in terms of normal stress contours and back-face velocity diagrams. Also, an in-depth examination of the induced damage was conducted in both models. From the implementation process and the numerical results, the following conclusions can be drawn.

Both models efficiently simulate the laser shock wave propagation and accurately describe the convergence of release waves in the 3D woven composite as it was revealed by the relatively good comparison between the computed back-face velocity diagrams and a corresponding experimental diagram. However, there are some considerable differences on the wave propagation between the two models. In the detailed model, more wave reflections take place, thus affecting the back-face velocity curve. Moreover, using the detailed model, the effect of spot location, which determines the shock propagation path, was revealed. When the wave intersects with more bundles more reflections take place, and the wave is delayed.

The detailed model predicted matrix cracking in the resin-rich areas and the fiber bundles and fiber failure, while the homogenized model predicted only matrix cracking whose exact location cannot be justified. Those differences are attributed to the capacity of the detailed model to compute more precisely the stresses at areas with geometrical variations, thus capturing the stress concentrations.

In terms of computational efficiency, the development of the homogenized model is much easier, while the homogenized model runs faster by 36% than the detailed model.

In conclusion, if the objective is to obtain an overview of the laser shock wave propagation process in order to adjust the laser parameters so as to create a customized delamination or a debonding, the homogenized model is sufficient. On the other hand, if the objective is to ensure that no damage is developed within the 3D woven composite and to study how the composite’s architecture affects the shock propagation process, then the detailed should be used.