$$E_{\text{mod}} = \underbrace{E_c^0}_{E_M} + \frac{\zeta E_c^\zeta + \zeta^2 E_c^{\zeta\zeta} + E_e}{E_B} + E_c^{\text{hg}}. \quad (1)$$

Here, $E_c^0$ is the compatible strain in the element-center which is used as the membrane strain $E_M$. The tensors $E_c^\zeta$ and $E_c^{\zeta\zeta}$ denote the first and second derivatives of $E_c$ with respect to the isoparametric thickness coordinate $\zeta$. They are used with the enhanced strain $E_e$ as the bending strain $E_B$. The remaining parts of the Taylor expansion $E_c^{\text{hg}}$ describe the in-plane deformation and are used to prevent hourglass deformation.

Analogously, the stress $S$ can be split into (cf. [13–15])

$$S = S_M^*(E_M) + S_B^*(E_B) + S^{\text{hg}}(E_c^{\text{hg}}). \quad (2)$$

The stresses $S_M^*$ and $S_B^*$ describe the membrane and the bending stresses that are related to the out-of-plane integration. The hourglass stress $S^{\text{hg}}$ solely depends on $E_c^{\text{hg}}$ and is used for the hourglass stabilization. The combination of Eq. 1 and Eq. 2 enables the membrane-bending decoupling within the solid-shell element. All stresses in Eq. 2 are modeled purely elastic by Schäfer et al [13].

*Viscoelastic material law.* To model the rate-dependency, the Voigt-Kelvin approach is considered. It is based on the parallel connection of the elastic part with a viscous part. Thus, the second Piola-Kirchhoff stress $S$ using the Voigt-Kelvin model is calculated by [23]

$$S = S_{\text{elast}} + S_{\text{visc}}. \quad (3)$$

The elastic part $S_{\text{elast}}$ is evaluated using the hyperelastic St. Venant-Kirchhoff model [24]

$$S_{\text{elast}} = \mathbb{C}[E]. \quad (4)$$

Here $\mathbb{C}$ is the stiffness tensor and $E$ the Green-Lagrange strain. To evaluate the viscous part $S_{\text{visc}}$, the isotropic Newton model is used. The viscous law is given by [25]

$$\sigma_{\text{visc}} = 2\eta D. \quad (5)$$

Here $\eta$ is the viscosity, $D$ the symmetric part of the velocity gradient, respectively the deformation rate tensor, and $\sigma_{\text{visc}}$ the resulting Cauchy stress. The deformation rate tensor $D$ can be calculated dependent on $\dot{E}$ using [23]

$$D = F^{-\text{T}} \dot{E} F^{\text{T}}. \quad (6)$$

To apply Eqs. 3 to 6 on the solid-shell element, the Green-Lagrange strain tensor $E$ is used for the elastic part (Eq. 3) and the rate $\dot{E}$ for the viscous part (Eq. 5 and Eq. 6). In Eq. 6 the deformation gradient $F$ is required to obtain $D$. However, within the solid-shell element the total Green-Lagrange strain tensor is modified using the EAS- and the ANS-method. Consequently, the deformation gradient $F$ that is derived solely from the nodal displacements $u_i$ does not necessarily describe the same deformation state as $E_{\text{mod}}$. Due to the symmetry of strain tensors, it is impossible to derive $F$ directly from $E$. Thus, a suitable approximation $F_{\text{tec}}$ for $F$ in dependency of $E_M$, respectively $E_B$, has to be defined to evaluate Eq. 6. Hauptmann et al. [26] presented a procedure to calculate a deformation gradient $F_{\text{tec}}$ from a modified Green-Lagrange strain tensor $E_{\text{tec}}$ which is adopted in this work. For this, the polar decomposition of $F_{\text{tec}}$ is performed
$$E_{\text{mod}} = \underbrace{E_c^0}_{E_M} + \frac{\zeta E_c^\zeta + \zeta^2 E_c^{\zeta\zeta} + E_e}{E_B} + E_c^{\text{hg}}. \quad (1)$$

Here, $E_c^0$ is the compatible strain in the element-center which is used as the membrane strain $E_M$. The tensors $E_c^\zeta$ and $E_c^{\zeta\zeta}$ denote the first and second derivatives of $E_c$ with respect to the isoparametric thickness coordinate $\zeta$. They are used with the enhanced strain $E_e$ as the bending strain $E_B$. The remaining parts of the Taylor expansion $E_c^{\text{hg}}$ describe the in-plane deformation and are used to prevent hourglass deformation.

Analogously, the stress $S$ can be split into (cf. [13–15])

$$S = S_{\text{M}}^{*}(E_{\text{M}}) + S_{\text{B}}^{*}(E_{\text{B}}) + S^{\text{hg}}(E_{\text{c}}^{\text{hg}}). \quad (2)$$

The stresses $S_M^*$ and $S_B^*$ describe the membrane and the bending stresses that are related to the out-of-plane integration. The hourglass stress $S^{\text{hg}}$ solely depends on $E_c^{\text{hg}}$ and is used for the hourglass stabilization. The combination of Eq. 1 and Eq. 2 enables the membrane-bending decoupling within the solid-shell element. All stresses in Eq. 2 are modeled purely elastic by Schäfer et al [13].

*Viscoelastic material law.* To model the rate-dependency, the Voigt-Kelvin approach is considered. It is based on the parallel connection of the elastic part with a viscous part. Thus, the second Piola-Kirchhoff stress $S$ using the Voigt-Kelvin model is calculated by [23]

$$S = S_{\text{elast}} + S_{\text{visc}}. \quad (3)$$

The elastic part $S_{\text{elast}}$ is evaluated using the hyperelastic St. Venant-Kirchhoff model [24]

$$S_{\text{elast}} = \mathbb{C}[E]. \quad (4)$$

Here $\mathbb{C}$ is the stiffness tensor and $E$ the Green-Lagrange strain. To evaluate the viscous part $S_{\text{visc}}$, the isotropic Newton model is used. The viscous law is given by [25]

$$\sigma_{\text{visc}} = 2\eta D. \quad (5)$$

Here $\eta$ is the viscosity, $D$ the symmetric part of the velocity gradient, respectively the deformation rate tensor, and $\sigma_{\text{visc}}$ the resulting Cauchy stress. The deformation rate tensor $D$ can be calculated dependent on $\dot{E}$ using [23]

$$D = F^{-\text{T}} \dot{E} F^{\text{T}}. \quad (6)$$

To apply Eqs. 3 to 6 on the solid-shell element, the Green-Lagrange strain tensor $E$ is used for the elastic part (Eq. 3) and the rate $\dot{E}$ for the viscous part (Eq. 5 and Eq. 6). In Eq. 6 the deformation gradient $F$ is required to obtain $D$. However, within the solid-shell element the total Green-Lagrange strain tensor is modified using the EAS- and the ANS-method. Consequently, the deformation gradient $F$ that is derived solely from the nodal displacements $u_i$ does not necessarily describe the same deformation state as $E_{\text{mod}}$. Due to the symmetry of strain tensors, it is impossible to derive $F$ directly from $E$. Thus, a suitable approximation $F_{\text{tec}}$ for $F$ in dependency of $E_M$, respectively $E_B$, has to be defined to evaluate Eq. 6. Hauptmann et al. [26] presented a procedure to calculate a deformation gradient $F_{\text{tec}}$ from a modified Green-Lagrange strain tensor $E_{\text{tec}}$ which is adopted in this work. For this, the polar decomposition of $F_{\text{tec}}$ is performed
$$F_{\text{tec}} = RU_{\text{tec}}. \qquad (7)$$

The stretch tensor $U_{\text{tec}}$ is a symmetric second-order tensor. It can be derived from $E_{\text{tec}}$ using

$$(U_{\text{tec}})^2 = 2E_{\text{tec}} + I. \qquad (8)$$

Here, $I$ is the identity for second-order tensors. The tensor $R$ in Eq. 7 is an orthogonal second-order tensor. It cannot be obtained from $E_{\text{tec}}$. Thus, it is calculated from the polar decomposition of the deformation gradient $F$ that corresponds to the displacements $u_i$ of the nodes within the element. Finally, $F_{\text{tec}}$ can be calculated using Eq. 7, and subsequently the viscous stress $\sigma_{\text{visc}}$ can be obtained from Eq. 5 and Eq. 6 with $F = F_{\text{tec}}$.

To account for rate-dependent material behavior within the solid-shell element, Eqs. 3 to 8 are applied separately as constitutive equations for the membrane and the bending behavior. In summary, Eq. 2 reads

$$S = \frac{S_M^*(E_M, \dot{E}_M) + S_B^*(E_B, \dot{E}_B)}{S^*} + S^{\text{hg}}(E_c^{\text{hg}}) \qquad (9)$$

for the considered viscoelastic material law. Note that the hourglass stabilization is modeled purely elastic and thus can be adopted from Schäfer et al. [13].

*Estimation of timestep increment.* The implementation in ABAQUS/EXPLICIT requires the definition of an approximation for the upper bound of the stable time increment $\Delta t$. Using the stability limit $\Delta t_c$ obtained from Courant, Friedrichs and Lewy [27] and the safety factor $\alpha \in (0,1]$ the upper bound for $\Delta t$ is obtained by

$$\Delta t = \alpha \Delta t_c \le \alpha \frac{L_e}{c_d}, \qquad (10)$$

where $L_e$ remarks the smallest geometric element dimension and $c_d$ the dilatational wave speed. The dilatational wave speed $c_d$ is dependent on the material law. Since solid-shell elements are thin structures, $L_e$ and subsequently $\Delta t_c$ are very small. For an isotropic elastic material law $c_d$ can be computed using [28]

$$c_d = \sqrt{\frac{\lambda+2\mu}{\rho}}, \qquad (11)$$

where $\lambda$ and $\mu$ denote the Lamé constants and $\rho$ the mass density of an isotropic linear elastic material. However, in this work a viscoelastic material law is considered which is described by the Voigt-Kelvin model in Eq. 3. As this material is rate-dependent, it is not appropriate to assume that the corresponding dilatational wave speed $c_d$ from Eq. 11 is constant. Thus, it is possible that the calculated upper bound $\Delta t$ for the stable timestep is estimated too large which leads to unstable iteration schemes. Hence, it is beneficial to calculate an effective dilatational wave speed $\hat{c}_d$ which is dependent on resulting stresses $S$ and strains $E$. This enables the calculation of a stable time increment $\Delta t$ for arbitrary material laws. To calculate an effective dilatational wave speed $\hat{c}_d$, the compression modulus $\hat{R}$ and shear modulus $\hat{\mu}$ for an effective isotropic linear elastic material are calculated using the hyperelastic St. Venant-Kirchhoff model [24]

$$3\hat{R}(E^\circ \cdot E^\circ) = S^\circ \cdot E^\circ \quad \text{and} \quad 2\hat{\mu}(E' \cdot E') = S' \cdot E'. \qquad (12)$$
$$F_{\text{tec}} = RU_{\text{tec}}. \qquad (7)$$

The stretch tensor $U_{\text{tec}}$ is a symmetric second-order tensor. It can be derived from $E_{\text{tec}}$ using

$$(U_{\text{tec}})^2 = 2E_{\text{tec}} + I. \qquad (8)$$

Here, $I$ is the identity for second-order tensors. The tensor $R$ in Eq. 7 is an orthogonal second-order tensor. It cannot be obtained from $E_{\text{tec}}$. Thus, it is calculated from the polar decomposition of the deformation gradient $F$ that corresponds to the displacements $u_i$ of the nodes within the element. Finally, $F_{\text{tec}}$ can be calculated using Eq. 7, and subsequently the viscous stress $\sigma_{\text{visc}}$ can be obtained from Eq. 5 and Eq. 6 with $F = F_{\text{tec}}$.

To account for rate-dependent material behavior within the solid-shell element, Eqs. 3 to 8 are applied separately as constitutive equations for the membrane and the bending behavior. In summary, Eq. 2 reads

$$S = \frac{S_M^*(E_M, \dot{E}_M) + S_B^*(E_B, \dot{E}_B)}{S^*} + S^{\text{hg}}(E_c^{\text{hg}}) \qquad (9)$$

for the considered viscoelastic material law. Note that the hourglass stabilization is modeled purely elastic and thus can be adopted from Schäfer et al. [13].

*Estimation of timestep increment.* The implementation in ABAQUS/EXPLICIT requires the definition of an approximation for the upper bound of the stable time increment $\Delta t$. Using the stability limit $\Delta t_c$ obtained from Courant, Friedrichs and Lewy [27] and the safety factor $\alpha \in (0,1]$ the upper bound for $\Delta t$ is obtained by

$$\Delta t = \alpha \Delta t_c \le \alpha \frac{L_e}{c_d}, \qquad (10)$$

where $L_e$ remarks the smallest geometric element dimension and $c_d$ the dilatational wave speed. The dilatational wave speed $c_d$ is dependent on the material law. Since solid-shell elements are thin structures, $L_e$ and subsequently $\Delta t_c$ are very small. For an isotropic elastic material law $c_d$ can be computed using [28]

$$c_d = \sqrt{\frac{\lambda+2\mu}{\rho}}, \qquad (11)$$

where $\lambda$ and $\mu$ denote the Lamé constants and $\rho$ the mass density of an isotropic linear elastic material. However, in this work a viscoelastic material law is considered which is described by the Voigt-Kelvin model in Eq. 3. As this material is rate-dependent, it is not appropriate to assume that the corresponding dilatational wave speed $c_d$ from Eq. 11 is constant. Thus, it is possible that the calculated upper bound $\Delta t$ for the stable timestep is estimated too large which leads to unstable iteration schemes. Hence, it is beneficial to calculate an effective dilatational wave speed $\hat{c}_d$ which is dependent on resulting stresses $S$ and strains $E$. This enables the calculation of a stable time increment $\Delta t$ for arbitrary material laws. To calculate an effective dilatational wave speed $\hat{c}_d$, the compression modulus $\hat{R}$ and shear modulus $\hat{\mu}$ for an effective isotropic linear elastic material are calculated using the hyperelastic St. Venant-Kirchhoff model [24]

$$3\hat{R}(E^\circ \cdot E^\circ) = S^\circ \cdot E^\circ \quad \text{and} \quad 2\hat{\mu}(E' \cdot E') = S' \cdot E'. \qquad (12)$$
The operators (·)° and (·)′ denote the spherical and deviatoric parts of second-order tensors. Finally, the effective Lamé constants $\hat{\lambda}$ and $\hat{\mu}$ can be obtained from $\hat{K}$ and $\hat{\mu}$ in Eq. 12 and subsequently be substituted into Eq. 11 to calculate the effective dilatational wave speed $\hat{c}_d$. To calculate the effective dilatational wave speed $\hat{c}_d$ for the solid-shell element, the stress $\mathbf{S}^* = \mathbf{S}_M^* + \mathbf{S}_B^*$ and the Green-Lagrange strain $\mathbf{E}^* = \mathbf{E}_M + \mathbf{E}_B$ are used.

Note, that no effective material constants $\hat{K}$ and $\hat{\mu}$ can be estimated from Eq. 12 if e.g., $\mathbf{E}^* = \mathbf{0}$. Hence, the provided time step estimation is applied only when the estimated time step using Eq. 12 is smaller than the time step derived from the purely elastic case. Otherwise, the estimated time step from the purely elastic case is applied.

## Numerical Study

To demonstrate the rate-dependent material law based on Eq. 9, the solid-shell element is loaded with simple tensile tests and cantilever-bending tests. In both tests the beam has the length $l = 30$ mm, the width $b = 5$ mm and the thickness $t = 0.3$ mm. In total, six elements are used to discretize the beam with only one element layer in thickness-direction $e_3$. The geometry of the beam and the load cases are displayed in Fig. 1. The beams for both load cases are fixed on the left-hand side while the load is applied on the tip of the beam. The tensile test (Fig. 1a) is performed displacement-controlled. The maximum tension applied in the tensile test (Fig. 1a) is $u_{M,max} = 10$ mm. The cantilever-bending test (Fig. 1b) is performed load-controlled. The maximum load applied in the cantilever-bending test (Fig. 1b) is $F_{B,max} = 2.5$ mN. A smoothed amplitude is used to ramp-up the load-amplitude from zero to the maximum value for both coupon-tests.

<figure><img src="image_435.png" /><figcaption>Figure 1: Geometries and load cases for the tensile test (a) and the cantilever test (b).</figcaption></figure>

Since this work focuses on the rate-dependent (viscous) part of the implemented material law, the elastic part of the material law (Eq. 4) is kept constant and equal for both the membrane and bending behavior throughout this work. It is described by the anisotropic stiffness tensor in Voigt notation

$$ C_M = C_B = \begin{bmatrix} 100 & 0 & 0 & 0 & 0 & 0 \\ 0 & 10 & 0 & 0 & 0 & 0 \\ 0 & 0 & 10 & 0 & 0 & 0 \\ 0 & 0 & 0 & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 10 & 0 \\ 0 & 0 & 0 & 0 & 0 & 10 \end{bmatrix} \mathbf{V}_{\xi} \otimes \mathbf{V}_{\eta} \text{ MPa.} \qquad (12) $$

where $\mathbf{V}_{\xi}$ denote the basis dyads of the Voigt notation.

*Membrane-bending decoupling.* To demonstrate that the membrane-bending decoupling still holds for the viscoelastic material law, simple tensile and cantilever tests are performed. For the elastic part, Schäfer et al. [13] showed that a membrane-bending decoupling is considered correctly with the presented solid-shell element. In the present work a similar approach is considered to demonstrate the membrane-bending decoupling for the viscous part of the viscoelastic material
The operators (·)° and (·)′ denote the spherical and deviatoric parts of second-order tensors. Finally, the effective Lamé constants $\hat{\lambda}$ and $\hat{\mu}$ can be obtained from $\hat{K}$ and $\hat{\mu}$ in Eq. 12 and subsequently be substituted into Eq. 11 to calculate the effective dilatational wave speed $\hat{c}_d$. To calculate the effective dilatational wave speed $\hat{c}_d$ for the solid-shell element, the stress $\mathbf{S}^* = \mathbf{S}_M^* + \mathbf{S}_B^*$ and the Green-Lagrange strain $\mathbf{E}^* = \mathbf{E}_M + \mathbf{E}_B$ are used.

Note, that no effective material constants $\hat{K}$ and $\hat{\mu}$ can be estimated from Eq. 12 if e.g., $\mathbf{E}^* = \mathbf{0}$. Hence, the provided time step estimation is applied only when the estimated time step using Eq. 12 is smaller than the time step derived from the purely elastic case. Otherwise, the estimated time step from the purely elastic case is applied.

## Numerical Study

To demonstrate the rate-dependent material law based on Eq. 9, the solid-shell element is loaded with simple tensile tests and cantilever-bending tests. In both tests the beam has the length $l = 30$ mm, the width $b = 5$ mm and the thickness $t = 0.3$ mm. In total, six elements are used to discretize the beam with only one element layer in thickness-direction $e_3$. The geometry of the beam and the load cases are displayed in Fig. 1. The beams for both load cases are fixed on the left-hand side while the load is applied on the tip of the beam. The tensile test (Fig. 1a) is performed displacement-controlled. The maximum tension applied in the tensile test (Fig. 1a) is $u_{M,max} = 10$ mm. The cantilever-bending test (Fig. 1b) is performed load-controlled. The maximum load applied in the cantilever-bending test (Fig. 1b) is $F_{B,max} = 2.5$ mN. A smoothed amplitude is used to ramp-up the load-amplitude from zero to the maximum value for both coupon-tests.

<figure><img src="image_446.png" /><figcaption>Figure 1: Geometries and load cases for the tensile test (a) and the cantilever test (b).</figcaption></figure>

Since this work focuses on the rate-dependent (viscous) part of the implemented material law, the elastic part of the material law (Eq. 4) is kept constant and equal for both the membrane and bending behavior throughout this work. It is described by the anisotropic stiffness tensor in Voigt notation

$$ C_M = C_B = \begin{bmatrix} 100 & 0 & 0 & 0 & 0 & 0 \\ 0 & 10 & 0 & 0 & 0 & 0 \\ 0 & 0 & 10 & 0 & 0 & 0 \\ 0 & 0 & 0 & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 10 & 0 \\ 0 & 0 & 0 & 0 & 0 & 10 \end{bmatrix} \mathbf{V}_{\xi} \otimes \mathbf{V}_{\eta} \text{ MPa.} \qquad (12) $$

where $\mathbf{V}_{\xi}$ denote the basis dyads of the Voigt notation.

*Membrane-bending decoupling.* To demonstrate that the membrane-bending decoupling still holds for the viscoelastic material law, simple tensile and cantilever tests are performed. For the elastic part, Schäfer et al. [13] showed that a membrane-bending decoupling is considered correctly with the presented solid-shell element. In the present work a similar approach is considered to demonstrate the membrane-bending decoupling for the viscous part of the viscoelastic material
law (Eq. 5). Virtual tensile and cantilever tests (Fig. 1) are performed with a constant velocity. The viscosities $\eta_M$ and $\eta_B$ (cf. Eq. 5) for the membrane and bending behavior, respectively, are both varied between a low viscosity $\eta^{\text{low}} = 0.1 \text{ MPa} \cdot \text{s}$ and a high viscosity $\eta^{\text{high}} = 1 \text{ MPa} \cdot \text{s}$.

<figure><img src="image_436.png" /><figcaption>Figure 2: Resulting virtual material response for different combinations of high and low viscosities for both the membrane and the bending behavior.</figcaption></figure>

The results for both the tensile tests and the cantilever tests are depicted in Fig. 2. The different combinations of high and low viscosities for the membrane and bending behavior are symbolized using different colors. Since the trajectories partially overlap, colored dots are used to distinguish the different curves. The results of the tensile test in Fig. 2 indicate that the reaction force in deformation direction undergoes changes solely when there is a variation in membrane viscosity $\eta_M$. A change of the bending viscosity $\eta_B$ does not affect the result of the tensile test. Analogously, the cantilever bending tests in Fig. 2 shows that the resulting trajectory is only dependent on the bending viscosity $\eta_B$ and a change of the membrane viscosity $\eta_M$ does not affect the result of the cantilever-bending test. In summary it can be concluded that the membrane-bending decoupling still holds for the viscous material behavior.

