**1. Introduction**

Composite materials with excellent performance are becoming more attractive in the field of high technology. The RTM manufacturing process of textile composites is widely used, which includes two main stages: the forming stage and the injection stage. The forming phase is important because it strongly influences the mechanical behavior of the final composite parts. In order to better predict the geometrical and mechanical characteristics of textile composites, different types of simulations can be distinguished at three different scales: macroscopic (the fabric), mesoscopic (the yarns), and microscopic (the fibers).

To simulate the forming process of textile composites, approaches are generally carried out at the macroscopic scale [1–4], in which the reinforcement is considered as a continuous medium. At this scale, it is possible to predict the appearance of wrinkles, which is the major defect that appears during this phase [5,6]. However, some important phenomena may appear on a smaller scale during the forming process of reinforcements [7–11]. Mesoscopic models, which consider reinforcements as a set of yarns in contact with their neighbors, are able to predict possible defects, such as the rupture, the gaping, and the local buckling of yarns. In addition, the information about deformations and orientations of yarns allows for the determination of the mechanical behavior of the final parts and the calculation of the permeability tensor of the deformed reinforcement, which can be used to simulate the

**Citation:** Wang, J.; Wang, P.; Hamila, N.; Boisse, P. Meso-Macro Simulations of the Forming of 3D Non-Crimp Woven Fabrics. *Textiles* **2022**, *2*, 112–123. https://doi.org/ 10.3390/textiles2010006

Academic Editors: Jiri Militky, Raul Fangueiro and Rajesh Mishra

Received: 15 December 2021 Accepted: 2 February 2022 Published: 11 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

resin flow in the injection stage [12–14]. Nevertheless, it is difficult to simulate the shaping of woven reinforcements at this mesoscopic scale taking into account the large number of yarns and their complex interactions, especially for 3D fabrics. In these models, one of the most critical problems is the computing time.

In this context, the objective of this paper is to develop an efficient multi-scale method that allows one to obtain the deformations and orientations of yarns for a whole composite structure during the shaping process of woven reinforcements, with a reasonable computation time. The multi-scale methods that have been developed in the literature consist of creating a passage of information between the macro and meso scales [15–20]. This passage makes it possible to transform the deformed information between the macroscopic and mesoscopic scales. In other words, the deformations and the orientations of yarns, as well as the variations of the transverse sections of yarns for the entire composite part, can be calculated from macroscopic forming simulations of woven reinforcements. This proposed method is divided into three main parts: (1) the macroscopic simulation of the woven reinforcements, (2) macro-meso embedded analysis, and (3) the mesoscopic results enriched by a local mesoscopic simulation of RVE. The objective is to present the method in a simple way, postponing the presentation of the laws of macroscopic and mesoscopic behavior to the appendix. The analysis of the number of RVEs in the mesoscopic analysis is presented as well as a comparison of the method with the experiment. The material studied here is an orthogonal 3D glass woven reinforcement (see Figure 1) [21].

**Figure 1.** Structure of the 3D non-crimp orthogonal woven fabric by X-ray tomography and its mesh.

### **2. Macroscopic Analysis of Woven Reinforcements**