*Demonstration of the rate-dependent membrane behavior.* To investigate the rate-dependency of the viscoelastic material law inside the membrane behavior, displacement-controlled tensile tests (Fig. 1a) are performed. The membrane viscosity is chosen equally for all tensile tests $\eta_M = 1 \text{ MPa} \cdot \text{s}$. The bending viscosity is omitted for the tensile tests since the above results (cf. Fig. 2) indicated that the bending viscosity $\eta_B$ does not affect the result of the tensile test. The results of the viscoelastic tensile tests are compared with the purely elastic case where the viscous part of the membrane behavior is neglected ($\eta_M = 0 \text{ MPa} \cdot \text{s}$).

The results of the tensile tests are displayed in Fig. 3. Here, the evolution of the reaction force $R_1$ in load direction are displayed over the displacement $u_M$ in load direction. The different colors represent different load times respectively load rates. The blue line represents the elastic case which is time independent. In the elastic case, the reaction force $R_1$ increases monotonically with increasing displacement $u_M$. Due to the geometric nonlinearity of the strain $E$, the trajectory of $R_1$ is not linear over $u_M$. In case of the viscoelastic material law, the reaction force $R_1$ increases with increasing deformation speed for a certain displacement $u_M$. This is due to the larger deformation rate tensor $D$ for larger deformation rates. This shows that the general concept of viscoelasticity is captured accurately, and a rate-dependent material law can be displayed using the viscoelastic solid-shell element. The viscous curves do not follow a monotonic pattern. However, they approach the same reaction force $R_1$ as observed in the elastic case. This is due to the smoothed load amplitude which starts and ends with very small deformation rates. Since the deformation
law (Eq. 5). Virtual tensile and cantilever tests (Fig. 1) are performed with a constant velocity. The viscosities $\eta_M$ and $\eta_B$ (cf. Eq. 5) for the membrane and bending behavior, respectively, are both varied between a low viscosity $\eta^{\text{low}} = 0.1 \text{ MPa} \cdot \text{s}$ and a high viscosity $\eta^{\text{high}} = 1 \text{ MPa} \cdot \text{s}$.

<figure><img src="image_447.png" /><figcaption>Figure 2: Resulting virtual material response for different combinations of high and low viscosities for both the membrane and the bending behavior.</figcaption></figure>

The results for both the tensile tests and the cantilever tests are depicted in Fig. 2. The different combinations of high and low viscosities for the membrane and bending behavior are symbolized using different colors. Since the trajectories partially overlap, colored dots are used to distinguish the different curves. The results of the tensile test in Fig. 2 indicate that the reaction force in deformation direction undergoes changes solely when there is a variation in membrane viscosity $\eta_M$. A change of the bending viscosity $\eta_B$ does not affect the result of the tensile test. Analogously, the cantilever bending tests in Fig. 2 shows that the resulting trajectory is only dependent on the bending viscosity $\eta_B$ and a change of the membrane viscosity $\eta_M$ does not affect the result of the cantilever-bending test. In summary it can be concluded that the membrane-bending decoupling still holds for the viscous material behavior.

*Demonstration of the rate-dependent membrane behavior.* To investigate the rate-dependency of the viscoelastic material law inside the membrane behavior, displacement-controlled tensile tests (Fig. 1a) are performed. The membrane viscosity is chosen equally for all tensile tests $\eta_M = 1 \text{ MPa} \cdot \text{s}$. The bending viscosity is omitted for the tensile tests since the above results (cf. Fig. 2) indicated that the bending viscosity $\eta_B$ does not affect the result of the tensile test. The results of the viscoelastic tensile tests are compared with the purely elastic case where the viscous part of the membrane behavior is neglected ($\eta_M = 0 \text{ MPa} \cdot \text{s}$).

The results of the tensile tests are displayed in Fig. 3. Here, the evolution of the reaction force $R_1$ in load direction are displayed over the displacement $u_M$ in load direction. The different colors represent different load times respectively load rates. The blue line represents the elastic case which is time independent. In the elastic case, the reaction force $R_1$ increases monotonically with increasing displacement $u_M$. Due to the geometric nonlinearity of the strain $E$, the trajectory of $R_1$ is not linear over $u_M$. In case of the viscoelastic material law, the reaction force $R_1$ increases with increasing deformation speed for a certain displacement $u_M$. This is due to the larger deformation rate tensor $D$ for larger deformation rates. This shows that the general concept of viscoelasticity is captured accurately, and a rate-dependent material law can be displayed using the viscoelastic solid-shell element. The viscous curves do not follow a monotonic pattern. However, they approach the same reaction force $R_1$ as observed in the elastic case. This is due to the smoothed load amplitude which starts and ends with very small deformation rates. Since the deformation
rate tensor **D** nearly vanishes for $u_M \approx u_{M,max}$ the resulting reaction force is only caused by the elastic contribution of the material law. Therefore, it can be inferred that the variation in force between the elastic curve (blue) and the viscoelastic curves (red, green, orange) is solely attributed to viscoelastic properties.

<figure><img src="image_437.png" /><figcaption>Figure 3: Resulting reaction forces for different deformation rates for the membrane behavior of the viscoelastic solid-shell element.</figcaption></figure>

*Demonstration of the rate-dependent bending behavior.* To investigate the rate-dependency of the viscoelastic material law of the bending behavior the cantilever-bending test is performed force-controlled (Fig. 1b). The bending viscosity is chosen equally for all cantilever-bending tests $\eta_B = 1 \text{ MPa} \cdot \text{s}$. The membrane viscosity is omitted for the bending tests since Fig. 2 showed that the membrane viscosity $\eta_M$ does not affect the result of the cantilever-bending test. The results of the viscoelastic cantilever-bending tests are compared with the purely elastic case. Here the viscous part of the bending behavior is also neglected ($\eta_B = 0 \text{ MPa} \cdot \text{s}$).

The results of the different cantilever-tests are displayed in Fig. 4. Here, the curves of the applied bending force $F_B$ over the resulting displacement of the tip of the beam $u_3$ are displayed. The different colors represent different load times respectively load speeds. The blue line represents the elastic case which is time independent. For elastic material behavior, the flattest curve is attained. Furthermore, an increase in load speed corresponds to a steeper trajectory of the bending force $F_B$ for the deformation $u_3$. Conversely, an increased load speed demands more bending force $F_B$ for the same tip displacement $u_3$. This is attributed to a greater deformation rate tensor **D** for higher load speeds. Consequently, it can be deduced that the material's rate-dependent behavior can be captured under bending load.

Note that none of the simulations in Fig. 4 reach a stationary state and all trajectories become nearly horizontal towards the end of the simulation time. For the viscous simulations (red, green, orange), the condition $\mathbf{D} \approx \mathbf{0}$ does not hold immediately at $F_B = F_{B,max}$ because the bending test is performed force-controlled (Fig. 1b). Thus, the elastic equilibrium, where all trajectories meet in the same point, is only reached for large simulation times. Since the major rate-dependent effects of the viscoelastic bending behavior are observable using the initial slope of the trajectories, larger simulation times are not considered in this work.
rate tensor **D** nearly vanishes for $u_M \approx u_{M,max}$ the resulting reaction force is only caused by the elastic contribution of the material law. Therefore, it can be inferred that the variation in force between the elastic curve (blue) and the viscoelastic curves (red, green, orange) is solely attributed to viscoelastic properties.

<figure><img src="image_448.png" /><figcaption>Figure 3: Resulting reaction forces for different deformation rates for the membrane behavior of the viscoelastic solid-shell element.</figcaption></figure>

*Demonstration of the rate-dependent bending behavior.* To investigate the rate-dependency of the viscoelastic material law of the bending behavior the cantilever-bending test is performed force-controlled (Fig. 1b). The bending viscosity is chosen equally for all cantilever-bending tests $\eta_B = 1 \text{ MPa} \cdot \text{s}$. The membrane viscosity is omitted for the bending tests since Fig. 2 showed that the membrane viscosity $\eta_M$ does not affect the result of the cantilever-bending test. The results of the viscoelastic cantilever-bending tests are compared with the purely elastic case. Here the viscous part of the bending behavior is also neglected ($\eta_B = 0 \text{ MPa} \cdot \text{s}$).

The results of the different cantilever-tests are displayed in Fig. 4. Here, the curves of the applied bending force $F_B$ over the resulting displacement of the tip of the beam $u_3$ are displayed. The different colors represent different load times respectively load speeds. The blue line represents the elastic case which is time independent. For elastic material behavior, the flattest curve is attained. Furthermore, an increase in load speed corresponds to a steeper trajectory of the bending force $F_B$ for the deformation $u_3$. Conversely, an increased load speed demands more bending force $F_B$ for the same tip displacement $u_3$. This is attributed to a greater deformation rate tensor **D** for higher load speeds. Consequently, it can be deduced that the material's rate-dependent behavior can be captured under bending load.

Note that none of the simulations in Fig. 4 reach a stationary state and all trajectories become nearly horizontal towards the end of the simulation time. For the viscous simulations (red, green, orange), the condition $\mathbf{D} \approx \mathbf{0}$ does not hold immediately at $F_B = F_{B,max}$ because the bending test is performed force-controlled (Fig. 1b). Thus, the elastic equilibrium, where all trajectories meet in the same point, is only reached for large simulation times. Since the major rate-dependent effects of the viscoelastic bending behavior are observable using the initial slope of the trajectories, larger simulation times are not considered in this work.
<figure><img src="image_438.png" /><figcaption>Figure 4: Resulting displacements for different load rates for the bending behavior of the viscoelastic solid-shell element.</figcaption></figure>

## Conclusion and outlook

This work presents the enhancement of a locking-free linear elastic and anisotropic solid-shell element to a viscoelastic solid-shell element for a future application in thermoforming simulation of thermoplastic tapes. The considered solid-shell element is implemented as a user-element (VUEL) in ABAQUS/EXPLICIT. The Voigt-Kelvin model is used to model both the membrane and the bending behavior of the solid-shell element. It is shown that the presented approach can preserve the membrane-bending decoupling for the viscous stresses in coupon-tests. Thus, the viscoelastic material law can be applied to both the membrane and the bending behavior independently. An alternative estimation of the stable time-increment is presented which takes viscoelastic stresses into account. Numerical studies are carried out to demonstrate the viscoelastic material law in coupon-tests. The results show that the rate-dependent material behavior is working for both the membrane and the bending-part of the solid-shell element. Thus, the general approach of the viscoelastic material behavior within the solid-shell element is suitable as a further step addressing out-of-plane phenomena in thermoforming simulations with a solid-shell element.

In this study only generic material properties were used to test the rate-dependent material laws. Future studies will extend the presented viscoelastic approach to more realistic material behavior to model the forming behavior of preimpregnated textile composites. A comparison of results of the conventional shell-elements and the solid-shell elements will be carried out to investigate the advantages of constitutive modeling the out-of-plane forming behavior.

## Acknowledgements

This work has been carried out in the DFG AI Research Unit 5339, funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 459291153. The work is also part of the Heisenberg project ‘Digitalization of fiber-reinforced polymer processes for resource-efficient manufacturing of lightweight components’, funded by the DFG (project no. 798455807141). The authors thank the German Research Foundation for its financial support.

## References

[1] F. Henning, L. Kärger, D. Dörr, F.J. Schirmaier, J. Seuffert, A. Bernath, Fast processing and continuous simulation of automotive structural composite components, Composites Science and Technology 171 (2019) 261–279. https://doi.org/10.1016/j.compscitech.2018.12.007
<figure><img src="image_449.png" /><figcaption>Figure 4: Resulting displacements for different load rates for the bending behavior of the viscoelastic solid-shell element.</figcaption></figure>

## Conclusion and outlook

This work presents the enhancement of a locking-free linear elastic and anisotropic solid-shell element to a viscoelastic solid-shell element for a future application in thermoforming simulation of thermoplastic tapes. The considered solid-shell element is implemented as a user-element (VUEL) in ABAQUS/EXPLICIT. The Voigt-Kelvin model is used to model both the membrane and the bending behavior of the solid-shell element. It is shown that the presented approach can preserve the membrane-bending decoupling for the viscous stresses in coupon-tests. Thus, the viscoelastic material law can be applied to both the membrane and the bending behavior independently. An alternative estimation of the stable time-increment is presented which takes viscoelastic stresses into account. Numerical studies are carried out to demonstrate the viscoelastic material law in coupon-tests. The results show that the rate-dependent material behavior is working for both the membrane and the bending-part of the solid-shell element. Thus, the general approach of the viscoelastic material behavior within the solid-shell element is suitable as a further step addressing out-of-plane phenomena in thermoforming simulations with a solid-shell element.

In this study only generic material properties were used to test the rate-dependent material laws. Future studies will extend the presented viscoelastic approach to more realistic material behavior to model the forming behavior of preimpregnated textile composites. A comparison of results of the conventional shell-elements and the solid-shell elements will be carried out to investigate the advantages of constitutive modeling the out-of-plane forming behavior.

## Acknowledgements

This work has been carried out in the DFG AI Research Unit 5339, funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 459291153. The work is also part of the Heisenberg project ‘Digitalization of fiber-reinforced polymer processes for resource-efficient manufacturing of lightweight components’, funded by the DFG (project no. 798455807141). The authors thank the German Research Foundation for its financial support.

## References

[1] F. Henning, L. Kärger, D. Dörr, F.J. Schirmaier, J. Seuffert, A. Bernath, Fast processing and continuous simulation of automotive structural composite components, Composites Science and Technology 171 (2019) 261–279. https://doi.org/10.1016/j.compscitech.2018.12.007
[2] P. Bussetta, N. Correia, Numerical forming of continuous fibre reinforced composite material: A review, Composites Part A: Applied Science and Manufacturing 113 (2018) 12-31.
https://doi.org/10.1016/j.compositesa.2018.07.010

[3] B. Schäfer, D. Dörr, L. Kärger, Reduced-Integrated 8-Node Hexahedral Solid-Shell Element for the Macroscopic Forming Simulation of Continuous Fibre-Reinforced Polymers, Procedia Manufacturing 47 (2020) 134-139. https://doi.org/10.1016/j.promfg.2020.04.154

[4] U. Sachs, R. Akkerman, Viscoelastic bending model for continuous fiber-reinforced thermoplastic composites in melt, Composites Part A: Applied Science and Manufacturing.
https://doi.org/10.1016/j.compositesa.2017.05.032

[5] A. Margossian, S. Bel, R. Hinterhoelzl, Bending characterisation of a molten unidirectional carbon fibre reinforced thermoplastic composite using a Dynamic Mechanical Analysis system, Composites Part A: Applied Science and Manufacturing 77 (2015) 154-163.
https://doi.org/10.1016/j.compositesa.2015.06.015

[6] S. Ropers, M. Kardos, T.A. Osswald, A thermo-viscoelastic approach for the characterization and modeling of the bending behavior of thermoplastic composites, Composites Part A: Applied Science and Manufacturing 90 (2016) 22-32. https://doi.org/10.1016/j.compositesa.2016.06.016

[7] D. Dörr, F. Henning, L. Kärger, Nonlinear hyperviscoelastic modelling of intra-ply deformation behaviour in Finite element forming simulation of continuously fibre-reinforced thermoplastics, Composites Part A: Applied Science and Manufacturing.
https://doi.org/10.1016/j.compositesa.2018.03.037

[8] D. Dörr, F.J. Schirmaier, F. Henning, L. Kärger, A viscoelastic approach for modeling bending behavior in finite element forming simulation of continuously fiber reinforced composites, Composites Part A: Applied Science and Manufacturing 94 (2017) 113-123.
https://doi.org/10.1016/j.compositesa.2016.11.027

[9] D. Dörr, T. Joppich, D. Kugele, F. Henning, L. Kärger, A coupled thermomechanical approach for finite element forming simulation of continuously fiber-reinforced semi-crystalline thermoplastics, Composites Part A: Applied Science and Manufacturing 125 (2019).
https://doi.org/10.1016/j.compositesa.2019.105508

[10] Hu Xiong, Simulation of forming, compaction and consolidation of thermoplastic composites based on solid shell elements, 2018.

[11] H. Xiong, N. Hamila, P. Boisse, Consolidation Modeling during Thermoforming of Thermoplastic Composite Prepregs, Materials (Basel, Switzerland) 12 (2019).
https://doi.org/10.3390/ma12182853

[12] P.D. Mulye, L. Morançay, C. Binetruy, S. Comas-Cardona, A. Leygue, D. Guillon, Unified numerical process modeling of forming and consolidation for thermoplastic composites with prepreg patches, Front. Mater. 10 (2023). https://doi.org/10.3389/fmats.2023.1176482

[13] B. Schäfer, D. Dörr, L. Kärger, Potential and challenges of a solid-shell element for the macroscopic forming simulation of engineering textiles, ESAFORM 2021.
https://doi.org/10.25518/esaform21.883

[14] M. Schwarze, S. Reese, A reduced integration solid-shell finite element based on the EAS and the ANS concept-Geometrically linear problems, Int. J. Numer. Meth. Engng. 80 (2009) 1322-1355. https://doi.org/10.1002/nme.2653
[2] P. Bussetta, N. Correia, Numerical forming of continuous fibre reinforced composite material: A review, Composites Part A: Applied Science and Manufacturing 113 (2018) 12-31.
https://doi.org/10.1016/j.compositesa.2018.07.010

[3] B. Schäfer, D. Dörr, L. Kärger, Reduced-Integrated 8-Node Hexahedral Solid-Shell Element for the Macroscopic Forming Simulation of Continuous Fibre-Reinforced Polymers, Procedia Manufacturing 47 (2020) 134-139. https://doi.org/10.1016/j.promfg.2020.04.154

[4] U. Sachs, R. Akkerman, Viscoelastic bending model for continuous fiber-reinforced thermoplastic composites in melt, Composites Part A: Applied Science and Manufacturing.
https://doi.org/10.1016/j.compositesa.2017.05.032

[5] A. Margossian, S. Bel, R. Hinterhoelzl, Bending characterisation of a molten unidirectional carbon fibre reinforced thermoplastic composite using a Dynamic Mechanical Analysis system, Composites Part A: Applied Science and Manufacturing 77 (2015) 154-163.
https://doi.org/10.1016/j.compositesa.2015.06.015

[6] S. Ropers, M. Kardos, T.A. Osswald, A thermo-viscoelastic approach for the characterization and modeling of the bending behavior of thermoplastic composites, Composites Part A: Applied Science and Manufacturing 90 (2016) 22-32. https://doi.org/10.1016/j.compositesa.2016.06.016

[7] D. Dörr, F. Henning, L. Kärger, Nonlinear hyperviscoelastic modelling of intra-ply deformation behaviour in Finite element forming simulation of continuously fibre-reinforced thermoplastics, Composites Part A: Applied Science and Manufacturing.
https://doi.org/10.1016/j.compositesa.2018.03.037

[8] D. Dörr, F.J. Schirmaier, F. Henning, L. Kärger, A viscoelastic approach for modeling bending behavior in finite element forming simulation of continuously fiber reinforced composites, Composites Part A: Applied Science and Manufacturing 94 (2017) 113-123.
https://doi.org/10.1016/j.compositesa.2016.11.027

[9] D. Dörr, T. Joppich, D. Kugele, F. Henning, L. Kärger, A coupled thermomechanical approach for finite element forming simulation of continuously fiber-reinforced semi-crystalline thermoplastics, Composites Part A: Applied Science and Manufacturing 125 (2019).
https://doi.org/10.1016/j.compositesa.2019.105508

[10] Hu Xiong, Simulation of forming, compaction and consolidation of thermoplastic composites based on solid shell elements, 2018.

[11] H. Xiong, N. Hamila, P. Boisse, Consolidation Modeling during Thermoforming of Thermoplastic Composite Prepregs, Materials (Basel, Switzerland) 12 (2019).
https://doi.org/10.3390/ma12182853

[12] P.D. Mulye, L. Morançay, C. Binetruy, S. Comas-Cardona, A. Leygue, D. Guillon, Unified numerical process modeling of forming and consolidation for thermoplastic composites with prepreg patches, Front. Mater. 10 (2023). https://doi.org/10.3389/fmats.2023.1176482

[13] B. Schäfer, D. Dörr, L. Kärger, Potential and challenges of a solid-shell element for the macroscopic forming simulation of engineering textiles, ESAFORM 2021.
https://doi.org/10.25518/esaform21.883

[14] M. Schwarze, S. Reese, A reduced integration solid-shell finite element based on the EAS and the ANS concept-Geometrically linear problems, Int. J. Numer. Meth. Engng. 80 (2009) 1322-1355. https://doi.org/10.1002/nme.2653
[15] M. Schwarze, S. Reese, A reduced integration solid-shell finite element based on the EAS and the ANS concept-Large deformation problems, Int. J. Numer. Meth. Engng. 85 (2011) 289-329.
https://doi.org/10.1002/nme.2966

[16] M. Pagani, S. Reese, U. Perego, Computationally efficient explicit nonlinear analyses using reduced integration-based solid-shell finite elements, Computer Methods In Applied Mechanics And Engineering 268 (2014) 141-159. https://doi.org/10.1016/j.cma.2013.09.005

[17] K. Jože, Multi-language and Multi-environment Generation of Nonlinear Finite Element Codes, Engineering with Computers 18 (2002) 312-327. https://doi.org/10.1007/s003660200028

[18] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Meth. Engng. 29 (1990) 1595-1638.
https://doi.org/10.1002/nme.1620290802

[19] J.C. Simo, F. Armero, Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes, Int. J. Numer. Meth. Engng. 33 (1992) 1413-1449.
https://doi.org/10.1002/nme.1620330705

[20] J. Simo, F. Armero, R. Taylor, Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems, Computer Methods In Applied Mechanics And Engineering 110 (1993) 359-386.

[21] T.J.R. Hughes, T.E. Tezduyar, Finite Elements Based Upon Mindlin Plate Theory With Particular Reference to the Four-Node Bilinear Isoparametric Element, J. Appl. Mech 48 (1981) 587-596. https://doi.org/10.1115/1.3157679

[22] M. Schwarze, I.N. Vladimirov, S. Reese, On the implementation of the EAS and ANS concept into a reduced integration continuum shell element and applications to sheet forming, Int J Mater Form 2 (2009) 919-922. https://doi.org/10.1007/s12289-009-0634-2.

[23] A. Bertram, Elasticity and Plasticity of Large Deformations, Springer International Publishing, Cham, 2021.

[24] T. Belytschko, W.K. Liu, B. Moran, K.I. Elkhodary, Nonlinear finite elements for continua and structures, 2nd ed., Wiley, Chichester, 2014.

[25] C.W. Macosko, Rheology: Principles, measurements, and applications, Wiley-VCH, New York NY u.a., 1994.

[26] R. Hauptmann, K. Schweizerhof, S. Doll, Extension of the 'solid-shell' concept for application to large elastic and large elastoplastic deformations, Int. J. Numer. Meth. Engng. 49 (2000) 1121-1141. https://doi.org/10.1002/1097-0207(20001130)49:9<1121:AID-NME130>3.0.CO;2-F

[27] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. 100 (1928) 32-74. https://doi.org/10.1007/BF01448839

[28] G. Cocchetti, M. Pagani, U. Perego, Selective mass scaling and critical time-step estimate for explicit dynamics analyses with solid-shell elements, Computers & Structures 127 (2013) 39-52.
https://doi.org/10.1016/j.compstruc.2012.10.021
[15] M. Schwarze, S. Reese, A reduced integration solid-shell finite element based on the EAS and the ANS concept-Large deformation problems, Int. J. Numer. Meth. Engng. 85 (2011) 289-329.
https://doi.org/10.1002/nme.2966

[16] M. Pagani, S. Reese, U. Perego, Computationally efficient explicit nonlinear analyses using reduced integration-based solid-shell finite elements, Computer Methods In Applied Mechanics And Engineering 268 (2014) 141-159. https://doi.org/10.1016/j.cma.2013.09.005

[17] K. Jože, Multi-language and Multi-environment Generation of Nonlinear Finite Element Codes, Engineering with Computers 18 (2002) 312-327. https://doi.org/10.1007/s003660200028

[18] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Meth. Engng. 29 (1990) 1595-1638.
https://doi.org/10.1002/nme.1620290802

[19] J.C. Simo, F. Armero, Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes, Int. J. Numer. Meth. Engng. 33 (1992) 1413-1449.
https://doi.org/10.1002/nme.1620330705

[20] J. Simo, F. Armero, R. Taylor, Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems, Computer Methods In Applied Mechanics And Engineering 110 (1993) 359-386.

[21] T.J.R. Hughes, T.E. Tezduyar, Finite Elements Based Upon Mindlin Plate Theory With Particular Reference to the Four-Node Bilinear Isoparametric Element, J. Appl. Mech 48 (1981) 587-596. https://doi.org/10.1115/1.3157679

[22] M. Schwarze, I.N. Vladimirov, S. Reese, On the implementation of the EAS and ANS concept into a reduced integration continuum shell element and applications to sheet forming, Int J Mater Form 2 (2009) 919-922. https://doi.org/10.1007/s12289-009-0634-2.

[23] A. Bertram, Elasticity and Plasticity of Large Deformations, Springer International Publishing, Cham, 2021.

[24] T. Belytschko, W.K. Liu, B. Moran, K.I. Elkhodary, Nonlinear finite elements for continua and structures, 2nd ed., Wiley, Chichester, 2014.

[25] C.W. Macosko, Rheology: Principles, measurements, and applications, Wiley-VCH, New York NY u.a., 1994.

[26] R. Hauptmann, K. Schweizerhof, S. Doll, Extension of the 'solid-shell' concept for application to large elastic and large elastoplastic deformations, Int. J. Numer. Meth. Engng. 49 (2000) 1121-1141. https://doi.org/10.1002/1097-0207(20001130)49:9<1121:AID-NME130>3.0.CO;2-F

[27] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. 100 (1928) 32-74. https://doi.org/10.1007/BF01448839

[28] G. Cocchetti, M. Pagani, U. Perego, Selective mass scaling and critical time-step estimate for explicit dynamics analyses with solid-shell elements, Computers & Structures 127 (2013) 39-52.
https://doi.org/10.1016/j.compstruc.2012.10.021
# Strain gradient calculation as a basis for localized roving slip prediction in macroscopic forming simulation of non-crimp fabrics

WANK Jan Paul¹,ª*, SCHÄFER Bastian¹, MITSCH Johannes¹ and KÄRGER Luise¹

¹Karlsruhe Institute of Technology (KIT), Institute of Vehicle System Technology (FAST) – Lightweight Engineering, 76131 Karlsruhe, Germany

ªjan.wank@kit.edu

**Keywords:** Non-Crimp Fabrics, Forming Simulation, Local Effects

**Abstract.** The utilization of the Finite Element Method (FEM) in forming simulation presents the possibility for a thorough examination of the deformation behavior exhibited by engineering textiles during the draping processes. In macroscopic forming simulations the relevant forming effects are depicted in a homogenized way. Slippage of fibers is an essential deformation mechanism of non-crimp fabrics (NCF). Experimentally, this is already observed at coupon level performing the bias-extension test (BET). Significant slippage occurs locally in the transition areas between shear zones with deviating shear angles. In existing macroscopic simulation approaches, roving slippage is only considered homogenized over the shear zones. A localized slip between individual fiber rovings cannot be modelled. Therefore, in the present work a neighboring element method for ABAQUS/EXPLICIT is introduced. This method uses multiple subroutines to transfer information between elements. The functionality of the neighboring method is confirmed by calculating a cross element gradient of the shear angle. The calculation of the shear angle gradient is shown in the simulation of the BET, giving rise to the transition zones which have been experimentally highlighted.

**Introduction**

Lightweight construction is a holistic development concept that aims to reduce the mass of a system while taking economic and ecological constraints into account [1]. Due to their fiber architecture continuous fiber-reinforced composites offer great potential for lightweight construction. A deep understanding of the draping process of engineering textiles is necessary to achieve optimal component design. The FEM is a proven method for simulating the forming behavior and process-induced defects of fiber-reinforced composites.

Due to their difference in fiber architecture, woven fabrics and NCF exhibit different draping behavior. While woven fabrics have an intrinsic cohesion, for NCF this is realized extrinsically through the stitching, resulting in straight fibers without undulations. Therefore, the shear behavior and thus the forming behavior of NCF is strongly dependent on the position of the stitching. Compared to bidirectional NCF the forming of unidirectional (UD-) NCF is predominantly influenced by simple shear rather than pure shear [2], leading to a more challenging draping process. The very high longitudinal stiffness of UD-NCF brings a greater degree of susceptibility for localized draping effects such as fiber gapping, slippage of fiber rovings and waviness or fiber curvature [3–6].

To determine the forming behavior of fiber-reinforced composites, commonly experimental tests are carried out. The BET is widely used to characterize the in-plane shear behavior of fabrics. For the BET, the pin-jointed net assumption is often made, which assumes that fiber crossings act as rotational knots and that fibers are inextensible [7]. This assumption is generally recognized for woven fabrics and results in a deformation with three distinct shear zones [4,7,8]. For woven fabrics, in one of the zones pure shear occurs. In the second zones only half of the shear from the

<figure><img src="image_439.png" /></figure>
# Strain gradient calculation as a basis for localized roving slip prediction in macroscopic forming simulation of non-crimp fabrics

WANK Jan Paul¹,ᵃ*, SCHÄFER Bastian¹, MITSCH Johannes¹ and KÄRGER Luise¹

¹Karlsruhe Institute of Technology (KIT), Institute of Vehicle System Technology (FAST) – Lightweight Engineering, 76131 Karlsruhe, Germany

ᵃjan.wank@kit.edu

**Keywords:** Non-Crimp Fabrics, Forming Simulation, Local Effects

**Abstract.** The utilization of the Finite Element Method (FEM) in forming simulation presents the possibility for a thorough examination of the deformation behavior exhibited by engineering textiles during the draping processes. In macroscopic forming simulations the relevant forming effects are depicted in a homogenized way. Slippage of fibers is an essential deformation mechanism of non-crimp fabrics (NCF). Experimentally, this is already observed at coupon level performing the bias-extension test (BET). Significant slippage occurs locally in the transition areas between shear zones with deviating shear angles. In existing macroscopic simulation approaches, roving slippage is only considered homogenized over the shear zones. A localized slip between individual fiber rovings cannot be modelled. Therefore, in the present work a neighboring element method for ABAQUS/EXPLICIT is introduced. This method uses multiple subroutines to transfer information between elements. The functionality of the neighboring method is confirmed by calculating a cross element gradient of the shear angle. The calculation of the shear angle gradient is shown in the simulation of the BET, giving rise to the transition zones which have been experimentally highlighted.

**Introduction**

Lightweight construction is a holistic development concept that aims to reduce the mass of a system while taking economic and ecological constraints into account [1]. Due to their fiber architecture continuous fiber-reinforced composites offer great potential for lightweight construction. A deep understanding of the draping process of engineering textiles is necessary to achieve optimal component design. The FEM is a proven method for simulating the forming behavior and process-induced defects of fiber-reinforced composites.

Due to their difference in fiber architecture, woven fabrics and NCF exhibit different draping behavior. While woven fabrics have an intrinsic cohesion, for NCF this is realized extrinsically through the stitching, resulting in straight fibers without undulations. Therefore, the shear behavior and thus the forming behavior of NCF is strongly dependent on the position of the stitching. Compared to bidirectional NCF the forming of unidirectional (UD-) NCF is predominantly influenced by simple shear rather than pure shear [2], leading to a more challenging draping process. The very high longitudinal stiffness of UD-NCF brings a greater degree of susceptibility for localized draping effects such as fiber gapping, slippage of fiber rovings and waviness or fiber curvature [3–6].

To determine the forming behavior of fiber-reinforced composites, commonly experimental tests are carried out. The BET is widely used to characterize the in-plane shear behavior of fabrics. For the BET, the pin-jointed net assumption is often made, which assumes that fiber crossings act as rotational knots and that fibers are inextensible [7]. This assumption is generally recognized for woven fabrics and results in a deformation with three distinct shear zones [4,7,8]. For woven fabrics, in one of the zones pure shear occurs. In the second zones only half of the shear from the

<figure><img src="image_450.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under license by Materials Research Forum LLC.
first zone occurs, while in the third zones no shear deformation arises. Although NCF do not have fiber crossings and therefore the kinematic assumption of pure rotation is not valid, the BET is still used for characterizing the in-plane shear behavior. While applying the BET on bidirectional NCF it can be observed that roving slippage occurs in the transition of the shear zones defining areas with a deviation of the shear angle [9]. It has been shown experimentally that when performing the BET on UD-NCF the different zones are less distinct from each other and more zones between high shear and zero shear occur [6,10].

To simulate technically relevant parts in forming accurately, simulation models at a macroscopic level of detail are necessary. Due to better drapability, the focus in macroscopic simulation models has historically been set on woven fabrics [11–15]. Simulations of NCF have mostly been carried out on the mesoscopic level [3,16–18]. In mesoscopic simulations, the fibers and the stitching are modeled discretely and thus locally occurring forming effects can be explicitly resolved. More recently, simulation models for NCF have been developed also on a macroscopic scale [19–22]. In order to depict the relevant forming mechanisms occurring on the meso-scale, macroscopic approaches need to capture them in a homogenized way. This results in predictions of areas with a high likelihood of the occurrence of local effects. These approaches locally assume constant material properties for each element. Therefore, effects induced by their surroundings like strain gradients are not considered.

The consideration of localized forming effects in macroscopic simulations induces the need for a neighboring element method. This paper will present the implemented framework of such a method within the commercial FEM solver ABAQUS/EXPLICIT by using multiple user subroutines and allocatable arrays. It has been experimentally documented that roving slippage occurs while conducting the BET for NCF [9]. To the best of the authors’ knowledge, there is no macroscopic model that considers localized roving slippage for NCF. Since the area of occurrence is limited to areas with a large gradient of the shear angle, a method is developed to determine a spatial gradient from a scalar quantity inside a FORTRAN subroutine. The cross element gradient calculation is used to demonstrate the functionality of the framework for the neighboring method and is applied to simulation of the BET.

## Gradient calculation based on a neighboring element method

To consider the modelling of localized draping effects in macroscopic forming simulations, non-local information in the onset area of the effects needs to be known. For woven fabrics such a non-local method has already been utilized by Steer et al. [23]. They calculate the in-plane curvature of the fibers by passing a fictive Bernoulli beam through a reference element and its neighbor element by using a neighboring element method.

ABAQUS/EXPLICIT has been utilized as a simulation tool for draping simulation of various fabrics [6,12,13,20,24]. Since there is no built-in method for determining the neighboring relationships between elements, the implemented method will be presented in the following [25]. Additionally, it will also be shown how a cross element spatial gradient is calculated from a scalar quantity based on the non-local framework.

*Neighboring element method.* The implemented non-local framework combines user subroutines and regular FORTRAN subroutines. Since external subroutines are used, the user subroutine VEXTERNALDB is needed. This subroutine is called multiple times during the analysis and enables the call of external programs as well as the data transfer between external and user subroutines. While submitting an analysis within ABAQUS the so-called input file is written out. This file establishes the foundation of the neighboring method, as it contains the discretization and thereby the element-node relationship. The input file is read in and processed in external subroutines. After processing the mesh information, the element-to-element connection is determined and stored in global arrays to make them accessible within other subroutines.
first zone occurs, while in the third zones no shear deformation arises. Although NCF do not have fiber crossings and therefore the kinematic assumption of pure rotation is not valid, the BET is still used for characterizing the in-plane shear behavior. While applying the BET on bidirectional NCF it can be observed that roving slippage occurs in the transition of the shear zones defining areas with a deviation of the shear angle [9]. It has been shown experimentally that when performing the BET on UD-NCF the different zones are less distinct from each other and more zones between high shear and zero shear occur [6,10].

To simulate technically relevant parts in forming accurately, simulation models at a macroscopic level of detail are necessary. Due to better drapability, the focus in macroscopic simulation models has historically been set on woven fabrics [11–15]. Simulations of NCF have mostly been carried out on the mesoscopic level [3,16–18]. In mesoscopic simulations, the fibers and the stitching are modeled discretely and thus locally occurring forming effects can be explicitly resolved. More recently, simulation models for NCF have been developed also on a macroscopic scale [19–22]. In order to depict the relevant forming mechanisms occurring on the meso-scale, macroscopic approaches need to capture them in a homogenized way. This results in predictions of areas with a high likelihood of the occurrence of local effects. These approaches locally assume constant material properties for each element. Therefore, effects induced by their surroundings like strain gradients are not considered.

The consideration of localized forming effects in macroscopic simulations induces the need for a neighboring element method. This paper will present the implemented framework of such a method within the commercial FEM solver ABAQUS/EXPLICIT by using multiple user subroutines and allocatable arrays. It has been experimentally documented that roving slippage occurs while conducting the BET for NCF [9]. To the best of the authors’ knowledge, there is no macroscopic model that considers localized roving slippage for NCF. Since the area of occurrence is limited to areas with a large gradient of the shear angle, a method is developed to determine a spatial gradient from a scalar quantity inside a FORTRAN subroutine. The cross element gradient calculation is used to demonstrate the functionality of the framework for the neighboring method and is applied to simulation of the BET.

## Gradient calculation based on a neighboring element method

To consider the modelling of localized draping effects in macroscopic forming simulations, non-local information in the onset area of the effects needs to be known. For woven fabrics such a non-local method has already been utilized by Steer et al. [23]. They calculate the in-plane curvature of the fibers by passing a fictive Bernoulli beam through a reference element and its neighbor element by using a neighboring element method.

ABAQUS/EXPLICIT has been utilized as a simulation tool for draping simulation of various fabrics [6,12,13,20,24]. Since there is no built-in method for determining the neighboring relationships between elements, the implemented method will be presented in the following [25]. Additionally, it will also be shown how a cross element spatial gradient is calculated from a scalar quantity based on the non-local framework.

*Neighboring element method.* The implemented non-local framework combines user subroutines and regular FORTRAN subroutines. Since external subroutines are used, the user subroutine VEXTERNALDB is needed. This subroutine is called multiple times during the analysis and enables the call of external programs as well as the data transfer between external and user subroutines. While submitting an analysis within ABAQUS the so-called input file is written out. This file establishes the foundation of the neighboring method, as it contains the discretization and thereby the element-node relationship. The input file is read in and processed in external subroutines. After processing the mesh information, the element-to-element connection is determined and stored in global arrays to make them accessible within other subroutines.
Gradient calculation. The construction of a cross element spatial gradient from an element-based scalar field variable is commonly applied for solving computational fluid dynamic problems. An extensively investigated and commonly used method is the cell-average weighted least-squares (WLSQ) gradient method [26,27]. The WLSQ gradient method will be derived below and the adaptions for the draping simulation will be shown. The derivation is performed for an arbitrary number of neighbor elements while the implementation will only consider the stencil of first edge neighbors (en1) as shown in Figure 1.

Starting point for the WLSQ gradient method yields the first-order Taylor series expansion of the shear angle $\gamma$

$$ \gamma_J = \gamma_I + \nabla\gamma_I \cdot d_{IJ} + H \quad (1) $$

with $\gamma_I$ and $\gamma_J$ being the shear angle on the respective element I, J and $\nabla\gamma_I$ being the spatial gradient of the element I in x-, y- and z-direction. The higher-order terms of the Taylor series H are neglected in the following. The variable $d_{IJ}$ is the direction vector from centroid of element I to centroid of element J.

<figure><img src="image_440.png" /><figcaption>Fig. 1: Stencil of first edge neighbors (en1) showing the reference element I (grey) and its neighboring elements J.</figcaption></figure>

It is worth noting that the indexes in capital letters are not indexes in the context of Einstein's summation convention but serve to identify the reference element I and the neighboring elements J. For an arbitrary number of elements N this leads to an overdetermined linear system of the form

$$ \begin{bmatrix} d_{I1}^T \\ d_{I2}^T \\ \vdots \\ d_{IN}^T \end{bmatrix} \begin{bmatrix} \partial\gamma/\partial x \\ \partial\gamma/\partial y \\ \partial\gamma/\partial z \end{bmatrix}_I = \begin{bmatrix} \Delta x_{I1} & \Delta y_{I1} & \Delta z_{I1} \\ \Delta x_{I2} & \Delta y_{I2} & \Delta z_{I2} \\ \vdots & \vdots & \vdots \\ \Delta x_{IN} & \Delta y_{IN} & \Delta z_{IN} \end{bmatrix} \begin{bmatrix} \partial\gamma/\partial x \\ \partial\gamma/\partial y \\ \partial\gamma/\partial z \end{bmatrix}_I = \begin{bmatrix} \Delta\gamma_{I1} \\ \Delta\gamma_{I2} \\ \vdots \\ \Delta\gamma_{IN} \end{bmatrix} \quad (2) $$

which can be summarized as $D\nabla\gamma_I = b$. The operator $\Delta$ can be interpreted as $\Delta(\cdot)_{IJ} = (\cdot)_J - (\cdot)_I$. As large element aspect ratios can be attained during the draping simulation, stabilization in the form of weights is introduced to Eq. (2), resulting in

$$ W D \nabla \gamma_I = W b. \quad (3) $$

The weighting matrix $W$ is a diagonal matrix containing the weights $w_{iIJ} = 1/\|d_{IJ}\|^p$ with $i = 1, ..., N$ and the parameter $p$. Solving the system of Eq. (3) is done by multiplying both sides from the left with $D^T W$. This results in the left-hand side having a quadratic form and therefore it can be inverted. This leads to
Gradient calculation. The construction of a cross element spatial gradient from an element-based scalar field variable is commonly applied for solving computational fluid dynamic problems. An extensively investigated and commonly used method is the cell-average weighted least-squares (WLSQ) gradient method [26,27]. The WLSQ gradient method will be derived below and the adaptions for the draping simulation will be shown. The derivation is performed for an arbitrary number of neighbor elements while the implementation will only consider the stencil of first edge neighbors (en1) as shown in Figure 1.

Starting point for the WLSQ gradient method yields the first-order Taylor series expansion of the shear angle $\gamma$

$$ \gamma_J = \gamma_I + \nabla\gamma_I \cdot d_{IJ} + H \quad (1) $$

with $\gamma_I$ and $\gamma_J$ being the shear angle on the respective element $I, J$ and $\nabla\gamma_I$ being the spatial gradient of the element $I$ in x-, y- and z-direction. The higher-order terms of the Taylor series $H$ are neglected in the following. The variable $d_{IJ}$ is the direction vector from centroid of element $I$ to centroid of element $J$.

<figure><img src="image_451.png" /><figcaption>Fig. 1: Stencil of first edge neighbors (en1) showing the reference element $I$ (grey) and its neighboring elements $J$.</figcaption></figure>

It is worth noting that the indexes in capital letters are not indexes in the context of Einstein's summation convention but serve to identify the reference element $I$ and the neighboring elements $J$. For an arbitrary number of elements $N$ this leads to an overdetermined linear system of the form

$$ \begin{bmatrix} d_{I1}^T \\ d_{I2}^T \\ \vdots \\ d_{IN}^T \end{bmatrix} \begin{bmatrix} \partial\gamma/\partial x \\ \partial\gamma/\partial y \\ \partial\gamma/\partial z \end{bmatrix}_I = \begin{bmatrix} \Delta x_{I1} & \Delta y_{I1} & \Delta z_{I1} \\ \Delta x_{I2} & \Delta y_{I2} & \Delta z_{I2} \\ \vdots & \vdots & \vdots \\ \Delta x_{IN} & \Delta y_{IN} & \Delta z_{IN} \end{bmatrix} \begin{bmatrix} \partial\gamma/\partial x \\ \partial\gamma/\partial y \\ \partial\gamma/\partial z \end{bmatrix}_I = \begin{bmatrix} \Delta\gamma_{I1} \\ \Delta\gamma_{I2} \\ \vdots \\ \Delta\gamma_{IN} \end{bmatrix} \quad (2) $$

which can be summarized as $D\nabla\gamma_I = b$. The operator $\Delta$ can be interpreted as $\Delta(\cdot)_{IJ} = (\cdot)_J - (\cdot)_I$. As large element aspect ratios can be attained during the draping simulation, stabilization in the form of weights is introduced to Eq. (2), resulting in