This part consisted of implementing the numerical simulations of the woven reinforcements at the macroscopic scale using the PlasFib software [22], which calculation codes based on the finite element method (FEM), developed by the LaMCoS laboratory. This finite element code uses a temporal scheme in an explicit dynamic. In macroscopic approaches, woven reinforcements are considered as a continuous medium with a continuous constitutive law. In order to describe the mechanical characteristics of thick textile reinforcements during draping, a 3D hyperelastic constitutive law is applied in the macroscopic simulations [23–25]. In this hyperelastic constitutive law, six independent deformation modes of the reinforcement are considered: elongation in the warp and weft direction, transverse compaction, in-plane shear, and transverse shear in the warp and weft direction. The stress–strain relationship is determined through the derivative of the strain energy density *w*, which is the sum of the strain energy of every deformation mode (Equation (1)). Then, the second Piola–Kirchhoff stress tensor *S* = can be calculated by the differentiation *∂wk*/*∂I<sup>k</sup>* and the right Cauchy–Green tensor *C* = (Equation (2)). *I<sup>k</sup>* , namely *Ielong*<sup>1</sup> , *Ielong*<sup>2</sup> , *Icomp*, *Ish*, *IshT*<sup>1</sup> , and *IshT*<sup>2</sup> , is the physical invariant associated with each deformation mode. In addition, the contribution of curvature in the hexahedral finite elements is taken into account by adding a local flexural stiffness [25].

$$w = w\_{\rm along1} \left( I\_{\rm along1} \right) + w\_{\rm along2} \left( I\_{\rm long} \right) + w\_{\rm comp} \left( I\_{\rm comp} \right) + w\_{\rm sh} \left( I\_{\rm sh} \right) + w\_{\rm sltT1} \left( I\_{\rm sltT1} \right) + w\_{\rm sltT2} \left( I\_{\rm sltT2} \right) \tag{1}$$

$$\underset{\stackrel{\circ}{=}}{S} = 2\frac{\partial w}{\partial \underline{\mathbf{C}}} = 2\frac{\partial w\_k}{\partial I\_k}\frac{\partial I\_k}{\partial \underline{\mathbf{C}}}\tag{2}$$

The characterization of the material parameters was carried out by different classical tests: tensile test, transverse compaction, in-plane shear, transverse shear, and bending test, by means of the inverse method with the Levenberg–Marquardt algorithm. The details of these tests have been described in [17], and the identified material parameters are given in Appendix A.

Once all the material parameters in the hyperelastic constitutive law were identified, the forming simulation of reinforcements was carried out in the PlasFib software. Macroscopic simulations were realized in two cases: three-point bending and hemispherical forming (Figures 2 and 3). Finally, in comparison with the corresponding experiments, the results obtained by the simulation were in good agreement. This was the case for the deformed mid-surface of the bending specimen in Figure 2b and for the shear angles after forming the hemisphere in Figure 3.

**Figure 2.** Comparisons between the numerical and experimental results for three-point bending (**a**) deformed geometry and (**b**) the deformation of the middle line of the reinforcement.

**Figure 3.** Comparisons between the experimental and numerical results for hemispherical forming. (**a**,**b**) In-plane shear angle in the same position. (**c**) Deformed geometry.

#### **3. Mesoscopic Analysis by Macro-Meso Embedded Approach**

Some mesoscopic models of the entire reinforcement have been used in different studies [8,16,26] to obtain the mesoscopic deformed information of woven composites; however, the size of the numerical model is large and the computational time is long. Therefore, the objective of the developed macro-meso embedded approach is to estimate the deformation of yarns for the entire reinforcement during the forming process with a reasonable computational time.

In this macro-meso embedded approach, the numerical results obtained by the macroscopic simulations of reinforcement are the basis. If the macroscopic results are relevant, they should be able to provide some basic and useful information at the mesoscopic scale. In short, the embedded mesoscopic geometry is determined from the result of the macroscopic simulation.

The macro-meso embedded approach allows one to link macroscopic simulations of thick reinforcements and mesoscopic modeling of RVEs (representative volume elements). For periodic textile reinforcements, the geometrical modeling of the RVE is first reconstructed according to X-ray tomography images based on the real geometry of the reinforcement [27]. The modeling of the entire reinforcement at the mesoscopic scale can then be generated in the initial configuration by the repetitions and translations of the RVE. Then, each mesoscopic node of the yarns in the woven reinforcement can be located or integrated into a macroscopic element (Figure 4). During deformation, the integrated mesoscopic element deforms together with the macroscopic element where this mesoscopic element is located, and its mesoscopic nodes have constant natural coordinates in the reference and real macroscopic element.

During deformation, the natural coordinates of a mesoscopic node in the corresponding macroscopic element remain identical. In other words, the material position of an embedded mesoscopic node remains constant in the macroscopic element. Therefore, the deformation obtained by the macroscopic simulations can provide a first geometry of reinforcements at the mesoscopic scale in the deformed configuration. As shown in Figures 5 and 6, the mesoscopic deformed geometries of the reinforcement can be obtained by the macro-meso embedded approach.

The mesoscopic results based on macroscopic simulations allow for the determination of the deformations and the deformed orientations of the yarns, and the voids in the elementary cell. The calculation time is very fast, from a few minutes to one hour depending on the reinforcement size.

**Figure 4.** Meso node of yarns embedded in a macro element. Reference and real frames.

**Figure 5.** Three-point bending. (**a**) Deformed geometry in experimental result. (**b**) Elongation of yarns in mesoscopic analysis. (**c**) Elongation of an extracted yarn.

**Figure 6.** Hemispherical stamping. (**a**) Deformed geometry in experimental result. (**b**) Compaction of yarns in mesoscopic analysis.

### **4. Mesoscopic Analysis Enriched by Local Mesoscopic Simulation**

*4.1. Local Mesoscopic Simulation*

As shown in Section 3, the mesoscopic analyses of the woven reinforcement obtained by the macro-meso embedded approach are efficient in providing a mesoscopic scale solution during shaping. However, this approach is based on the assumptions of continuous medium and the finite element method. The mesoscopic deformed configuration is directly obtained from the macroscopic analysis, without solving a mechanical problem. From the point of view of the constitutive law, the local equilibrium of the stresses is not attained, for example in the red zones in Figure 5. That is to say, local slippages between the yarns is not been taken into account, which can, in some cases, lead to an excessive elongation of yarns. For example, as shown in Figure 5, the local elongation in the yarns obtained by the macro-meso embedded approach can reach 11%, which is much higher than the reality. In order to overcome this difficulty, a local mesoscopic simulation based on macroscopic simulations and macro-meso embedded mesoscopic analysis is proposed.

The local mesoscopic simulation is performed locally, usually on one or more representative volume elements (RVEs) based on the results of the macro-meso embedded approach. The deformed configuration obtained by the embedded mesoscopic simulation constitutes the initial state of the local mesoscopic simulation. A relatively simple transverse isotropic hyperelastic behavior law (neo-Hookean) was used to describe the mechanical behavior of the yarns [28,29]. The material parameters used are given in Appendix B. Consequently, stress equilibrium was achieved within the yarns and the local mesoscopic simulation eliminated the phenomena of the excessive elongation of yarns. The results obtained by the local mesoscopic simulation are compared with the results obtained by the macro-meso embedded approach in Figures 7 and 8. In both cases, the spurious elongations obtained in the yarns by the embedded calculation were reduced to small values after the mesoscopic calculation on a RVE both in the bending case (Figure 7) and in the hemispherical forming case (Figure 8). In addition, the calculation time of each step in this multi-scale method is shown in Table 1.

**Figure 7.** Influence of the enriched local mesoscopic simulation on the elongation of yarns in threepoint bending. (**a**) Macro-meso embedded analysis. (**b**) Macro-meso embedded method with local mesoscopic simulation analysis. (**c**) Change in elongation for the same elements.

**Figure 8.** Influence of the enriched local mesoscopic simulation on the elongation of yarns in hemispherical forming. (**a**) Macro-meso embedded analysis. (**b**) Macro-meso embedded method with local mesoscopic simulation analysis. (**c**) Change in elongation for the same elements.

### *4.2. Comparison with Experimental Results*

Finally, the mesoscopic numerical results are compared with the experimental results and the numerical results with three RVEs for three-point bending. The deformations of the longitudinal and transverse sections of the yarns were observed and analyzed by a microscope.

In comparison with the experimental results (Figure 9b), it can be seen that the numerical result (Figure 9a) based on the multi-scale method was in good agreement with the experimental result relating to the deformed geometry and the outlines of the yarns.

**Figure 9.** Deformed geometry of RVE in the same position obtained by: (**a**) simulation, (**b**) experiment, and (**c**) comparison of the outlines of binder yarns and warp yarns.

#### *4.3. Influence of the Number of RVEs*

The local mesoscopic simulation was performed on a single RVE, and the boundary conditions from the macro-meso embedded approach were applied on the RVE edges. The influence of the number of RVEs was not taken into account. Therefore, to analyze the edge effects on the mesoscopic results, another simulation was carried out on three successive RVEs for three-point bending, in which the middle RVE had the same location as the single RVE. The results shown in in Figures 10 and 11 indicate that the number of RVEs had only a very weak influence on the mesoscopic results.

**Figure 10.** Influence of the number of RVEs on the meso results.

**Figure 11.** Comparisons of yarn elongation for ten elements in the same positions between the meso results of one RVE and threes RVEs: (**a**) with the result of Macro-Meso embedded approach, (**b**) without the result of Macro-Meso embedded approach.


**Table 1.** Calculation time of three-point bending and hemispherical forming.

#### **5. Conclusions**

The deformation information of woven reinforcements on the mesoscopic scale is essential to evaluate the mechanical properties of the final composite part, as well as to simulate resin flow during the next injection step. In this paper, the developed macro-meso embedded approach allows one to obtain the deformed configuration at the mesoscopic scale of woven reinforcements during the forming process. However, it has drawbacks, especially the excessive elongation of the yarns. To overcome this problem, a local mesoscopic simulation that takes into account the local slippage of yarns can be carried out by defining specific boundary conditions. This step requires a reasonable computational cost. The results of one RVE were satisfactory in comparison with the experimental results and the numerical results with three RVEs. Furthermore, the comparison with the full meso simulation at a large scale was interesting. Another perspective to validate the quality of the multi-scale approach consists in using the Hill–Mandell macro-homogeneity condition to ensure that the same amount of energy is dissipated in the RVE and the corresponding part of material in the macro-model.

**Author Contributions:** J.W.: investigation, methodology, software, validation, visualization, and writing—original draft. P.W.: investigation. N.H.: software and supervision. P.B.: conceptualization, methodology, supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The 3D hyperelastic anisotropic constitutive law is applied in the macroscopic simulations. The second Piola–Kirchhoff stress tensor *S* = can be calculated by the differentiation *∂wk*/*∂I<sup>k</sup>* and the right Cauchy–Green tensor *C* = (see Equation (2) in Section 2). For each deformation mode, the physical invariant *I<sup>k</sup>* described in Section 2 is defined by the classical invariants (Equation (A1)), and the strain energy density *w<sup>k</sup>* can be expressed in polynomial form (Equation (A2)). In order to take into consideration the curvature of the fibers, an independent bending stiffness of fibers χ is added to calculate the bending moment *Mbend* by Equation A3. The material parameters *ki*, *D*0, and *D*<sup>1</sup> are characterized by the corresponding mechanical tests (Table A1).