$$ W D \nabla \gamma_I = W b. \quad (3) $$

The weighting matrix $W$ is a diagonal matrix containing the weights $w_{iIJ} = 1/\|d_{IJ}\|^p$ with $i = 1, ..., N$ and the parameter $p$. Solving the system of Eq. (3) is done by multiplying both sides from the left with $D^T W$. This results in the left-hand side having a quadratic form and therefore it can be inverted. This leads to
$$ \nabla\gamma_I = (D^TWWD)^{-1}D^TWWb. \qquad (4) $$

From Eq. (4) it can be seen that through the matrix $D$ and its transpose the distance between the centroids is squared in the first inverted part while it is only linear in the second part. This yields a dependency of the dimension of the gradient on the element size of the mesh. To circumvent this, the matrix $D$ is normed as follows

$$ D_{\text{norm}} = \begin{bmatrix} d_{I1}^T / \|d_{I1}\| \\ d_{I2}^T / \|d_{I2}\| \\ \vdots \\ d_{IN}^T / \|d_{IN}\| \end{bmatrix} = \begin{bmatrix} d_{\text{norm},I1}^T \\ d_{\text{norm},I2}^T \\ \vdots \\ d_{\text{norm},IN}^T \end{bmatrix}. \qquad (5) $$

The vector $d_{\text{norm},Ij}$ and the matrix $D_{\text{norm}}$ therefore only contain directional information. The matrix $D_{\text{norm}}$ is substituted into Eq. (4) instead of $D$.

Calculating the gradient using the WLSQ method leads to the gradient components expressed in the global configuration {$e_i$}. To utilize the gradient in the context of forming simulation its components need to be expressed in the current configuration {$a_i$}. Figure 2 exemplifies the transformation between the different configurations, resulting in the non-orthogonal principal material orientations {$a_i$} in the current configuration. Since the initial fiber orientation {$A_i$} does not have to be aligned with the global configuration, a preceding transformation is necessary. Based on the assumption of a cross-ply NCF this transformation is a simple orthogonal rotation $R$. The deformation gradient $F$ relates the initial {$A_i$} and the current configuration {$a_i$}. Therefore, the gradient components with respect to the current configuration are obtained through

$$ \nabla\gamma_{i,\{a_i\}} = (FR)^{-1}\nabla\gamma_{j,\{e_j\}}. \qquad (6) $$

<figure><img src="image_441.png" /><figcaption>Fig. 2: Schematic illustration of the transformations between global, initial, and current configuration.</figcaption></figure>

## Numerical verification

To investigate the WLSQ gradient method for forming simulation of NCF a planar small-scale example is examined. The two applied element-based fields of a synthetic scalar variable are shown in Figure 3. For investigational purposes, in both depicted cases the initial and current configurations coincide. The applied field in Figure 3 a) will be referred to as diagonal distribution. The rotation matrix $R_{\text{diag}}$ corresponds to a rotation of 45° about the $e_3$-axis. The scalar field shown in Figure 3 b) will be referred to as horizontal distribution. All configurations coincide and it follows that $R_{\text{horiz}} = I$.
$$ \nabla \gamma_I = (D^T WWD)^{-1} D^T WWb. \qquad (4) $$

From Eq. (4) it can be seen that through the matrix $D$ and its transpose the distance between the centroids is squared in the first inverted part while it is only linear in the second part. This yields a dependency of the dimension of the gradient on the element size of the mesh. To circumvent this, the matrix $D$ is normed as follows

$$ D_{\text{norm}} = \begin{bmatrix} d_{I1}^T / \|d_{I1}\| \\ d_{I2}^T / \|d_{I2}\| \\ \vdots \\ d_{IN}^T / \|d_{IN}\| \end{bmatrix} = \begin{bmatrix} d_{\text{norm},I1}^T \\ d_{\text{norm},I2}^T \\ \vdots \\ d_{\text{norm},IN}^T \end{bmatrix}. \qquad (5) $$

The vector $d_{\text{norm},Ij}$ and the matrix $D_{\text{norm}}$ therefore only contain directional information. The matrix $D_{\text{norm}}$ is substituted into Eq. (4) instead of $D$.

Calculating the gradient using the WLSQ method leads to the gradient components expressed in the global configuration {$e_i$}. To utilize the gradient in the context of forming simulation its components need to be expressed in the current configuration {$a_i$}. Figure 2 exemplifies the transformation between the different configurations, resulting in the non-orthogonal principal material orientations {$a_i$} in the current configuration. Since the initial fiber orientation {$A_i$} does not have to be aligned with the global configuration, a preceding transformation is necessary. Based on the assumption of a cross-ply NCF this transformation is a simple orthogonal rotation $R$. The deformation gradient $F$ relates the initial {$A_i$} and the current configuration {$a_i$}. Therefore, the gradient components with respect to the current configuration are obtained through

$$ \nabla \gamma_{i,\{a_i\}} = (FR)^{-1} \nabla \gamma_{j,\{e_j\}}. \qquad (6) $$

<figure><img src="image_452.png" /><figcaption>Fig. 2: Schematic illustration of the transformations between global, initial, and current configuration.</figcaption></figure>

## Numerical verification

To investigate the WLSQ gradient method for forming simulation of NCF a planar small-scale example is examined. The two applied element-based fields of a synthetic scalar variable are shown in Figure 3. For investigational purposes, in both depicted cases the initial and current configurations coincide. The applied field in Figure 3 a) will be referred to as diagonal distribution. The rotation matrix $R_{\text{diag}}$ corresponds to a rotation of 45° about the $e_3$-axis. The scalar field shown in Figure 3 b) will be referred to as horizontal distribution. All configurations coincide and it follows that $R_{\text{horiz}} = I$.
<figure><img src="image_442.png" /><figcaption>Fig. 3: Visualization of the two considered distributions of a synthetic element-based scalar variable $\gamma$.</figcaption></figure>

The calculated gradient of the synthetic scalar variable is shown in Figure 4. The elements located at the outer edge of the meshes are not fully surrounded by neighbor elements. Whenever a reference element is missing one neighboring element and if the remaining neighboring element has a different value of the scalar field variable than the reference element ($\Delta\gamma_{IJ} \neq 0$), this leads to so-called boundary effects. These boundary effects are recognizable in the bottom left and top right corner in Figure 4. Due to the undeformed state of the mesh and the fact that triangle elements are used, the missing element leads to one of the two components being zero (cf. Fig. 4 a) and b)). This depends on whether the vector $d_{norm,IJ}$ of the missing element $J$ contains the $e_1$- or the $e_2$-direction. If $d_{norm,IJ} = |e_1|$ does not appear in the matrix $D_{norm}$ it results in the gradient component in $e_2$-direction being zero. The non-zero components then have a value that exhibits the value of the scalar field variable $\gamma$. The gradient components of the inner elements in Figures 4 a) and b) in the global configuration are equal in amount. Rotating the gradient to the current configuration using $R_{diag}$ results in disappearing of the gradient component in the direction of $a_1$ for the inner elements (cf. Figure 4 c)). Only the gradient in the $a_2$-direction exists in case of neglecting the elements with a boundary influence. The gradient component's magnitude of the inner elements is less than that of the element-based scalar variable. However, the value of the gradient is also dependent on the choice of the weights. In this work a weight with $p = 2$ is used. Comparing the gradient components with the given diagonal distribution, this is a sufficient result.

Applying the WLSQ gradient method on the horizontal distribution leads to the result depicted in Figure 5. The Figure shows that a gradient occurs in the $e_1$-direction. Comparing this to the given horizontal distribution, the existence of the gradient component in the direction of $e_1$ is not desirable. However, the occurrence of a gradient in the $e_1$-direction cannot be prevented. The reason for this is the WLSQ method, as the direction between the reference and neighboring element $d_{norm,IJ}$ in the matrix $D_{norm}$ also influences the gradient although the difference between the scalar variable $\Delta\gamma_{IJ}$ vanishes. However, this effect is a local error that disappears in a global consideration when integrating the incorrect component over the domain. It can be stated that the gradient in the $e_2$-direction is 1.5 times greater in terms of the magnitude.
<figure><img src="image_453.png" /><figcaption>Fig. 3: Visualization of the two considered distributions of a synthetic element-based scalar variable $\gamma$.</figcaption></figure>

The calculated gradient of the synthetic scalar variable is shown in Figure 4. The elements located at the outer edge of the meshes are not fully surrounded by neighbor elements. Whenever a reference element is missing one neighboring element and if the remaining neighboring element has a different value of the scalar field variable than the reference element ($\Delta\gamma_{IJ} \neq 0$), this leads to so-called boundary effects. These boundary effects are recognizable in the bottom left and top right corner in Figure 4. Due to the undeformed state of the mesh and the fact that triangle elements are used, the missing element leads to one of the two components being zero (cf. Fig. 4 a) and b). This depends on whether the vector $d_{norm,IJ}$ of the missing element $J$ contains the $e_1$- or the $e_2$-direction. If $d_{norm,IJ} = |e_1|$ does not appear in the matrix $D_{norm}$ it results in the gradient component in $e_2$-direction being zero. The non-zero components then have a value that exhibits the value of the scalar field variable $\gamma$. The gradient components of the inner elements in Figures 4 a) and b) in the global configuration are equal in amount. Rotating the gradient to the current configuration using $R_{diag}$ results in disappearing of the gradient component in the direction of $a_1$ for the inner elements (cf. Figure 4 c)). Only the gradient in the $a_2$-direction exists in case of neglecting the elements with a boundary influence. The gradient component's magnitude of the inner elements is less than that of the element-based scalar variable. However, the value of the gradient is also dependent on the choice of the weights. In this work a weight with $p = 2$ is used. Comparing the gradient components with the given diagonal distribution, this is a sufficient result.

Applying the WLSQ gradient method on the horizontal distribution leads to the result depicted in Figure 5. The Figure shows that a gradient occurs in the $e_1$-direction. Comparing this to the given horizontal distribution, the existence of the gradient component in the direction of $e_1$ is not desirable. However, the occurrence of a gradient in the $e_1$-direction cannot be prevented. The reason for this is the WLSQ method, as the direction between the reference and neighboring element $d_{norm,IJ}$ in the matrix $D_{norm}$ also influences the gradient although the difference between the scalar variable $\Delta\gamma_{IJ}$ vanishes. However, this effect is a local error that disappears in a global consideration when integrating the incorrect component over the domain. It can be stated that the gradient in the $e_2$-direction is 1.5 times greater in terms of the magnitude.
<figure><img src="image_443.png" /><figcaption>Fig. 4: Results of the gradient component 1 a) and component 2 b) in the global configuration {$e_i$} and component 1 c) and component 2 d) in the current configuration {$a_i$} of the diagonal distribution.</figcaption></figure>

<figure><img src="image_444.png" /><figcaption>Fig. 5: Results of the gradient component 1 a) and component 2 b) in the global configuration {$e_i$} of the horizontal distribution. The current configuration {$a_i$} equals the global configuration {$e_i$} in this case.</figcaption></figure>

As shown in the above small-scale examples, the WLSQ gradient method as presented here is not yet able to make a quantitative, but only a qualitative statement on the occurrence of a gradient of an element-based scalar variable. It is sufficient to localize the transition zones between areas with a high deviation in a scalar variable such as the shear angle. The potential of the method can be seen, but challenges still need to be overcome. In Figure 4 the problem of elements at the
<figure><img src="image_454.png" /><figcaption>Fig. 4: Results of the gradient component 1 a) and component 2 b) in the global configuration {$e_i$} and component 1 c) and component 2 d) in the current configuration {$a_i$} of the diagonal distribution.</figcaption></figure>

<figure><img src="image_455.png" /><figcaption>Fig. 5: Results of the gradient component 1 a) and component 2 b) in the global configuration {$e_i$} of the horizontal distribution. The current configuration {$a_i$} equals the global configuration {$e_i$} in this case.</figcaption></figure>

As shown in the above small-scale examples, the WLSQ gradient method as presented here is not yet able to make a quantitative, but only a qualitative statement on the occurrence of a gradient of an element-based scalar variable. It is sufficient to localize the transition zones between areas with a high deviation in a scalar variable such as the shear angle. The potential of the method can be seen, but challenges still need to be overcome. In Figure 4 the problem of elements at the
boundary is pointed out, while an incorrect component is predicted at local scale, which disappears in a global consideration (cf. Fig. 5).

## Application on the bias extension test

The model for simulating the BET is created in ABAQUS/EXPLICIT. The dimensions of the considered UD-NCF are 320 mm × 160 mm. The geometry is discretized using 4096 triangular elements, each with an edge length of 5 mm. The nodes at $x_2 = 0$ mm are pinned in their degrees of freedom. The simulation is carried out displacement-controlled, with the displacement of $u_2 = 80$ mm being applied to the nodes at $x_2 = 320$ mm. The fibers and stitching are initially orthogonal, which is reflected by the {$A_i$}-frame. The fibers are oriented along the $A_1$-direction and the stitching along the $A_2$-direction. To transfer the {$A_i$}-frame to the global frame {$e_i$} the rotation matrix $R_{diag}$ is used.

<figure><img src="image_445.png" /><figcaption>Fig. 6: Simulation results of the BET showing the shear angle distribution in a) and the components of the gradients in the global b), c) and current d), e) configuration.</figcaption></figure>

The material model chosen is the hyperelastic St. Venant-Kirchhoff model [28]. The model is implemented in the user subroutine VUMAT. The formulation for the second Piola-Kirchhoff stress tensor $S$ reads

$$S = C[E], \qquad (7)$$

where $C$ is the stiffness tensor and $E$ is the Green-Lagrangian strain tensor. The stiffnesses remain constant during the simulation and their values are selected to qualitatively represent the forming behavior of the NCF for simulating the BET. Since the fibers are significantly stiffer than the stitching, an order of magnitude of three is defined between their stiffnesses. Consequently, for the Young's modulus of the fiber follows $E_1 = 1000$ MPa and for Young's modulus of the stitching $E_2 = 1$ MPa. For the shear stiffness a shear modulus $G_{12} = 0.5$ MPa is used. If the shear modulus is set too low, the middle section of the NCF will collapse. The Poisson's ratio is omitted. Since
boundary is pointed out, while an incorrect component is predicted at local scale, which disappears in a global consideration (cf. Fig. 5).

## Application on the bias extension test

The model for simulating the BET is created in ABAQUS/EXPLICIT. The dimensions of the considered UD-NCF are 320 mm × 160 mm. The geometry is discretized using 4096 triangular elements, each with an edge length of 5 mm. The nodes at $x_2 = 0$ mm are pinned in their degrees of freedom. The simulation is carried out displacement-controlled, with the displacement of $u_2 = 80$ mm being applied to the nodes at $x_2 = 320$ mm. The fibers and stitching are initially orthogonal, which is reflected by the {$A_i$}-frame. The fibers are oriented along the $A_1$-direction and the stitching along the $A_2$-direction. To transfer the {$A_i$}-frame to the global frame {$e_i$} the rotation matrix $R_{diag}$ is used.

<figure><img src="image_456.png" /><figcaption>Fig. 6: Simulation results of the BET showing the shear angle distribution in a) and the components of the gradients in the global b), c) and current d), e) configuration.</figcaption></figure>

The material model chosen is the hyperelastic St. Venant-Kirchhoff model [28]. The model is implemented in the user subroutine VUMAT. The formulation for the second Piola-Kirchhoff stress tensor $S$ reads

$$S = C[E], \qquad (7)$$

where $C$ is the stiffness tensor and $E$ is the Green-Lagrangian strain tensor. The stiffnesses remain constant during the simulation and their values are selected to qualitatively represent the forming behavior of the NCF for simulating the BET. Since the fibers are significantly stiffer than the stitching, an order of magnitude of three is defined between their stiffnesses. Consequently, for the Young's modulus of the fiber follows $E_1 = 1000$ MPa and for Young's modulus of the stitching $E_2 = 1$ MPa. For the shear stiffness a shear modulus $G_{12} = 0.5$ MPa is used. If the shear modulus is set too low, the middle section of the NCF will collapse. The Poisson's ratio is omitted. Since
the BET is a planar problem these material parameters are sufficient to describe and construct the orthotropic stiffness tensor $C$ in Eq. (7).

In Figure 6 the simulation results of the BET are visualized. The outermost row of elements is neglected in the figure because the gradient has a significantly higher value there than in the interior of the mesh due to the boundary effects. Since the fiber and stitching stiffnesses are three orders of magnitude apart, the three zones with constant shear deformation, as characteristic for woven fabrics, do not occur [7]. The deformation in the middle zone is not a pure shear deformation for UD-NCF with a maximal shear angle $γ ≈ 49°$ (cf. Figure 6 a)). This behavior was demonstrated experimentally on the BET [6,10]. Along the connection of $\overline{AB}$ as well as $\overline{CD}$ there is a steep decrease in the shear angle. In addition, a smaller deviation can be seen in the vicinity of E and F respectively. Applying the introduced WLSQ method leads to the gradient components shown in Figure 6 b)-e). The gradient components along the transition zones $\overline{AB}$ and $\overline{CD}$ stand out particularly. In the remaining mesh except the areas around E and F, the gradient components show very small values $∇γ_i < 1.0°$. The component in the $e_2$-direction is more pronounced than in the $e_1$-direction. The rotation and application of the deformation gradient leads to the components in the material configuration {$a_i$} in Figure 6 d) and e). Similarly, it is observable that the component along the orientation of the fibers $a_1$ is smaller than in the direction of the stitching $a_2$. The qualitative agreement with the shear angle distribution in Figure 6 a) is thus given.

## Conclusion and outlook

A general framework in ABAQUS/EXPLICIT which allows to consider localized draping effects in macroscopic forming simulations is presented. For this purpose, a combination of external FORTRAN subroutines and user subroutines is used to implement a neighboring element method. The neighbor information stored in global arrays can be retrieved by any other subroutine at runtime. The functionality of the neighboring element method is demonstrated in the calculation of a cross element gradient of a scalar element-based variable utilizing the cell-average weighted least-squares (WLSQ) gradient method. The WLSQ gradient method is first examined on two small-scale examples. It is shown that a qualitative statement about the occurrence of gradients can be made. Applied on macroscopic simulation of the bias-extension test (BET), the transition zones occurring at deviating shear angles are detected.

To improve of the WLSQ gradient method further extensions are necessary. These include the consideration of additional stencils, the mitigation of the boundary effects, and addressing the local error of the gradient components. Subsequently the gradient of the shear angle computed using the WLSQ method is to be integrated into a hyperelastic material approach based on pseudo-invariants by using a suitable coupling mechanism in order to model the localized fiber roving slippage in the transition zones.

## Acknowledgement

The authors would like to thank the German Federal Ministry for Economic Affairs and Climate Action (BMWK) for funding the LuFo research project Electra (20W1912D), for which the work has been carried out. Furthermore, the authors would like to thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) and the French National Research Agency (ANR) for funding the collaborative project "Composite forming simulation for non-crimp fabrics based on generalized continuum approaches" (DFG: 431354059, ANR: ANR-19-CE06-0031-012), in which essential preliminary work had been conducted. The work is also part of the Heisenberg project "Digitalization of fiber-reinforced polymer processes for resource-efficient manufacturing of lightweight components", funded by the DFG (798455807141).
the BET is a planar problem these material parameters are sufficient to describe and construct the orthotropic stiffness tensor $C$ in Eq. (7).

In Figure 6 the simulation results of the BET are visualized. The outermost row of elements is neglected in the figure because the gradient has a significantly higher value there than in the interior of the mesh due to the boundary effects. Since the fiber and stitching stiffnesses are three orders of magnitude apart, the three zones with constant shear deformation, as characteristic for woven fabrics, do not occur [7]. The deformation in the middle zone is not a pure shear deformation for UD-NCF with a maximal shear angle $γ ≈ 49°$ (cf. Figure 6 a)). This behavior was demonstrated experimentally on the BET [6,10]. Along the connection of $\overline{AB}$ as well as $\overline{CD}$ there is a steep decrease in the shear angle. In addition, a smaller deviation can be seen in the vicinity of E and F respectively. Applying the introduced WLSQ method leads to the gradient components shown in Figure 6 b)-e). The gradient components along the transition zones $\overline{AB}$ and $\overline{CD}$ stand out particularly. In the remaining mesh except the areas around E and F, the gradient components show very small values $∇γ_i < 1.0°$. The component in the $e_2$-direction is more pronounced than in the $e_1$-direction. The rotation and application of the deformation gradient leads to the components in the material configuration {$a_i$} in Figure 6 d) and e). Similarly, it is observable that the component along the orientation of the fibers $a_1$ is smaller than in the direction of the stitching $a_2$. The qualitative agreement with the shear angle distribution in Figure 6 a) is thus given.

## Conclusion and outlook

A general framework in ABAQUS/EXPLICIT which allows to consider localized draping effects in macroscopic forming simulations is presented. For this purpose, a combination of external FORTRAN subroutines and user subroutines is used to implement a neighboring element method. The neighbor information stored in global arrays can be retrieved by any other subroutine at runtime. The functionality of the neighboring element method is demonstrated in the calculation of a cross element gradient of a scalar element-based variable utilizing the cell-average weighted least-squares (WLSQ) gradient method. The WLSQ gradient method is first examined on two small-scale examples. It is shown that a qualitative statement about the occurrence of gradients can be made. Applied on macroscopic simulation of the bias-extension test (BET), the transition zones occurring at deviating shear angles are detected.

To improve of the WLSQ gradient method further extensions are necessary. These include the consideration of additional stencils, the mitigation of the boundary effects, and addressing the local error of the gradient components. Subsequently the gradient of the shear angle computed using the WLSQ method is to be integrated into a hyperelastic material approach based on pseudo-invariants by using a suitable coupling mechanism in order to model the localized fiber roving slippage in the transition zones.

## Acknowledgement

The authors would like to thank the German Federal Ministry for Economic Affairs and Climate Action (BMWK) for funding the LuFo research project Electra (20W1912D), for which the work has been carried out. Furthermore, the authors would like to thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) and the French National Research Agency (ANR) for funding the collaborative project "Composite forming simulation for non-crimp fabrics based on generalized continuum approaches" (DFG: 431354059, ANR: ANR-19-CE06-0031-012), in which essential preliminary work had been conducted. The work is also part of the Heisenberg project "Digitalization of fiber-reinforced polymer processes for resource-efficient manufacturing of lightweight components", funded by the DFG (798455807141).
## References

[1] F. Henning, L. Kärger, D. Dörr, F.J. Schirmaier, J. Seuffert, A. Bernath, Fast processing and continuous simulation of automotive structural composite components, Composites Science and Technology 171 (2019) 261–279. https://doi.org/10.1016/j.compscitech.2018.12.007

[2] C. Krogh, J.A. Kepler, J. Jakobsen, Pure and simple: investigating the in-plane shear kinematics of a quasi-unidirectional glass fiber non-crimp fabric using the bias-extension test, Int J Mater Form 19 (2021) 1483–1495. https://doi.org/10.1007/s12289-021-01642-8

[3] P. Böhler, F. Härtel, P. Middendorf, Identification of Forming Limits for Unidirectional Carbon Textiles in Reality and Mesoscopic Simulation, Key Engineering Materials 554-557 (2013) 423–432. https://doi.org/10.4028/www.scientific.net/KEM.554-557.423

[4] N. Hamila, P. Boisse, Locking in simulation of composite reinforcement deformations. Analysis and treatment, Composites Part A: Applied Science and Manufacturing 53 (2013) 109–117. https://doi.org/10.1016/j.compositesa.2013.06.001

[5] E. Kunze, S. Galkin, R. Böhm, M. Gude, L. Kärger, The Impact of Draping Effects on the Stiffness and Failure Behavior of Unidirectional Non-Crimp Fabric Fiber Reinforced Composites, Materials 13 (2020). https://doi.org/10.3390/ma13132959

[6] F.J. Schirmaier, K.A. Weidenmann, L. Kärger, F. Henning, Characterisation of the draping behaviour of unidirectional non-crimp fabrics (UD-NCF), Composites Part A: Applied Science and Manufacturing 80 (2016) 28–38. https://doi.org/10.1016/j.compositesa.2015.10.004

[7] G. Lebrun, M.N. Bureau, J. Denault, Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics, Selected Papers from the Symposium on Design and Manufacturing of Composites 61 (2003) 341–352. https://doi.org/10.1016/S0263-8223(03)00057-6

[8] J. Launay, G. Hivet, A.V. Duong, P. Boisse, Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements, Composites Science and Technology 68 (2008) 506–515. https://doi.org/10.1016/j.compscitech.2007.06.021

[9] S. Bel, P. Boisse, F. Dumont, Analyses of the Deformation Mechanisms of Non-Crimp Fabric Composite Reinforcements during Preforming, Applied Composite Materials 19 (2012) 513–528. https://doi.org/10.1007/s10443-011-9207-x

[10] B. Schäfer, R. Zheng, N. Naouar, L. Kärger, Membrane behavior of uni- and bidirectional non-crimp fabrics in off-axis-tension tests, Int J Mater Form 16 (2023) 1–15. https://doi.org/10.1007/s12289-023-01792-x

[11] M.A. Khan, T. Mabrouki, E. Vidal-Sallé, P. Boisse, Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark, Journal of Materials Processing Technology 210 (2010) 378–388. https://doi.org/10.1016/j.jmatprotec.2009.09.027

[12] M. Machado, M. Fischlschweiger, Z. Major, A rate-dependent non-orthogonal constitutive model for describing shear behaviour of woven reinforced thermoplastic composites, Composites Part A: Applied Science and Manufacturing 80 (2016) 194–203. https://doi.org/10.1016/j.compositesa.2015.10.028

[13] F. Schäfer, H.O. Werner, F. Henning, L. Kärger, A hyperelastic material model considering biaxial coupling of tension–compression and shear for the forming simulation of woven fabrics, Composites Part A: Applied Science and Manufacturing 165 (2023). https://doi.org/10.1016/j.compositesa.2022.107323
## References

[1] F. Henning, L. Kärger, D. Dörr, F.J. Schirmaier, J. Seuffert, A. Bernath, Fast processing and continuous simulation of automotive structural composite components, Composites Science and Technology 171 (2019) 261–279. https://doi.org/10.1016/j.compscitech.2018.12.007

[2] C. Krogh, J.A. Kepler, J. Jakobsen, Pure and simple: investigating the in-plane shear kinematics of a quasi-unidirectional glass fiber non-crimp fabric using the bias-extension test, Int J Mater Form 19 (2021) 1483–1495. https://doi.org/10.1007/s12289-021-01642-8

[3] P. Böhler, F. Härtel, P. Middendorf, Identification of Forming Limits for Unidirectional Carbon Textiles in Reality and Mesoscopic Simulation, Key Engineering Materials 554-557 (2013) 423–432. https://doi.org/10.4028/www.scientific.net/KEM.554-557.423

[4] N. Hamila, P. Boisse, Locking in simulation of composite reinforcement deformations. Analysis and treatment, Composites Part A: Applied Science and Manufacturing 53 (2013) 109–117. https://doi.org/10.1016/j.compositesa.2013.06.001

[5] E. Kunze, S. Galkin, R. Böhm, M. Gude, L. Kärger, The Impact of Draping Effects on the Stiffness and Failure Behavior of Unidirectional Non-Crimp Fabric Fiber Reinforced Composites, Materials 13 (2020). https://doi.org/10.3390/ma13132959

[6] F.J. Schirmaier, K.A. Weidenmann, L. Kärger, F. Henning, Characterisation of the draping behaviour of unidirectional non-crimp fabrics (UD-NCF), Composites Part A: Applied Science and Manufacturing 80 (2016) 28–38. https://doi.org/10.1016/j.compositesa.2015.10.004

[7] G. Lebrun, M.N. Bureau, J. Denault, Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics, Selected Papers from the Symposium on Design and Manufacturing of Composites 61 (2003) 341–352. https://doi.org/10.1016/S0263-8223(03)00057-6

[8] J. Launay, G. Hivet, A.V. Duong, P. Boisse, Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements, Composites Science and Technology 68 (2008) 506–515. https://doi.org/10.1016/j.compscitech.2007.06.021

[9] S. Bel, P. Boisse, F. Dumont, Analyses of the Deformation Mechanisms of Non-Crimp Fabric Composite Reinforcements during Preforming, Applied Composite Materials 19 (2012) 513–528. https://doi.org/10.1007/s10443-011-9207-x

[10] B. Schäfer, R. Zheng, N. Naouar, L. Kärger, Membrane behavior of uni- and bidirectional non-crimp fabrics in off-axis-tension tests, Int J Mater Form 16 (2023) 1–15. https://doi.org/10.1007/s12289-023-01792-x

[11] M.A. Khan, T. Mabrouki, E. Vidal-Sallé, P. Boisse, Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark, Journal of Materials Processing Technology 210 (2010) 378–388. https://doi.org/10.1016/j.jmatprotec.2009.09.027

[12] M. Machado, M. Fischlschweiger, Z. Major, A rate-dependent non-orthogonal constitutive model for describing shear behaviour of woven reinforced thermoplastic composites, Composites Part A: Applied Science and Manufacturing 80 (2016) 194–203. https://doi.org/10.1016/j.compositesa.2015.10.028

[13] F. Schäfer, H.O. Werner, F. Henning, L. Kärger, A hyperelastic material model considering biaxial coupling of tension–compression and shear for the forming simulation of woven fabrics, Composites Part A: Applied Science and Manufacturing 165 (2023). https://doi.org/10.1016/j.compositesa.2022.107323
[14] Y. Gong, D. Yan, Y. Yao, R. Wei, H. Hu, P. Xu, X. Peng, An Anisotropic Hyperelastic Constitutive Model with Tension-Shear Coupling for Woven Composite Reinforcements, Int. J. Appl. Mechanics 09 (2017). https://doi.org/10.1142/S1758825117500831

[15] Y. Yao, X. Huang, X. Peng, P. Liu, G. Youkun, An anisotropic hyperelastic constitutive model for plain weave fabric considering biaxial tension coupling, Textile Research Journal 89 (2019) 434-444. https://doi.org/10.1177/0040517517748495

[16] S. Lomov, D. Ivanov, I. Verpoest, M. Zako, T. Kurashiki, H. Nakai, S. Hirosawa, Meso-FE modelling of textile composites: Road map, data flow and algorithms, Composites Science and Technology 67 (2007) 1870-1891. https://doi.org/10.1016/j.compscitech.2006.10.017

[17] J. Sirtautas, A.K. Pickett, P. Lépicier, A mesoscopic model for coupled drape-infusion simulation of biaxial Non-Crimp Fabric, Composites Part B: Engineering 47 (2013) 48-57. https://doi.org/10.1016/j.compositesb.2012.09.088

[18] G. Creech, A. Pickett, Meso-modelling of Non-Crimp Fabric composites for coupled drape and failure analysis, Journal of Materials Science 41 (2006) 6725-6736. https://doi.org/10.1007/s10853-006-0213-6

[19] F.J. Schirmaier, D. Dörr, F. Henning, L. Kärger, A macroscopic approach to simulate the forming behaviour of stitched unidirectional non-crimp fabrics (UD-NCF), Composites Part A: Applied Science and Manufacturing 102 (2017) 322-335. https://doi.org/10.1016/j.compositesa.2017.08.009

[20] B. Schäfer, S. Haas, P. Boisse, L. Kärger, Investigation of the Membrane Behavior of UD-NCF in Macroscopic Forming Simulations, Key Engineering Materials 926 (2022) 1413-1422. https://doi.org/10.4028/p-2977b4

[21] V.N. Khiêm, H. Krieger, M. Itskov, T. Gries, S.E. Stapleton, An averaging based hyperelastic modeling and experimental analysis of non-crimp fabrics, International Journal of Solids and Structures. https://doi.org/10.1016/j.ijsolstr.2016.12.018

[22] B. Schäfer, D. Dörr, R. Zheng, N. Naouar, L. Kärge, A hyperelastic approach for modeling the membrane behavior in finite element forming simulation of unidirectional non-crimp fabrics (UD-NCF), Composites Part A: Applied Science and Manufacturing (2023).

[23] Q. Steer, J. Colmars, N. Naouar, P. Boisse, Modeling and analysis of in-plane bending in fibrous reinforcements with rotation-free shell finite elements., International Journal of Solids and Structures 222-223 (2021). https://doi.org/10.1016/j.ijsolstr.2021.03.001

[24] Y. Aimène, E. Vidal-Sallé, B. Hagège, F. Sidoroff, P. Boisse, A Hyperelastic Approach for Composite Reinforcement Large Deformation Analysis, Journal of Composite Materials 44 (2010) 5-26. https://doi.org/10.1177/0021998309345348

[25] C.T. Poppe, C. Krauß, F. Albrecht, L. Kärger, A 3D process simulation model for wet compression moulding, Composites Part A: Applied Science and Manufacturing 50 (2021) 106379. https://doi.org/10.1016/j.compositesa.2021.106379

[26] J.A. White, H. Nishikawa, R.A. Baurle, Weighted Least-squares Cell-Average Gradient Construction Methods For The VULCAN-CFD Second-Order Accurate Unstructured Grid Cell-Centered Finite-Volume Solver, AIAA Scitech 2019 Forum. https://doi.org/10.2514/6.2019-0127

[27] S. Seo, C. Lee, E. Kim, K. Yune, C. Kim, Least-square Switching Process for Accurate and Efficient Gradient Estimation on Unstructured Grid, Journal of the Korean Society for Industrial and Applied Mathematics 24 (2020) 1-22. https://doi.org/10.12941/JKSIAM.2020.24.001

[28] T. Belytschko, W.K. Liu, B. Moran, K.I. Elkhodary, Nonlinear finite elements for continua and structures, 2nd ed., Wiley, Chichester, 2014.
[14] Y. Gong, D. Yan, Y. Yao, R. Wei, H. Hu, P. Xu, X. Peng, An Anisotropic Hyperelastic Constitutive Model with Tension-Shear Coupling for Woven Composite Reinforcements, Int. J. Appl. Mechanics 09 (2017). https://doi.org/10.1142/S1758825117500831

[15] Y. Yao, X. Huang, X. Peng, P. Liu, G. Youkun, An anisotropic hyperelastic constitutive model for plain weave fabric considering biaxial tension coupling, Textile Research Journal 89 (2019) 434-444. https://doi.org/10.1177/0040517517748495

[16] S. Lomov, D. Ivanov, I. Verpoest, M. Zako, T. Kurashiki, H. Nakai, S. Hirosawa, Meso-FE modelling of textile composites: Road map, data flow and algorithms, Composites Science and Technology 67 (2007) 1870-1891. https://doi.org/10.1016/j.compscitech.2006.10.017

[17] J. Sirtautas, A.K. Pickett, P. Lépicier, A mesoscopic model for coupled drape-infusion simulation of biaxial Non-Crimp Fabric, Composites Part B: Engineering 47 (2013) 48-57. https://doi.org/10.1016/j.compositesb.2012.09.088

[18] G. Creech, A. Pickett, Meso-modelling of Non-Crimp Fabric composites for coupled drape and failure analysis, Journal of Materials Science 41 (2006) 6725-6736. https://doi.org/10.1007/s10853-006-0213-6

[19] F.J. Schirmaier, D. Dörr, F. Henning, L. Kärger, A macroscopic approach to simulate the forming behaviour of stitched unidirectional non-crimp fabrics (UD-NCF), Composites Part A: Applied Science and Manufacturing 102 (2017) 322-335. https://doi.org/10.1016/j.compositesa.2017.08.009

[20] B. Schäfer, S. Haas, P. Boisse, L. Kärger, Investigation of the Membrane Behavior of UD-NCF in Macroscopic Forming Simulations, Key Engineering Materials 926 (2022) 1413-1422. https://doi.org/10.4028/p-2977b4

[21] V.N. Khiêm, H. Krieger, M. Itskov, T. Gries, S.E. Stapleton, An averaging based hyperelastic modeling and experimental analysis of non-crimp fabrics, International Journal of Solids and Structures. https://doi.org/10.1016/j.ijsolstr.2016.12.018

[22] B. Schäfer, D. Dörr, R. Zheng, N. Naouar, L. Kärge, A hyperelastic approach for modeling the membrane behavior in finite element forming simulation of unidirectional non-crimp fabrics (UD-NCF), Composites Part A: Applied Science and Manufacturing (2023).

[23] Q. Steer, J. Colmars, N. Naouar, P. Boisse, Modeling and analysis of in-plane bending in fibrous reinforcements with rotation-free shell finite elements., International Journal of Solids and Structures 222-223 (2021). https://doi.org/10.1016/j.ijsolstr.2021.03.001

[24] Y. Aimène, E. Vidal-Sallé, B. Hagège, F. Sidoroff, P. Boisse, A Hyperelastic Approach for Composite Reinforcement Large Deformation Analysis, Journal of Composite Materials 44 (2010) 5-26. https://doi.org/10.1177/0021998309345348

[25] C.T. Poppe, C. Krauß, F. Albrecht, L. Kärger, A 3D process simulation model for wet compression moulding, Composites Part A: Applied Science and Manufacturing 50 (2021) 106379. https://doi.org/10.1016/j.compositesa.2021.106379

[26] J.A. White, H. Nishikawa, R.A. Baurle, Weighted Least-squares Cell-Average Gradient Construction Methods For The VULCAN-CFD Second-Order Accurate Unstructured Grid Cell-Centered Finite-Volume Solver, AIAA Scitech 2019 Forum. https://doi.org/10.2514/6.2019-0127

[27] S. Seo, C. Lee, E. Kim, K. Yune, C. Kim, Least-square Switching Process for Accurate and Efficient Gradient Estimation on Unstructured Grid, Journal of the Korean Society for Industrial and Applied Mathematics 24 (2020) 1-22. https://doi.org/10.12941/JKSIAM.2020.24.001

[28] T. Belytschko, W.K. Liu, B. Moran, K.I. Elkhodary, Nonlinear finite elements for continua and structures, 2nd ed., Wiley, Chichester, 2014.
# Improving precision in the cavity thickness measurement for fabric compaction

SOUSA Pedro<fnref n="1" />,³a*, LIU Xiao<fnref n="2" />,ᵇ, LOMOV Stepan V.<fnref n="1" />,ᶜ and IVENS Jan<fnref n="3" />,ᵈ

**Keywords:** Fabrics/Textiles, Compressibility, Mechanical Testing, Thickness Measurement Methods

**Abstract.** As seen in the first international compressibility benchmark, inaccuracies in the fabric stack thickness measurement, the approach to compliance correction and the non-parallelism between compaction plates resulted in highly inaccurate compression curves. The factors influencing the accuracy of thickness measurements in compaction tests were studied to provide a comprehensive guidance and measurement recommendations to improve upon the current procedures and enhance the accuracy of thickness measurements. These include the variability in the direct thickness measurements due to accuracy of the laser sensors measurements and the variability in the machine compliance measurements. In conclusion, while both thickness measurement methods yielded comparable results, it is noteworthy that the direct method exhibited greater variability in thickness measurements compared to the indirect method, which rely on the machine displacement. Minor changes in the rig's displacement or in the orientations between compaction plates due to differences in the pressure distribution during compaction, combined with other sources of variability such as external interferences or vibrations from the compaction plate led to variations in measurement's accuracy throughout the tests.

## Introduction

A common element to the majority of LCM (liquid composite molding) manufacturing processes is the compaction of a fabric reinforcement inside a rigid or flexible mold prior to the resin injection. Compaction increases the fiber volume fraction and defines the microstructure of the reinforcement for injection. This process subsequently impacts the permeability and impregnation behavior, ultimately determining the mechanical properties of the finished component [1]. Knowing the compaction behavior of a fabric reinforcement across a wide range of pressures provides valuable information for the manufacturing processes. Particularly in the design stage, knowing the compaction behavior of a fabric reinforcement – the stack thickness as a function of pressure – facilitates the selection of the manufacturing process parameters and mold design. The compaction behavior of a reinforcement is determined by compacting a fabric stack between two rigid surfaces in a UTM (universal testing machine) at a constant displacement rate while monitoring the compaction force and the distance between plates, the cavity gap. However, as was seen in the first international benchmark on textile compaction response promoted by NPL (National Physics Laboratory, UK) in the years 2017 and 2018, the variability in these measurements is considerable [2]. Alike the results obtained in the radial permeability benchmark [3] which ran simultaneously with the compressibility benchmark, the scatter in the compression

<figure><img src="image_446.png" /></figure>

<footnote n="1">¹ Department of Materials Engineering, KU Leuven, Kasteelpark Arenberg 44, BE-3001 Leuven, Belgium</footnote>

<footnote n="2">² School of Information Engineering, Wuhan University of Technology, Wuhan 430070, China</footnote>

<footnote n="3">³ Department of Materials Engineering, KU Leuven – Campus De Nayer, J. De Nayerlaan 5, B-2870 Sint-Katelijne Waver, Belgium</footnote>

<footnote>a pedro.sousa@kuleuven.be, b liuxiao001@whut.edu.cn, c stepan.lomov@kuleuven.be,</footnote>

<footnote>d jan.ivens@kuleuven.be</footnote>
# Improving precision in the cavity thickness measurement for fabric compaction

SOUSA Pedro<fnref n="1" />,³a*, LIU Xiao<fnref n="2" />,ᵇ, LOMOV Stepan V.<fnref n="1" />,ᶜ and IVENS Jan<fnref n="3" />,ᵈ

**Keywords:** Fabrics/Textiles, Compressibility, Mechanical Testing, Thickness Measurement Methods

**Abstract.** As seen in the first international compressibility benchmark, inaccuracies in the fabric stack thickness measurement, the approach to compliance correction and the non-parallelism between compaction plates resulted in highly inaccurate compression curves. The factors influencing the accuracy of thickness measurements in compaction tests were studied to provide a comprehensive guidance and measurement recommendations to improve upon the current procedures and enhance the accuracy of thickness measurements. These include the variability in the direct thickness measurements due to accuracy of the laser sensors measurements and the variability in the machine compliance measurements. In conclusion, while both thickness measurement methods yielded comparable results, it is noteworthy that the direct method exhibited greater variability in thickness measurements compared to the indirect method, which rely on the machine displacement. Minor changes in the rig's displacement or in the orientations between compaction plates due to differences in the pressure distribution during compaction, combined with other sources of variability such as external interferences or vibrations from the compaction plate led to variations in measurement's accuracy throughout the tests.

## Introduction

A common element to the majority of LCM (liquid composite molding) manufacturing processes is the compaction of a fabric reinforcement inside a rigid or flexible mold prior to the resin injection. Compaction increases the fiber volume fraction and defines the microstructure of the reinforcement for injection. This process subsequently impacts the permeability and impregnation behavior, ultimately determining the mechanical properties of the finished component [1]. Knowing the compaction behavior of a fabric reinforcement across a wide range of pressures provides valuable information for the manufacturing processes. Particularly in the design stage, knowing the compaction behavior of a fabric reinforcement – the stack thickness as a function of pressure – facilitates the selection of the manufacturing process parameters and mold design. The compaction behavior of a reinforcement is determined by compacting a fabric stack between two rigid surfaces in a UTM (universal testing machine) at a constant displacement rate while monitoring the compaction force and the distance between plates, the cavity gap. However, as was seen in the first international benchmark on textile compaction response promoted by NPL (National Physics Laboratory, UK) in the years 2017 and 2018, the variability in these measurements is considerable [2]. Alike the results obtained in the radial permeability benchmark [3] which ran simultaneously with the compressibility benchmark, the scatter in the compression

<figure><img src="image_457.png" /></figure>

<footnote n="1">¹ Department of Materials Engineering, KU Leuven, Kasteelpark Arenberg 44, BE-3001 Leuven, Belgium</footnote>

<footnote n="2">² School of Information Engineering, Wuhan University of Technology, Wuhan 430070, China</footnote>

<footnote n="3">³ Department of Materials Engineering, KU Leuven – Campus De Nayer, J. De Nayerlaan 5, B-2870 Sint-Katelijne Waver, Belgium</footnote>

<footnote>a pedro.sousa@kuleuven.be, b liuxiao001@whut.edu.cn, c stepan.lomov@kuleuven.be,</footnote>