$$\begin{aligned} I\_{\text{elongs}} &= \ln(\sqrt{I\_{4a}}) \text{ ( $a = 1, 2$ )} & I\_{\text{comp}} &= \frac{1}{2} \ln\left(\frac{I\_3}{I\_{41}I\_{42}(1 - I\_{cp}^2)}\right) \\\ I\_{\text{cp}} &= \frac{I\_{412}}{\sqrt{I\_{41}I\_{42}}} = \sin(\gamma\_{12}) & I\_{\text{clta}} &= \frac{I\_{43}}{\sqrt{I\_{4a}I\_{43}}} = \sin(\gamma\_{a3}) \end{aligned} \tag{A1}$$

$$w = \begin{cases} \begin{array}{c} \sum\_{i=1}^{6} k\_{i} I^{i} \\ i=1 \\ \sum\_{i=1}^{6} k\_{i+6} I^{i} & \text{if } I \le 0 \end{array} \end{cases} \tag{A2}$$

$$M\_{bend} = \begin{cases} (D\_0 - D\_1|\chi|)\chi & \text{if } \chi < \frac{D\_0}{2D\_1} \\ & \frac{D\_0}{2}\chi & \text{if not} \end{cases} \tag{A3}$$


**Table A1.** Material parameters identified in the macroscopic simulations (k<sup>i</sup> in MPa; D<sup>0</sup> in Nmm; in Nmm<sup>2</sup> ).

#### **Appendix B**

In the hyperelastic transverse isotropic (Neo-Hookean) constitutive law, the potential is considered to consist of two potentials [28]: the isotropic potential *wiso*, concerning the longitudinal behavior of a yarn, and the potential *wtrans*, concerning the mechanical behavior in the transverse plane of a yarn (Equation (A4)).

$$\begin{cases} w = w\_{\rm iso}(I\_1, I\_2, I\_3) + w\_{\rm trans}(I\_1, I\_2, I\_3, I\_4, I\_5) \\ \qquad w\_{\rm iso} = \frac{1}{2}\mu(I\_1 - 3) - \mu\ln(I) + \frac{1}{2}\lambda(I - 1)^2 \\ \quad w\_{\rm trans} = \left[\mathfrak{a} + \beta\ln(I) + \gamma(I\_4 - 1)\right](I\_4 - 1) - \frac{1}{2}\mathfrak{a}(I\_5 - 1) \end{cases} \tag{A4}$$

For simplification, the parameters λ, µ, α, β and γ are defined by four material parameters in Equation (A5): transverse modulus of elasticity E, longitudinal modulus of elasticity EA, Poisson's ratio ν, and longitudinal shear modulus GA, respectively. The material parameters used in the mesoscopic simulation are given in Table A2 (details in [29]).

$$\begin{aligned} \lambda &= \frac{E(\nu + m^2)}{m(1+\nu)} \\ \mu &= \frac{E}{2(1+\nu)} \end{aligned} $$

$$\begin{aligned} \alpha &= \mu - G\_A \\ \beta &= \frac{E\nu^2(1-n)}{4m(1+\nu)} \\ \gamma &= \frac{E\_A(1-\nu)}{8m} - \frac{\lambda + 2\mu}{8} + \frac{a}{2} - \beta \\ m &= 1 - \nu - 2m\nu^2 \end{aligned} \tag{A5}$$

$$n = E\_A / E$$


**Table A2.** Material parameters identified in the mesoscopic simulations (in MPa).