<footnote>d jan.ivens@kuleuven.be</footnote>
curves reported by different participants was much higher than the scatter in individual measurements. Though the variability in the compression curves reported by each participant was consistent, a scatter of 38% among participants was found for the compaction pressure at the target thickness of 3 mm. In saturated conditions, the scatter among participants was 40% for the Saertex fabric (NCF multiaxial) and 50% for the Hexcel fabric (2/2 twill woven) both E-glass fabrics. The scatter in the compaction pressure was attributed to inaccuracies in the thickness measurements, the approach to the compliance correction and the non-parallelism between compaction plates. Such results showed the need for improved compressibility guidelines regarding the thickness measurement technique and compliance corrections. These challenges were tackled in a new methodology developed by the benchmark consortium, with participation of the authors, for the second international benchmarking exercise in 2022 [4]. The second compressibility benchmark, also promoted by the NPL, intended to increase the accuracy in compression test and refine the test procedure for a future standardization through the implementation of recommendations for the measurements. In this benchmark, the usage of direct thickness measurement methods – laser sensors or LVDTs – was mandatory in addition to self-aligning plates with a lockable mechanism. However, several aspects such as the placement of distance sensors in the rig, the usage of self-aligning plates, and the approach to the compliance correction require a deeper study because these can still compromise the measurement’s accuracy.

This paper presents a methodology for obtaining and analyzing data in compression tests with simultaneous use of direct (laser sensors) and indirect (UTM) thickness measurement methods. This study shows how several variability sources related to the measurement's methodology can be decreased to narrow the differences in the accuracy of both measurement methods. Good measurement practices are proposed to improve the measurement's accuracy, aiming to the establishment of a standardized and consistent approach to fabric compaction measurements.

## The implementation of the testing protocol

### The compressibility rig

The compaction tests were performed on a radial-permeability rig [5] installed in a UTM, an Instron 4467 displacement-controlled machine with a load cell capacity of 30 kN (±0.5%). The compaction rig is comprised of an upper compaction plate, a cylindrical aluminum plate 250 mm in diameter, and a lower plate, a square glass plate (Herasil 102) with dimension of 300 mm x 300 mm. Both plates are 30-mm thick as shown in Figure 1. A ball-and-socket joint system attached to the upper compaction plate allows the free rotation of this plate, ensuring the correct parallelism between plates when the plates are pressed against each other. Once the desired load is achieved, the ball joint is locked in place by a threaded clamp. The distance between compaction plates, the cavity gap, is measured by three laser sensors (UR, UB and UF), model CP08MHT80 from Wenglor®, resolution of 8 µm, installed on the top of the upper compaction plate, 90° out of phase with each other. Red electrical and aluminum tapes were placed on the glass plate to improve the laser back reflection. The aluminum tape removes the plate's transparency, and the electrical tape facilitates the cleaning and its replacement during tests. In saturated tests, a 4-mm thick layer of masking tape was placed in front of the reflective tapes to prevent oil from covering these tapes and interfering with the laser's reflection. The sensors output is digitized by a USB-6210, a multifunction I/O Device from National Instruments programmed in LabVIEW. When the pre-set cavity thickness is detected by two of the three laser sensors, the UTM's crosshead movement direction is automatically reverted. The arrangement around the compaction plate results in three compression curves measured on the right, back, and front sides of the rig, allowing to define a plane on the compaction plate surface and monitor its parallelism with the glass plate below.
curves reported by different participants was much higher than the scatter in individual measurements. Though the variability in the compression curves reported by each participant was consistent, a scatter of 38% among participants was found for the compaction pressure at the target thickness of 3 mm. In saturated conditions, the scatter among participants was 40% for the Saertex fabric (NCF multiaxial) and 50% for the Hexcel fabric (2/2 twill woven) both E-glass fabrics. The scatter in the compaction pressure was attributed to inaccuracies in the thickness measurements, the approach to the compliance correction and the non-parallelism between compaction plates. Such results showed the need for improved compressibility guidelines regarding the thickness measurement technique and compliance corrections. These challenges were tackled in a new methodology developed by the benchmark consortium, with participation of the authors, for the second international benchmarking exercise in 2022 [4]. The second compressibility benchmark, also promoted by the NPL, intended to increase the accuracy in compression test and refine the test procedure for a future standardization through the implementation of recommendations for the measurements. In this benchmark, the usage of direct thickness measurement methods – laser sensors or LVDTs – was mandatory in addition to self-aligning plates with a lockable mechanism. However, several aspects such as the placement of distance sensors in the rig, the usage of self-aligning plates, and the approach to the compliance correction require a deeper study because these can still compromise the measurement's accuracy.

This paper presents a methodology for obtaining and analyzing data in compression tests with simultaneous use of direct (laser sensors) and indirect (UTM) thickness measurement methods. This study shows how several variability sources related to the measurement's methodology can be decreased to narrow the differences in the accuracy of both measurement methods. Good measurement practices are proposed to improve the measurement's accuracy, aiming to the establishment of a standardized and consistent approach to fabric compaction measurements.

## The implementation of the testing protocol

### The compressibility rig

The compaction tests were performed on a radial-permeability rig [5] installed in a UTM, an Instron 4467 displacement-controlled machine with a load cell capacity of 30 kN (±0.5%). The compaction rig is comprised of an upper compaction plate, a cylindrical aluminum plate 250 mm in diameter, and a lower plate, a square glass plate (Herasil 102) with dimension of 300 mm x 300 mm. Both plates are 30-mm thick as shown in Figure 1. A ball-and-socket joint system attached to the upper compaction plate allows the free rotation of this plate, ensuring the correct parallelism between plates when the plates are pressed against each other. Once the desired load is achieved, the ball joint is locked in place by a threaded clamp. The distance between compaction plates, the cavity gap, is measured by three laser sensors (UR, UB and UF), model CP08MHT80 from Wenglor®, resolution of 8 µm, installed on the top of the upper compaction plate, 90° out of phase with each other. Red electrical and aluminum tapes were placed on the glass plate to improve the laser back reflection. The aluminum tape removes the plate's transparency, and the electrical tape facilitates the cleaning and its replacement during tests. In saturated tests, a 4-mm thick layer of masking tape was placed in front of the reflective tapes to prevent oil from covering these tapes and interfering with the laser's reflection. The sensors output is digitized by a USB-6210, a multifunction I/O Device from National Instruments programmed in LabVIEW. When the pre-set cavity thickness is detected by two of the three laser sensors, the UTM's crosshead movement direction is automatically reverted. The arrangement around the compaction plate results in three compression curves measured on the right, back, and front sides of the rig, allowing to define a plane on the compaction plate surface and monitor its parallelism with the glass plate below.
## Materials and methods

The compression tests were performed with a Saertex multiaxial NCF E-glass reinforcement centrally distributed to the benchmark participants. The fabric stacks were compacted up to a target thickness of 3 mm at a constant speed rate of 1 mm/min in dry and saturated conditions on five specimens. At the target thickness of 3 mm, as measured by the direct thickness method, the compaction process was reverted, and the stack unloaded at the same displacement rate. Square fabric samples with side 150 mm (± 1.9 mm) were cut from the fabric roll with a fabric shears and stacks of 10 layers were formed. All layers in the stacks had the 0° (warp) direction coinciding and the same surface facing upwards. The fabrics dimensions are smaller than the compression surface area as requested in the measurement procedure; therefore, the load values measured by the UTM were converted to applied pressure by dividing them by the area of the specimen surface, which has a value of 150±1.9 cm². An areal density of 452 ± 3.5 g/m² was determined by measuring the area of each fabric layer with a ruler and weighed in a Mettler AE 24 scale (± 0.1 gram).

<figure><img src="image_447.png" /><figcaption>Figure 1: The rig's installation in an UTM allows to compress a fabric stack at a constant displacement rate up to a pre-set pressure or thickness. The upper compaction plate, a cylindrical aluminum plate 250 mm in diameter, is self-aligning with the glass plate below by a ball-and-socket joint. The cavity gap measured by the indirect thickness method is obtained from the crosshead position with a correction for the rig's compliance, while the direct thickness method measures the cavity gap by calculating the difference between a fixed and movable reference point measured outside the cavity.</figcaption></figure>

The compaction in saturated conditions was conducted by submersing an entire fabric stack in a silicone oil bath, Dow Corning® Xiameter PMX 200/100 cS, 0.96 kg/l, for 15 min at an oil bath temperature of 22.4 ± 0.3 °C and a room temperature of 23.1 ± 0.2 °C. The temperature within the rig cavity was measured before placing a fabric stack, revealing an average of 23.8 ± 0.3 °C. At these temperatures, the oil viscosity ranges between 93.7 ± 0.5 mPa.s and 94.2 ± 0.5 mPa.s at the lowest and highest registered temperature, respectively. The temperature in the rig cavity is slightly
## Materials and methods

The compression tests were performed with a Saertex multiaxial NCF E-glass reinforcement centrally distributed to the benchmark participants. The fabric stacks were compacted up to a target thickness of 3 mm at a constant speed rate of 1 mm/min in dry and saturated conditions on five specimens. At the target thickness of 3 mm, as measured by the direct thickness method, the compaction process was reverted, and the stack unloaded at the same displacement rate. Square fabric samples with side 150 mm (± 1.9 mm) were cut from the fabric roll with a fabric shears and stacks of 10 layers were formed. All layers in the stacks had the 0° (warp) direction coinciding and the same surface facing upwards. The fabrics dimensions are smaller than the compression surface area as requested in the measurement procedure; therefore, the load values measured by the UTM were converted to applied pressure by dividing them by the area of the specimen surface, which has a value of 150±1.9 cm². An areal density of 452 ± 3.5 g/m² was determined by measuring the area of each fabric layer with a ruler and weighed in a Mettler AE 24 scale (± 0.1 gram).

<figure><img src="image_458.png" /><figcaption>Figure 1: The rig's installation in an UTM allows to compress a fabric stack at a constant displacement rate up to a pre-set pressure or thickness. The upper compaction plate, a cylindrical aluminum plate 250 mm in diameter, is self-aligning with the glass plate below by a ball-and-socket joint. The cavity gap measured by the indirect thickness method is obtained from the crosshead position with a correction for the rig's compliance, while the direct thickness method measures the cavity gap by calculating the difference between a fixed and movable reference point measured outside the cavity.</figcaption></figure>

The compaction in saturated conditions was conducted by submersing an entire fabric stack in a silicone oil bath, Dow Corning® Xiameter PMX 200/100 cS, 0.96 kg/l, for 15 min at an oil bath temperature of 22.4 ± 0.3 °C and a room temperature of 23.1 ± 0.2 °C. The temperature within the rig cavity was measured before placing a fabric stack, revealing an average of 23.8 ± 0.3 °C. At these temperatures, the oil viscosity ranges between 93.7 ± 0.5 mPa.s and 94.2 ± 0.5 mPa.s at the lowest and highest registered temperature, respectively. The temperature in the rig cavity is slightly
higher than the room temperature due to heat being generated by the machine electronics, placed below the machine bench, and conducted upwards through the rig.

# Results

## Machine compliance measurements

The usage of compaction plates with locking mechanisms is essential to ensure the correct parallelism between plates during compression because non-locked compaction plates do not retain the parallelism during compaction. The plates were pressed together at a displacement rate of 0.5 mm/min up and locked at a load of 28 kN (0.57 MPa), the maximum pressure allowed by the system. The locking load is higher than the load encountered during compression tests to prevent a plate realignment. This load criterion is important because the compaction load lowers the clamping pressure in the ball-joint which allows the plate to readjust its position. Based on the observations made during the tests, it is recommended that the locking load be set at a minimum of 25% higher than the compaction load to prevent readjustments during compaction. However, locking the compaction plate at a greater load carries the risk of calibrating the rig at a higher deformation level, which can impact the zero-thickness point in direct methods and compromise the accuracy of the calibration. To ensure that the displacement measured by the laser sensors remains constant as the pressure increases, the machine compliance was measured by the laser sensors and by the UTM up to the locking pressure of 0.57 MPa. Before registering the first compliance curves with the laser sensors, the rig is compressed 10 times in a cyclic manner to eliminate the non-sample displacement in the loading system resulting from squeezing the self-aligning ball joint.

The machine compliance measured by the UTM is characterized by a nonlinearity below 0.1 MPa, followed by a linear region where the compliance increases at a constant factor, and again nonlinearity above 0.4 MPa, see Figure 2A. The initial nonlinearity is mainly due to small tolerances in the threads and pin joints attached to the UTM crosshead that slow down the load increase. Only when all parts are in their lowest position and the contact between plates is established, the load shows a linear increase with displacement. In terms of repeatability, the difference between successive compliance curves was below 2 µm. The maximum difference between the initial and the compliance after tests was of 3.5 ± 0.24 µm, registered under dry conditions. Overall, the compliance registered by the UTM remained constant through tests.

The machine compliance measured by the laser sensors showed a constant displacement across the pressure range (< 4 µm) for sensors UR and UB, while sensor UF measured a displacement decrease of 12 µm up to a pressure of 0.1 MPa (5 kN) before remaining constant as the two other sensors, see Figure 2B. The maximum increase in displacement between the lowest pressure of 0.03 MPa and a pressure of 0.1 MPa was 20 µm, remaining constant above 0.1 MPa up to the maximum pressure of 0.6 MPa. The difference among the three curves implies that the two compaction plates are parallel within a tolerance of 20 µm. A constant displacement with pressure was only achieved with the laser sensors installed on the top of the compaction plate, 10 mm from the plate, because the rig's deformation is smaller in this region in comparison to other regions of the rig further from the compaction plate. Additionally, it is important to place the sensors close to the cavity gap to minimize the measurement distance, thus reducing the linearity error associated with the measurement distance. Because the compliance curves registered by the laser sensors showed a constant displacement with pressure, the zero-thickness point necessary for direct thickness measurements can be registered at any pressure value. However, as can be seen in Figure 2B, sensor UF one of the compliance curves shows a major deviation in the measured displacement. In dry conditions, the maximum magnitude of these "shifts" was 25 µm, while in saturated conditions a magnitude of 50 µm was registered after a series of compression tests. These shifts in displacement can either be constant along pressure, decrease below the sensor's accuracy limit at a higher pressure, or also being below the sensor's accuracy limit for the entire pressure
higher than the room temperature due to heat being generated by the machine electronics, placed below the machine bench, and conducted upwards through the rig.

## Results

### Machine compliance measurements

The usage of compaction plates with locking mechanisms is essential to ensure the correct parallelism between plates during compression because non-locked compaction plates do not retain the parallelism during compaction. The plates were pressed together at a displacement rate of 0.5 mm/min up and locked at a load of 28 kN (0.57 MPa), the maximum pressure allowed by the system. The locking load is higher than the load encountered during compression tests to prevent a plate realignment. This load criterion is important because the compaction load lowers the clamping pressure in the ball-joint which allows the plate to readjust its position. Based on the observations made during the tests, it is recommended that the locking load be set at a minimum of 25% higher than the compaction load to prevent readjustments during compaction. However, locking the compaction plate at a greater load carries the risk of calibrating the rig at a higher deformation level, which can impact the zero-thickness point in direct methods and compromise the accuracy of the calibration. To ensure that the displacement measured by the laser sensors remains constant as the pressure increases, the machine compliance was measured by the laser sensors and by the UTM up to the locking pressure of 0.57 MPa. Before registering the first compliance curves with the laser sensors, the rig is compressed 10 times in a cyclic manner to eliminate the non-sample displacement in the loading system resulting from squeezing the self-aligning ball joint.

The machine compliance measured by the UTM is characterized by a nonlinearity below 0.1 MPa, followed by a linear region where the compliance increases at a constant factor, and again nonlinearity above 0.4 MPa, see Figure 2A. The initial nonlinearity is mainly due to small tolerances in the threads and pin joints attached to the UTM crosshead that slow down the load increase. Only when all parts are in their lowest position and the contact between plates is established, the load shows a linear increase with displacement. In terms of repeatability, the difference between successive compliance curves was below 2 µm. The maximum difference between the initial and the compliance after tests was of 3.5 ± 0.24 µm, registered under dry conditions. Overall, the compliance registered by the UTM remained constant through tests.

The machine compliance measured by the laser sensors showed a constant displacement across the pressure range (< 4 µm) for sensors UR and UB, while sensor UF measured a displacement decrease of 12 µm up to a pressure of 0.1 MPa (5 kN) before remaining constant as the two other sensors, see Figure 2B. The maximum increase in displacement between the lowest pressure of 0.03 MPa and a pressure of 0.1 MPa was 20 µm, remaining constant above 0.1 MPa up to the maximum pressure of 0.6 MPa. The difference among the three curves implies that the two compaction plates are parallel within a tolerance of 20 µm. A constant displacement with pressure was only achieved with the laser sensors installed on the top of the compaction plate, 10 mm from the plate, because the rig's deformation is smaller in this region in comparison to other regions of the rig further from the compaction plate. Additionally, it is important to place the sensors close to the cavity gap to minimize the measurement distance, thus reducing the linearity error associated with the measurement distance. Because the compliance curves registered by the laser sensors showed a constant displacement with pressure, the zero-thickness point necessary for direct thickness measurements can be registered at any pressure value. However, as can be seen in Figure 2B, sensor UF one of the compliance curves shows a major deviation in the measured displacement. In dry conditions, the maximum magnitude of these "shifts" was 25 µm, while in saturated conditions a magnitude of 50 µm was registered after a series of compression tests. These shifts in displacement can either be constant along pressure, decrease below the sensor's accuracy limit at a higher pressure, or also being below the sensor's accuracy limit for the entire pressure
range. The source of these shifts can either be related to small readjustments to the rig's components or in the orientation between compaction plates; however, the randomness of these deviations did not allow to pin-point its source. To comprehend this variability, further investigation will be necessary.

<figure><img src="image_448.png" /><figcaption>Figure 2: Machine compliance measurements A) Compliance measured by the UTM B) Compliance measured by the laser sensors in three locations around the compaction plate.</figcaption></figure>

The compliance curves registered before and after the compression tests were compared to access their variability. Sensor UB registered an average difference of 19.5 ± 2.4 µm and 35.4 ± 1.2 µm in dry and saturated conditions, respectively. Sensor UF, on the other hand, recorded a difference of 43.4 ± 10.5 µm in dry conditions and 1.6 ± 2 µm in saturated conditions. Finally, Sensor UR exhibited the lowest difference, with 1.3 ± 1 µm in dry conditions and 9.8 ± 1.7 µm in saturated conditions.

Overall, the compliance measurements taken with the laser sensors exhibited a larger variability than those recorded by the UTM. This discrepancy can be attributed to the sensors' higher sensitivity to minor changes in the plate's orientation, as well as their susceptibility to external interferences or vibrations from the compaction plate's movement. On the other hand, the UTM's displacement measurements are less susceptible to such interferences and are unaffected by the plate's orientation.

To account for the variability in the compliance measured by the laser sensors during compression, the zero-thickness point was registered and changed after each compression test. The zero-thickness point was obtained from compliance curves with a constant displacement along pressure to ensure its accuracy. In dry conditions, the zero-thickness points remained within the variability range of 20 µm, which was consistent with the variability observed in the compliance measurements. In saturated conditions a 40-µm deviation was observed immediately after the first zero-thickness point, remaining constant (within 10 µm) until the last measurement.

Currently, no guidelines exist to help reduce the variability associated with direct methods beyond frequent compliance measurements. As a result, it is important to keep in mind the potential limitations of using laser sensors for compliance measurements. In contrast, UTM displacement measurements may provide more stable results; however, it is worth noting that UTM measurements may not fully capture the changes in compliance that occur during compression, because the laser sensors have demonstrated higher sensitivity to such changes.

### Compression curves

The compression curves obtained with the laser sensors were processed with a curve-fit to smooth out oscillations in the curves, see Figure 3A. These oscillations have a maximum amplitude of 20
range. The source of these shifts can either be related to small readjustments to the rig's components or in the orientation between compaction plates; however, the randomness of these deviations did not allow to pin-point its source. To comprehend this variability, further investigation will be necessary.

<figure><img src="image_459.png" /><figcaption>Figure 2: Machine compliance measurements A) Compliance measured by the UTM B) Compliance measured by the laser sensors in three locations around the compaction plate.</figcaption></figure>

The compliance curves registered before and after the compression tests were compared to access their variability. Sensor UB registered an average difference of 19.5 ± 2.4 µm and 35.4 ± 1.2 µm in dry and saturated conditions, respectively. Sensor UF, on the other hand, recorded a difference of 43.4 ± 10.5 µm in dry conditions and 1.6 ± 2 µm in saturated conditions. Finally, Sensor UR exhibited the lowest difference, with 1.3 ± 1 µm in dry conditions and 9.8 ± 1.7 µm in saturated conditions.

Overall, the compliance measurements taken with the laser sensors exhibited a larger variability than those recorded by the UTM. This discrepancy can be attributed to the sensors' higher sensitivity to minor changes in the plate's orientation, as well as their susceptibility to external interferences or vibrations from the compaction plate's movement. On the other hand, the UTM's displacement measurements are less susceptible to such interferences and are unaffected by the plate's orientation.

To account for the variability in the compliance measured by the laser sensors during compression, the zero-thickness point was registered and changed after each compression test. The zero-thickness point was obtained from compliance curves with a constant displacement along pressure to ensure its accuracy. In dry conditions, the zero-thickness points remained within the variability range of 20 µm, which was consistent with the variability observed in the compliance measurements. In saturated conditions a 40-µm deviation was observed immediately after the first zero-thickness point, remaining constant (within 10 µm) until the last measurement.

Currently, no guidelines exist to help reduce the variability associated with direct methods beyond frequent compliance measurements. As a result, it is important to keep in mind the potential limitations of using laser sensors for compliance measurements. In contrast, UTM displacement measurements may provide more stable results; however, it is worth noting that UTM measurements may not fully capture the changes in compliance that occur during compression, because the laser sensors have demonstrated higher sensitivity to such changes.

### Compression curves

The compression curves obtained with the laser sensors were processed with a curve-fit to smooth out oscillations in the curves, see Figure 3A. These oscillations have a maximum amplitude of 20
µm, measured as the difference between the displacement in the curves and the displacement provided by the curve fitting, see Figure 3B. These oscillations are likely due to a wobble in the plate's downward movement because it disappears above a pressure of 0.14 MPa, appearing again in the unloading branch below a pressure of 0.12 MPa. For that reason, these oscillations are not visible in the compliance measurements when the plates are pressed against each other. Smoothing these oscillations using the curve-fit technique facilitates the further analysis and interpretation of the obtained results. All the subsequent results obtained with the laser sensors are presented in terms of the smoothed curves.

<figure><img src="image_449.png" /></figure>

<figure><img src="image_450.png" /><figcaption>Figure 3: The compression curves obtained by the laser sensors showed interferences in the form of oscillations A) These oscillations were smoothed-out with a curve fit B) Difference between the fit-curve and the compression curve in both loading and unloading curve branches.</figcaption></figure>

The direct thickness measurements are obtained by the thickness readings of three laser sensors, leading to three compaction curves. When these curves are superimposed, instead of three overlapped curves, the results show two of the curves overlapped while the third one is offset, see Figure 4A. The offset curve was either UB, UF, or UR, changing on a sample-to-sample basis which indicates that the compaction plate settles in a slightly different manner with each compression. The average overlap and offset between curves were estimated from the differences in thickness at the same compaction pressure. Under dry conditions, the two curves overlapped by an average of 14.7±8.9 µm during loading and 22.9±6.2 µm during unloading. Under saturated conditions an average overlapping of 21.6±13 µm was obtained during loading, while the overlapping during unloading was 29.2±14 µm. The overlapping between curves falls within the uncertainty found in the compliance measurements registered by the laser sensors, therefore within the achievable accuracy with this rig.

Regarding the offset curve, on average, it was shifted on average 51.2±16.4 µm from the furthest curve, while in saturated conditions, the offset was of 61.2±22 µm. A similar offset was detected in the compliance measurements registered by the laser sensors; however, the offset curve was constantly UF, and the magnitude of this deviation was just 12±1.9 µm. This offset suggests that the compaction plates have a certain level of mobility, resulting in more displacement when the glass plate is centrally loaded with fabrics.

Despite the common use of both direct and indirect thickness measurement methods, the results from these methods are often not reported together and directly compared in the literature. Consequently, the variability affecting these methods was not yet the target of a study, leading to inaccuracies when results obtained by these methods are compared. The average thickness determined by three laser sensors was compared to that obtained through the indirect method, aiming to eliminate the influence of the platen tilting and determine the difference between the two methods.
µm, measured as the difference between the displacement in the curves and the displacement provided by the curve fitting, see Figure 3B. These oscillations are likely due to a wobble in the plate's downward movement because it disappears above a pressure of 0.14 MPa, appearing again in the unloading branch below a pressure of 0.12 MPa. For that reason, these oscillations are not visible in the compliance measurements when the plates are pressed against each other. Smoothing these oscillations using the curve-fit technique facilitates the further analysis and interpretation of the obtained results. All the subsequent results obtained with the laser sensors are presented in terms of the smoothed curves.

<figure><img src="image_460.png" /></figure>

<figure><img src="image_461.png" /><figcaption>Figure 3: The compression curves obtained by the laser sensors showed interferences in the form of oscillations A) These oscillations were smoothed-out with a curve fit B) Difference between the fit-curve and the compression curve in both loading and unloading curve branches.</figcaption></figure>

The direct thickness measurements are obtained by the thickness readings of three laser sensors, leading to three compaction curves. When these curves are superimposed, instead of three overlapped curves, the results show two of the curves overlapped while the third one is offset, see Figure 4A. The offset curve was either UB, UF, or UR, changing on a sample-to-sample basis which indicates that the compaction plate settles in a slightly different manner with each compression. The average overlap and offset between curves were estimated from the differences in thickness at the same compaction pressure. Under dry conditions, the two curves overlapped by an average of 14.7±8.9 µm during loading and 22.9±6.2 µm during unloading. Under saturated conditions an average overlapping of 21.6±13 µm was obtained during loading, while the overlapping during unloading was 29.2±14 µm. The overlapping between curves falls within the uncertainty found in the compliance measurements registered by the laser sensors, therefore within the achievable accuracy with this rig.

Regarding the offset curve, on average, it was shifted on average 51.2±16.4 µm from the furthest curve, while in saturated conditions, the offset was of 61.2±22 µm. A similar offset was detected in the compliance measurements registered by the laser sensors; however, the offset curve was constantly UF, and the magnitude of this deviation was just 12±1.9 µm. This offset suggests that the compaction plates have a certain level of mobility, resulting in more displacement when the glass plate is centrally loaded with fabrics.

Despite the common use of both direct and indirect thickness measurement methods, the results from these methods are often not reported together and directly compared in the literature. Consequently, the variability affecting these methods was not yet the target of a study, leading to inaccuracies when results obtained by these methods are compared. The average thickness determined by three laser sensors was compared to that obtained through the indirect method, aiming to eliminate the influence of the platen tilting and determine the difference between the two methods.
<figure><img src="image_451.png" /></figure>

<figure><img src="image_452.png" /><figcaption>Figure 4: Compression in dry conditions with three laser sensors A) The compression curves as measured by the laser sensors showed two overlapped compression curves while one of them is offset. B) The three compression curves were averaged, resulting a single direct thickness curve, and compared with the thickness measured by the indirect measurement method corrected for the machine compliance.</figcaption></figure>

Overall, both methods provided similar measurements within the accuracy allowed by the laser sensors, see Figure 4B. The average differences between methods stayed below 35 µm for all tests except one test in dry conditions which showed an average difference of 38.6±4.2 µm between curves, see Figure 5. While the compliance curves showed greater variability under saturated conditions, this variability did not impact the thickness measurements as frequent compliance measurements were taken to mitigate potential discrepancies.

<figure><img src="image_453.png" /><figcaption>Figure 5: Difference between the thickness obtained through direct and indirect methods at the same compaction pressure.</figcaption></figure>

In dry conditions the curves obtained by the direct method are on average 22.4±10 µm shifted up in relation to the indirect method, see Figure 6. No specific trend was observed in the difference between curves because the difference between methods can either be constant along the pressure or decrease with pressure. The thickness measured at the maximum compaction pressure showed a comparable level of accuracy between the direct method (3.00±0.02 mm) and the indirect method (2.99±0.01 mm). Under saturated conditions, both methods provided the same thickness within an accuracy of 1.7%. The compression in saturated conditions showed little effect on the compaction process due to the presence of stitches. The presence of stitches reduces the mobility of the fiber bundles and prevents nesting [6]. For that reason, the compression tests performed under saturated conditions showed little difference to those performed in dry conditions. The compaction pressure at the target thickness of 3 mm is similar in both cases, 392.8±38.7 kPa in dry conditions and 359.6±40.9 kPa in saturated conditions.
<figure><img src="image_462.png" /></figure>

<figure><img src="image_463.png" /><figcaption>Figure 4: Compression in dry conditions with three laser sensors A) The compression curves as measured by the laser sensors showed two overlapped compression curves while one of them is offset. B) The three compression curves were averaged, resulting a single direct thickness curve, and compared with the thickness measured by the indirect measurement method corrected for the machine compliance.</figcaption></figure>

Overall, both methods provided similar measurements within the accuracy allowed by the laser sensors, see Figure 4B. The average differences between methods stayed below 35 µm for all tests except one test in dry conditions which showed an average difference of 38.6±4.2 µm between curves, see Figure 5. While the compliance curves showed greater variability under saturated conditions, this variability did not impact the thickness measurements as frequent compliance measurements were taken to mitigate potential discrepancies.

<figure><img src="image_464.png" /><figcaption>Figure 5: Difference between the thickness obtained through direct and indirect methods at the same compaction pressure.</figcaption></figure>

In dry conditions the curves obtained by the direct method are on average 22.4±10 µm shifted up in relation to the indirect method, see Figure 6. No specific trend was observed in the difference between curves because the difference between methods can either be constant along the pressure or decrease with pressure. The thickness measured at the maximum compaction pressure showed a comparable level of accuracy between the direct method (3.00±0.02 mm) and the indirect method (2.99±0.01 mm). Under saturated conditions, both methods provided the same thickness within an accuracy of 1.7%. The compression in saturated conditions showed little effect on the compaction process due to the presence of stitches. The presence of stitches reduces the mobility of the fiber bundles and prevents nesting [6]. For that reason, the compression tests performed under saturated conditions showed little difference to those performed in dry conditions. The compaction pressure at the target thickness of 3 mm is similar in both cases, 392.8±38.7 kPa in dry conditions and 359.6±40.9 kPa in saturated conditions.
<figure><img src="image_454.png" /><figcaption>Figure 6: Compression curves obtained with direct and indirect thickness, Samples 1 to 5.</figcaption></figure>

The variability among compression curves was more significant in saturated conditions. At a pressure of 0.1 MPa, the direct method measured 3.34±0.02 mm, while the indirect method measured a thickness of 3.31±0.04 mm in dry conditions. In saturated conditions, a thickness of 3.32±0.05 mm was measured by the indirect method, while the direct method measured 3.33±0.06 mm. In both cases, the measured thickness was similarly independent of the measurement method.

### *The parallelism between compaction plates*

The parallelism between plates was checked with a Fujifilm Prescale® Ultra Low-pressure pressure-sensitive film. Squared sheets with side 260 mm were cut from the roll and placed in the cavity between compaction surfaces. The cavity was then closed with a pressure of 0.41 MPa, which is close to the specified pressure limit for this paper, set at 0.5 MPa. The cavity was then opened, and the pressure film revealed the pressure distribution profile that occurred between the two surfaces. The color intensity exhibited by the pressure film is directly proportional to the magnitude of the applied pressure, with a visual accuracy of approximately ±10% and a spatial resolution between 5 and 15 microns.

The impression left on the paper showed a uniform and consistent pressure only along the edges of the compaction plate, see Figure 7A. A significant portion of the compaction surface did not exhibit any contact, indicating a concave shape in the round compaction plate. This condition can be related to the surface finish which appears not to have been appropriately flattened and smoothed. The surface roughness was measured with a profilometer resulting in an average Ra of 1.00 µm, which indicates a certain roughness. The compression curves measured by the laser sensors indicated a certain level of mobility when fabrics are compacted. For that reason, the pressure distribution was also checked with fabric layers in the cavity. For this measurement, three layers of the twill 2/2 fabric used in the first compressibility benchmark were placed in the cavity and compressed up to a pressure of 0.41 MPa, see Figure 7B. The pressure distribution shows a slightly darker color on the left side of the paper indicating a different pressure distribution. This difference can be related to the concave shape of the compaction plate; therefore, this verification needs to be repeated once the plate is flattened.

### *Comparison with the benchmark results*

The accuracy of the compression measurements was assessed by a comparison with the results reported by other benchmark participants. While the results from the latest benchmark are yet not available, the fabric is of the same type as the fabric used in the first benchmark exercise [2]. Although this fabric is from a different batch, the results are still comparable to some degree.
<figure><img src="image_465.png" /><figcaption>Figure 6: Compression curves obtained with direct and indirect thickness, Samples 1 to 5.</figcaption></figure>

The variability among compression curves was more significant in saturated conditions. At a pressure of 0.1 MPa, the direct method measured 3.34±0.02 mm, while the indirect method measured a thickness of 3.31±0.04 mm in dry conditions. In saturated conditions, a thickness of 3.32±0.05 mm was measured by the indirect method, while the direct method measured 3.33±0.06 mm. In both cases, the measured thickness was similarly independent of the measurement method.

### *The parallelism between compaction plates*

The parallelism between plates was checked with a Fujifilm Prescale® Ultra Low-pressure pressure-sensitive film. Squared sheets with side 260 mm were cut from the roll and placed in the cavity between compaction surfaces. The cavity was then closed with a pressure of 0.41 MPa, which is close to the specified pressure limit for this paper, set at 0.5 MPa. The cavity was then opened, and the pressure film revealed the pressure distribution profile that occurred between the two surfaces. The color intensity exhibited by the pressure film is directly proportional to the magnitude of the applied pressure, with a visual accuracy of approximately ±10% and a spatial resolution between 5 and 15 microns.

The impression left on the paper showed a uniform and consistent pressure only along the edges of the compaction plate, see Figure 7A. A significant portion of the compaction surface did not exhibit any contact, indicating a concave shape in the round compaction plate. This condition can be related to the surface finish which appears not to have been appropriately flattened and smoothed. The surface roughness was measured with a profilometer resulting in an average Ra of 1.00 µm, which indicates a certain roughness. The compression curves measured by the laser sensors indicated a certain level of mobility when fabrics are compacted. For that reason, the pressure distribution was also checked with fabric layers in the cavity. For this measurement, three layers of the twill 2/2 fabric used in the first compressibility benchmark were placed in the cavity and compressed up to a pressure of 0.41 MPa, see Figure 7B. The pressure distribution shows a slightly darker color on the left side of the paper indicating a different pressure distribution. This difference can be related to the concave shape of the compaction plate; therefore, this verification needs to be repeated once the plate is flattened.

### *Comparison with the benchmark results*

The accuracy of the compression measurements was assessed by a comparison with the results reported by other benchmark participants. While the results from the latest benchmark are yet not available, the fabric is of the same type as the fabric used in the first benchmark exercise [2]. Although this fabric is from a different batch, the results are still comparable to some degree.
<figure><img src="image_455.png" /><figcaption>Figure 7: Contact (A) between compaction plates, and (B) between three layers of Twill 2/2 (Hexcel fabric)</figcaption></figure>

An average thickness of 3.45±0.1 mm at a pressure of 0.1 MPa was reported in the first benchmark exercise under dry conditions, while in saturated conditions, the average thickness was 3.39±0.14 mm. In this study, at a pressure of 0.1 MPa a thickness of 3.31±0.01 mm and 3.34±0.02 mm was measured by the indirect and direct thickness methods, respectively. In saturated conditions, an average thickness of 3.33±0.06 mm was measured at this pressure. No significant differences were observed between the direct and indirect thickness measurement methods. The target thickness of 3 mm was reached at a compaction pressure of 393±39 kPa under dry conditions, while the reported participant's average is 370±141 kPa. Under saturated conditions, the target thickness was achieved at a pressure of 360±41 kPa, which is also comparable to the reported pressure of 317±127 kPa.

## Summary

This study showed that both direct and indirect thickness measurement methods are reliable options for measuring thickness during compressibility measurements, provided that the potential sources of variability are carefully considered. Overall, the compliance measurements taken with laser sensors exhibited a larger variability in comparison with those recorded by the UTM. This discrepancy can be attributed to the sensors' higher sensitivity to minor changes in the plate's orientation, as well as their susceptibility to external interferences or vibrations from the compaction plate's movement.

The compression curves measured by the laser sensors showed a shifted curve changing on a sample-to-sample basis suggesting that the compaction plates have a certain level of mobility, resulting in more displacement when the glass plate is centrally loaded with fabrics. The parallelism between plates was checked with a pressure-sensitive paper, showing a nonuniform pressure distribution when a three-layer fabric stack was compressed. The surface finish of the compaction plate is likely causing a nonuniform pressure distribution. However, the obtained results are comparable to those reported in the first compressibility benchmark where the same fabric reinforcement was used.

Based on the findings reported in this paper, a list of practical recommendations and guidelines for calibration and accuracy assessment of both direct (laser sensors) and indirect (UTM) thickness measurement methods can be withdrawn:

* When using compaction plates with locking mechanisms, the compaction plate must be locked at a load at least 25% higher than the compression loads during testing to prevent any further readjustments during compaction.
<figure><img src="image_466.png" /><figcaption>Figure 7: Contact (A) between compaction plates, and (B) between three layers of Twill 2/2 (Hexcel fabric)</figcaption></figure>

An average thickness of 3.45±0.1 mm at a pressure of 0.1 MPa was reported in the first benchmark exercise under dry conditions, while in saturated conditions, the average thickness was 3.39±0.14 mm. In this study, at a pressure of 0.1 MPa a thickness of 3.31±0.01 mm and 3.34±0.02 mm was measured by the indirect and direct thickness methods, respectively. In saturated conditions, an average thickness of 3.33±0.06 mm was measured at this pressure. No significant differences were observed between the direct and indirect thickness measurement methods. The target thickness of 3 mm was reached at a compaction pressure of 393±39 kPa under dry conditions, while the reported participant's average is 370±141 kPa. Under saturated conditions, the target thickness was achieved at a pressure of 360±41 kPa, which is also comparable to the reported pressure of 317±127 kPa.

## Summary

This study showed that both direct and indirect thickness measurement methods are reliable options for measuring thickness during compressibility measurements, provided that the potential sources of variability are carefully considered. Overall, the compliance measurements taken with laser sensors exhibited a larger variability in comparison with those recorded by the UTM. This discrepancy can be attributed to the sensors' higher sensitivity to minor changes in the plate's orientation, as well as their susceptibility to external interferences or vibrations from the compaction plate's movement.

The compression curves measured by the laser sensors showed a shifted curve changing on a sample-to-sample basis suggesting that the compaction plates have a certain level of mobility, resulting in more displacement when the glass plate is centrally loaded with fabrics. The parallelism between plates was checked with a pressure-sensitive paper, showing a nonuniform pressure distribution when a three-layer fabric stack was compressed. The surface finish of the compaction plate is likely causing a nonuniform pressure distribution. However, the obtained results are comparable to those reported in the first compressibility benchmark where the same fabric reinforcement was used.

Based on the findings reported in this paper, a list of practical recommendations and guidelines for calibration and accuracy assessment of both direct (laser sensors) and indirect (UTM) thickness measurement methods can be withdrawn:

* When using compaction plates with locking mechanisms, the compaction plate must be locked at a load at least 25% higher than the compression loads during testing to prevent any further readjustments during compaction.
* To eliminate any non-sample displacement that may result from squeezing the self-aligning ball joint, it is advisable to compress the rig in a cyclic manner up to the same load as in the compliance curves after locking the compaction plate.

* Register the machine compliance up to the locking pressure using both distance sensors and the UTM. This will check for compliance differences caused by locking the compaction plate at a higher load than the load encountered during testing.

* To achieve accurate results and set the achievable accuracy, frequent compliance measurements should be taken to monitor changes and minimize their impact in both direct and indirect methods. Compliance curves can easily change due to common rig operations becoming important to monitor compliance changes regularly.

## Acknowledgements

The authors express their sincere gratitude to Dr. Ana Yong from NPL for the organization, support, and coordination of the second international benchmark. This work is part of a PhD funded by project PERMEA C2/16/24 accredited by the KU Leuven Research Council.

## References

[1] A. Endruweit and A. C. Long, "Analysis of compressibility and permeability of selected 3D woven reinforcements," J Compos Mater, vol. 44, no. 24, pp. 2833-2862, 2010. https://doi.org/10.1177/0021998310369586

[2] A. X. H. Yong et al., "Experimental characterisation of textile compaction response: A benchmark exercise," Compos Part A Appl Sci Manuf, vol. 142, 2021. https://doi.org/10.1016/j.compositesa.2020.106243

[3] D. May et al., "In-plane permeability characterization of engineering textiles based on radial flow experiments: A benchmark exercise," Compos Part A Appl Sci Manuf, vol. 121, pp. 100-114, 2019. https://doi.org/10.1016/j.compositesa.2019.03.006

[4] Y. Ana, M. David, and E. Andreas, "Measurement of textile compaction response and out-of-plane permeability: 2nd international benchmarking exercise," Proceedings of the 20th European Conference on Composite Materials - Composites Meet Sustainability (Vol 1-6). Lausanne, 26-30 June 2022.

[5] P. Sousa, S. V. Lomov, and J. Ivens, "Hurdles and limitations for design of a radial permeameter conforming to the benchmark requirements," Front Mater, vol. 9, Sep. 2022. https://doi.org/10.3389/fmats.2022.871235

[6] Z. Yousaf, P. Potluri, and P. J. Withers, "Influence of Tow Architecture on Compaction and Nesting in Textile Preforms," Applied Composite Materials, vol. 24, no. 2, 2017. https://doi.org/10.1007/s10443-016-9554-8
* To eliminate any non-sample displacement that may result from squeezing the self-aligning ball joint, it is advisable to compress the rig in a cyclic manner up to the same load as in the compliance curves after locking the compaction plate.

* Register the machine compliance up to the locking pressure using both distance sensors and the UTM. This will check for compliance differences caused by locking the compaction plate at a higher load than the load encountered during testing.

* To achieve accurate results and set the achievable accuracy, frequent compliance measurements should be taken to monitor changes and minimize their impact in both direct and indirect methods. Compliance curves can easily change due to common rig operations becoming important to monitor compliance changes regularly.

## Acknowledgements

The authors express their sincere gratitude to Dr. Ana Yong from NPL for the organization, support, and coordination of the second international benchmark. This work is part of a PhD funded by project PERMEA C2/16/24 accredited by the KU Leuven Research Council.

## References

[1] A. Endruweit and A. C. Long, "Analysis of compressibility and permeability of selected 3D woven reinforcements," J Compos Mater, vol. 44, no. 24, pp. 2833-2862, 2010. https://doi.org/10.1177/0021998310369586

[2] A. X. H. Yong et al., "Experimental characterisation of textile compaction response: A benchmark exercise," Compos Part A Appl Sci Manuf, vol. 142, 2021. https://doi.org/10.1016/j.compositesa.2020.106243

[3] D. May et al., "In-plane permeability characterization of engineering textiles based on radial flow experiments: A benchmark exercise," Compos Part A Appl Sci Manuf, vol. 121, pp. 100-114, 2019. https://doi.org/10.1016/j.compositesa.2019.03.006

[4] Y. Ana, M. David, and E. Andreas, "Measurement of textile compaction response and out-of-plane permeability: 2nd international benchmarking exercise," Proceedings of the 20th European Conference on Composite Materials - Composites Meet Sustainability (Vol 1-6). Lausanne, 26-30 June 2022.

[5] P. Sousa, S. V. Lomov, and J. Ivens, "Hurdles and limitations for design of a radial permeameter conforming to the benchmark requirements," Front Mater, vol. 9, Sep. 2022. https://doi.org/10.3389/fmats.2022.871235

[6] Z. Yousaf, P. Potluri, and P. J. Withers, "Influence of Tow Architecture on Compaction and Nesting in Textile Preforms," Applied Composite Materials, vol. 24, no. 2, 2017. https://doi.org/10.1007/s10443-016-9554-8
# Forming of thermoplastic polymer and magnesium alloy-based fiber metal laminates at elevated temperatures

LIU Zheng<sup>1,a*</sup>, SIMONETTO Enrico<sup>1,b</sup>, GHIOTTI Andrea<sup>1,c</sup>, BRUSCHI Stefania<sup>1,d</sup>

<sup>1</sup>Department of Industrial Engineering, University of Padova, Via Venezia 1, 35131 Padova, Italy

<sup>a</sup>zheng.liu@studenti.unipd.it, <sup>b</sup>enrico.simonetto.1@unipd.it, <sup>c</sup>andrea.ghiotti@unipd.it,

<sup>d</sup>stefania.bruschi@unipd.it

**Keywords:** Fiber Metal Laminates, Magnesium Alloy, Forming

**Abstract.** This paper illustrates the thermoforming process carried out on thermoplastic polymer and magnesium-based fiber metal laminates (FMLs). Flat laminates were formed at elevated temperatures into a hat-shape part. The forming force was acquired and, after forming, the thickness of each constituent of the FMLs in different zones was measured. A non-uniform thickness distribution was found in the formed parts, with a significant reduction of the prepreg thickness at the part bottom radii. Moreover, it was observed that the higher the blank-holder force the higher the forming force and the more significant the prepreg thickness variation.

## Introduction

Currently, the aerospace and automotive sectors are encountering increasing environmental and economic challenges in the production of parts that must satisfy high stringent quality standards. Specifically, improving energy efficiency in transportation can be achieved by reducing the weight of structural components. Fiber Metal Laminates (FMLs), which consist of a combination of metal and fiber-reinforced polymer (FRP) sheets also known as prepreg, exhibit significant potential to this regard. This is due to their advantageous strength-to-weight ratio, as well as their favorable performance in terms of fatigue and impact resistance [1]. In particular, hat-shaped parts (also known as hat channels or hat sections) can be used in both aerospace and automotive industries for structural and design applications such as hat-stringer-stiffened panels on the aircraft fuselage [2] and hat-channel components in vehicles.

The process for manufacturing FML components typically entails the individual shaping of the metal and prepreg sheets according to the desired configuration, subsequently followed by the adhesive bonding of the metal and prepreg components [3,4]. As alternative, thermoforming can be applied to form the FML as a whole at elevated temperature and pressure. The applied thermal field softens the polymer matrix, allowing its deformation and bonding with the metal layers, whereas the applied pressure contributes to ensure a good adhesion between the different layers. This method is often used for manufacturing FML parts with complex shapes or curved surfaces. Moreover, as a one-step forming technique, it facilitates the manufacturing process by offering a simultaneous forming and bonding process.

The polymer matrices utilized in FMLs can be classified as either thermosetting or thermoplastic. Thermosetting polymers, such as epoxy resin, have found extensive applications in the field: the so-called GLARE (glass-reinforced aluminum) is an example of FMLs currently used in the aerospace sector [5]. However, one of the drawbacks of using thermosetting polymers is given by the long duration of the forming cycle that must include the curing process of the polymer [6]. In contrast, thermoplastic polymers offer a viable alternative solution as they are characterized by excellent recyclability, remarkable toughness, and reduced forming time [7]. The thermoforming process employed for thermoplastic polymer-based FMLs necessitates to be carried out at temperatures above the thermoplastic resin's melting point. Under such conditions,

<figure><img src="image_456.png" /></figure>
# Forming of thermoplastic polymer and magnesium alloy-based fiber metal laminates at elevated temperatures

LIU Zheng<sup>1,a*</sup>, SIMONETTO Enrico<sup>1,b</sup>, GHIOTTI Andrea<sup>1,c</sup>, BRUSCHI Stefania<sup>1,d</sup>

<sup>1</sup>Department of Industrial Engineering, University of Padova, Via Venezia 1, 35131 Padova, Italy

<sup>a</sup>zheng.liu@studenti.unipd.it, <sup>b</sup>enrico.simonetto.1@unipd.it, <sup>c</sup>andrea.ghiotti@unipd.it,

<sup>d</sup>stefania.bruschi@unipd.it

**Keywords:** Fiber Metal Laminates, Magnesium Alloy, Forming

**Abstract.** This paper illustrates the thermoforming process carried out on thermoplastic polymer and magnesium-based fiber metal laminates (FMLs). Flat laminates were formed at elevated temperatures into a hat-shape part. The forming force was acquired and, after forming, the thickness of each constituent of the FMLs in different zones was measured. A non-uniform thickness distribution was found in the formed parts, with a significant reduction of the prepreg thickness at the part bottom radii. Moreover, it was observed that the higher the blank-holder force the higher the forming force and the more significant the prepreg thickness variation.

## Introduction

Currently, the aerospace and automotive sectors are encountering increasing environmental and economic challenges in the production of parts that must satisfy high stringent quality standards. Specifically, improving energy efficiency in transportation can be achieved by reducing the weight of structural components. Fiber Metal Laminates (FMLs), which consist of a combination of metal and fiber-reinforced polymer (FRP) sheets also known as prepreg, exhibit significant potential to this regard. This is due to their advantageous strength-to-weight ratio, as well as their favorable performance in terms of fatigue and impact resistance [1]. In particular, hat-shaped parts (also known as hat channels or hat sections) can be used in both aerospace and automotive industries for structural and design applications such as hat-stringer-stiffened panels on the aircraft fuselage [2] and hat-channel components in vehicles.

The process for manufacturing FML components typically entails the individual shaping of the metal and prepreg sheets according to the desired configuration, subsequently followed by the adhesive bonding of the metal and prepreg components [3,4]. As alternative, thermoforming can be applied to form the FML as a whole at elevated temperature and pressure. The applied thermal field softens the polymer matrix, allowing its deformation and bonding with the metal layers, whereas the applied pressure contributes to ensure a good adhesion between the different layers. This method is often used for manufacturing FML parts with complex shapes or curved surfaces. Moreover, as a one-step forming technique, it facilitates the manufacturing process by offering a simultaneous forming and bonding process.

The polymer matrices utilized in FMLs can be classified as either thermosetting or thermoplastic. Thermosetting polymers, such as epoxy resin, have found extensive applications in the field: the so-called GLARE (glass-reinforced aluminum) is an example of FMLs currently used in the aerospace sector [5]. However, one of the drawbacks of using thermosetting polymers is given by the long duration of the forming cycle that must include the curing process of the polymer [6]. In contrast, thermoplastic polymers offer a viable alternative solution as they are characterized by excellent recyclability, remarkable toughness, and reduced forming time [7]. The thermoforming process employed for thermoplastic polymer-based FMLs necessitates to be carried out at temperatures above the thermoplastic resin's melting point. Under such conditions,

<figure><img src="image_467.png" /></figure>
the prepreg can undergo dynamic thickness variations at varying process parameters due to the flow characteristics of the polymer matrix at elevated temperature.

It is important to emphasize that the FML, being a hybrid system, faces difficulties in quality control due to the heterogeneity in the metal and prepreg characteristics, both during the manufacturing process and throughout the service life of the formed component [8]. Hence, it is necessary to acquire a more extensive comprehension of the deformation mechanisms exhibited by the metal and prepreg sheets during the forming process. One of the deformation mechanisms that occurs during the forming process is the transverse flow of the prepreg, particularly when dealing with parts that have a curved shape and high thickness [9]. During the deformation process, the transverse flow of the prepreg can occur due to an uneven distribution of the normal pressure at the bending regions. This leads to the transverse squeezing of the prepreg. Consequently, a modification in the volume fraction of fibers occurs at the same time [10]. The investigation conducted in [11] examined the compaction behavior of cross-plied unidirectional (UD) carbon fiber reinforced thermoplastic polymers under varying temperature and normal pressure. At temperatures exceeding their melting points, thermoplastic polymers demonstrate significant deformations.

The variation in the part thickness during the FML forming process is also heavily influenced by the geometrical characteristics of the part. When stamping channel parts, it was observed that the thickness at the bottom radii decreased due to the high compaction forces in these areas [12]. The findings in [13] also demonstrated comparable outcomes when stamping cup-shaped components, where material accumulation was observed in the bottom and side walls adjacent to the radii areas.

To gain a more comprehensive understanding of how the magnesium alloy-based FMLs with a PA6 thermoplastic polymer matrix behave when they are thermoformed into hat-shaped parts, thermoforming tests were conducted at varying process parameters and the related thickness distribution of the different FML constituents was evaluated.

## Materials and methods

Materials. The FMLs investigated in this study consist of two layers of AZ31B magnesium alloy sheets (thickness 0.5 mm) as skins and one layer of glass fibre-reinforced PA6 consolidated prepreg (Tepex® 102-RG600 (2)/47% Type B, thickness 1 mm) as core. The mechanical characteristics of the magnesium alloy sheets and prepregs in the as-received condition were determined by carrying out standard tensile tests (ISO-6892) on a 50 kN MTS™ 322 hydraulic dynamometer, see Table 1. The main data of the prepregs are listed in Table 2.

<table><caption>Table 1. Mechanical properties of the AZ31B sheets and prepregs.</caption><thead><tr><th colspan="2">AZ31B sheets</th><th colspan="2">Prepregs</th></tr></thead><tbody><tr><td>Elastic Modulus</td><td>45 GPa</td><td>Tensile modulus</td><td>18 GPa</td></tr><tr><td>Yield Strength</td><td>158±2 MPa</td><td>Tensile strength</td><td>380±5 MPa</td></tr><tr><td>Tensile Strength</td><td>248±4 MPa</td><td>Strain at break</td><td>2.3%</td></tr><tr><td>Shear Modulus</td><td>16.7 GPa</td><td>Flexural Modulus</td><td>16 GPa</td></tr><tr><td>Poisson's Ratio</td><td>0.35</td><td>Flexural Strength</td><td>300 ±5 MPa</td></tr></tbody></table>
the prepreg can undergo dynamic thickness variations at varying process parameters due to the flow characteristics of the polymer matrix at elevated temperature.

It is important to emphasize that the FML, being a hybrid system, faces difficulties in quality control due to the heterogeneity in the metal and prepreg characteristics, both during the manufacturing process and throughout the service life of the formed component [8]. Hence, it is necessary to acquire a more extensive comprehension of the deformation mechanisms exhibited by the metal and prepreg sheets during the forming process. One of the deformation mechanisms that occurs during the forming process is the transverse flow of the prepreg, particularly when dealing with parts that have a curved shape and high thickness [9]. During the deformation process, the transverse flow of the prepreg can occur due to an uneven distribution of the normal pressure at the bending regions. This leads to the transverse squeezing of the prepreg. Consequently, a modification in the volume fraction of fibers occurs at the same time [10]. The investigation conducted in [11] examined the compaction behavior of cross-plied unidirectional (UD) carbon fiber reinforced thermoplastic polymers under varying temperature and normal pressure. At temperatures exceeding their melting points, thermoplastic polymers demonstrate significant deformations.

The variation in the part thickness during the FML forming process is also heavily influenced by the geometrical characteristics of the part. When stamping channel parts, it was observed that the thickness at the bottom radii decreased due to the high compaction forces in these areas [12]. The findings in [13] also demonstrated comparable outcomes when stamping cup-shaped components, where material accumulation was observed in the bottom and side walls adjacent to the radii areas.

To gain a more comprehensive understanding of how the magnesium alloy-based FMLs with a PA6 thermoplastic polymer matrix behave when they are thermoformed into hat-shaped parts, thermoforming tests were conducted at varying process parameters and the related thickness distribution of the different FML constituents was evaluated.

## Materials and methods

Materials. The FMLs investigated in this study consist of two layers of AZ31B magnesium alloy sheets (thickness 0.5 mm) as skins and one layer of glass fibre-reinforced PA6 consolidated prepreg (Tepex® 102-RG600 (2)/47% Type B, thickness 1 mm) as core. The mechanical characteristics of the magnesium alloy sheets and prepregs in the as-received condition were determined by carrying out standard tensile tests (ISO-6892) on a 50 kN MTS™ 322 hydraulic dynamometer, see Table 1. The main data of the prepregs are listed in Table 2.

<table><caption>Table 1. Mechanical properties of the AZ31B sheets and prepregs.</caption><thead><tr><th colspan="2">AZ31B sheets</th><th colspan="2">Prepregs</th></tr></thead><tbody><tr><td>Elastic Modulus</td><td>45 GPa</td><td>Tensile modulus</td><td>18 GPa</td></tr><tr><td>Yield Strength</td><td>158±2 MPa</td><td>Tensile strength</td><td>380±5 MPa</td></tr><tr><td>Tensile Strength</td><td>248±4 MPa</td><td>Strain at break</td><td>2.3%</td></tr><tr><td>Shear Modulus</td><td>16.7 GPa</td><td>Flexural Modulus</td><td>16 GPa</td></tr><tr><td>Poisson's Ratio</td><td>0.35</td><td>Flexural Strength</td><td>300 ±5 MPa</td></tr></tbody></table>
<table><caption>Table 2. Main characteristics of the prepregs.</caption><thead><tr><th rowspan="2">Layup</th><th colspan="2">Value</th><th rowspan="2">Unit</th></tr><tr><th>Longitudinal</th><th>Transversal</th></tr></thead><tbody><tr><td>Fibre</td><td colspan="2">E-Glass</td><td></td></tr><tr><td>Weaving style</td><td colspan="2">Twill 2/2</td><td></td></tr><tr><td>Area weight (dry)</td><td colspan="2">600</td><td>g/m<sup>2</sup></td></tr><tr><td>Yarn count</td><td colspan="2">1200</td><td>tex</td></tr><tr><td>Yarn density</td><td>2.5</td><td>2.5</td><td>1/cm</td></tr><tr><td>Weight rate</td><td>50</td><td>50</td><td>%</td></tr><tr><td>Fibre content</td><td colspan="2">47</td><td>Vol-%</td></tr><tr><td>Thickness per layer</td><td colspan="2">0.5</td><td>mm</td></tr><tr><td>Laminate density</td><td colspan="2">1.8</td><td>g/cm<sup>3</sup></td></tr></tbody></table>

Thermoforming tests. To assess the possible thickness variation of the FML when it is thermoformed into a hat shape part, thermoforming tests were conducted by employing the testing device and FML specimen shown in Fig. 1 (a) and (b), which provides the configuration before and after processing. FML specimens with the dimension of 80x160 mm were positioned as shown in Fig.1 (a). They were heated in an MTS™ 651 environmental chamber to the designed temperature before being tested on the MTS™ 322 hydraulic dynamometer. The FML sheet was stamped by using a punch speed of 1 mm/s with the holding force exerted by the left and right blank holders, while the lower die was fixed on the MTS frame. The FML specimens were formed at 255°C (35°C higher than the melting point of PA6) under different holding forces. Before assembling the specimen, the testing set-up was preheated to the designed temperature for 1 hour to ensure an even temperature distribution and avoid the difference in the temperature between the set-up and specimen. After the stamping process, the forming and holding force were kept for 240 s, followed by cooling down the specimens and dies to room temperature with still the force applied to finish the consolidation process.

<figure><img src="image_457.png" /><figcaption>Fig. 1 Thermoforming set-up for hat-shaped parts: before (a) and after (b) forming.</figcaption></figure>
<table><caption>Table 2. Main characteristics of the prepregs.</caption><thead><tr><th rowspan="2">Layup</th><th colspan="2">Value</th><th rowspan="2">Unit</th></tr><tr><th>Longitudinal</th><th>Transversal</th></tr></thead><tbody><tr><td>Fibre</td><td colspan="2">E-Glass</td><td></td></tr><tr><td>Weaving style</td><td colspan="2">Twill 2/2</td><td></td></tr><tr><td>Area weight (dry)</td><td colspan="2">600</td><td>g/m²</td></tr><tr><td>Yarn count</td><td colspan="2">1200</td><td>tex</td></tr><tr><td>Yarn density</td><td>2.5</td><td>2.5</td><td>1/cm</td></tr><tr><td>Weight rate</td><td>50</td><td>50</td><td>%</td></tr><tr><td>Fibre content</td><td colspan="2">47</td><td>Vol-%</td></tr><tr><td>Thickness per layer</td><td colspan="2">0.5</td><td>mm</td></tr><tr><td>Laminate density</td><td colspan="2">1.8</td><td>g/cm³</td></tr></tbody></table>

Thermoforming tests. To assess the possible thickness variation of the FML when it is thermoformed into a hat shape part, thermoforming tests were conducted by employing the testing device and FML specimen shown in Fig. 1 (a) and (b), which provides the configuration before and after processing. FML specimens with the dimension of 80x160 mm were positioned as shown in Fig.1 (a). They were heated in an MTS™ 651 environmental chamber to the designed temperature before being tested on the MTS™ 322 hydraulic dynamometer. The FML sheet was stamped by using a punch speed of 1 mm/s with the holding force exerted by the left and right blank holders, while the lower die was fixed on the MTS frame. The FML specimens were formed at 255°C (35°C higher than the melting point of PA6) under different holding forces. Before assembling the specimen, the testing set-up was preheated to the designed temperature for 1 hour to ensure an even temperature distribution and avoid the difference in the temperature between the set-up and specimen. After the stamping process, the forming and holding force were kept for 240 s, followed by cooling down the specimens and dies to room temperature with still the force applied to finish the consolidation process.

<figure><img src="image_468.png" /><figcaption>Fig. 1 Thermoforming set-up for hat-shaped parts: before (a) and after (b) forming.</figcaption></figure>
Two levels of holding force, namely 500 N ($Hf_1$) and 2000 N($Hf_2$), were selected to assess the effect of the blank-holder force on the FML thickness distribution. For each testing condition, three tests were performed to ensure the repeatability of the results. After the tests the samples were cut in A-A plane as shown in Fig. 2 and the through-thickness of each layer of the FML at different areas was measured by employing the Dino-Lite digital microscope. The hat-shaped part can be divided into five zones (from Z1 to Z5) according to its geometrical characteristics and the difference in loading conditions during the forming process.

<figure><img src="image_458.png" /><figcaption>Fig. 2 Different zones of the hat-shaped FML for thickness measurement</figcaption></figure>

## Result and discussion

Thickness distribution. Fig. 3 presents the average thickness distribution of the two AZ31B layers. Overall, the thickness remains constant around 0.5 mm in each zone, regardless of the blank-holding force, which means no significant thickening and thinning occur on the AZ31B sheets during the thermoforming process. This can be explained that the thickness of the metal sheets during the forming process is mainly governed by the in-plain stretching deformation. The low holding force in this research allows the relative slide between metal sheets and dies, resulting in no significant stretching of the AZ31B sheets. The thickness of both the AZ31B layers at Z2 and Z4 zones, where bending occurs, is around 5% lower than in the other zones. Moreover, the larger bending radius of layer 1 compared to layer 2 in Z2 zone contributes to the lower thickness of AZ31B at layer 1. Similarly, a thicker value can be observed at layer 2 compared to that of layer 1 as a result of the lower bending radius of layer 2 on Z4 zone.

<figure><img src="image_459.png" /><figcaption>Fig. 3 AZ31B thickness distribution at upper layer (a) and lower layer (b).</figcaption></figure>

As expected, non-uniform thickness can be observed in the different zones of the prepreg (see Fig. 4). The higher the interfacial pressure the lower the thickness of the prepreg due to the high mobility and flow of the prepreg when the temperature is above the melting point of the polymer matrix, resulting in the thickness of the prepreg layer governed by the interfacial contact pressure. The blank-holding force applied at Z1 zone and the contact between the AZ31B layers and dies at Z4 and Z5 zones both contribute to the higher interfacial pressure in these zones compared to that in Z2 and Z3 zones: this contributes to the prepreg being squeezed from Z1, Z4 and Z5 zones to
Two levels of holding force, namely 500 N ($Hf_1$) and 2000 N($Hf_2$), were selected to assess the effect of the blank-holder force on the FML thickness distribution. For each testing condition, three tests were performed to ensure the repeatability of the results. After the tests the samples were cut in A-A plane as shown in Fig. 2 and the through-thickness of each layer of the FML at different areas was measured by employing the Dino-Lite digital microscope. The hat-shaped part can be divided into five zones (from Z1 to Z5) according to its geometrical characteristics and the difference in loading conditions during the forming process.

<figure><img src="image_469.png" /><figcaption>Fig. 2 Different zones of the hat-shaped FML for thickness measurement</figcaption></figure>

## Result and discussion

Thickness distribution. Fig. 3 presents the average thickness distribution of the two AZ31B layers. Overall, the thickness remains constant around 0.5 mm in each zone, regardless of the blank-holding force, which means no significant thickening and thinning occur on the AZ31B sheets during the thermoforming process. This can be explained that the thickness of the metal sheets during the forming process is mainly governed by the in-plain stretching deformation. The low holding force in this research allows the relative slide between metal sheets and dies, resulting in no significant stretching of the AZ31B sheets. The thickness of both the AZ31B layers at Z2 and Z4 zones, where bending occurs, is around 5% lower than in the other zones. Moreover, the larger bending radius of layer 1 compared to layer 2 in Z2 zone contributes to the lower thickness of AZ31B at layer 1. Similarly, a thicker value can be observed at layer 2 compared to that of layer 1 as a result of the lower bending radius of layer 2 on Z4 zone.

<figure><img src="image_470.png" /><figcaption>Fig. 3 AZ31B thickness distribution at upper layer (a) and lower layer (b).</figcaption></figure>

As expected, non-uniform thickness can be observed in the different zones of the prepreg (see Fig. 4). The higher the interfacial pressure the lower the thickness of the prepreg due to the high mobility and flow of the prepreg when the temperature is above the melting point of the polymer matrix, resulting in the thickness of the prepreg layer governed by the interfacial contact pressure. The blank-holding force applied at Z1 zone and the contact between the AZ31B layers and dies at Z4 and Z5 zones both contribute to the higher interfacial pressure in these zones compared to that in Z2 and Z3 zones: this contributes to the prepreg being squeezed from Z1, Z4 and Z5 zones to