# Multiscale modeling of mechanically recycled glass fiber reinforced polyamide 6 composites accounting for viscoelasticity, viscoplasticity, and anisotropic damage

Saif Eddine Sekkal, Fodil Meraghni, George Chatzigeorgiou, Francis Praud, Nelly Durand

► To cite this version:

Saif Eddine Sekkal, Fodil Meraghni, George Chatzigeorgiou, Francis Praud, Nelly Durand. Multi-scale modeling of mechanically recycled glass fiber reinforced polyamide 6 composites accounting for viscoelasticity, viscoplasticity, and anisotropic damage. Composite Structures, 2025, 361, pp.119016. 10.1016/j.compstruct.2025.119016. hal-05008164

HAL Id: hal-05008164

https://artsetmetiers.hal.science/hal-05008164v1

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Multiscale modeling of mechanically recycled glass fiber reinforced polyamide 6 composites accounting for viscoelasticity, viscoplasticity, and anisotropic damage

S.E. Sekkal<sup>a,b</sup>, F. Meraghni<sup>a ID,*</sup>, G. Chatzigeorgiou<sup>a ID</sup>, F. Praud<sup>a ID</sup>, N. Durand<sup>b</sup>

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

## ARTICLE INFO

### Keywords:

Mechanically recycled composites

Multiscale nonlinear modeling

Microstructure generation

Compression molding

## ABSTRACT

Fiber-reinforced thermoplastic composites are valued for their strength-to-weight ratio, cost-effectiveness, and recyclability, highlighting the need for efficient recycling technologies amid environmental concerns. This study addresses these challenges by examining the mechanical response of recycled glass fiber reinforced polyamide 6 composites and modeling their nonlinear, time-dependent behavior under complex loading conditions. Advanced nonlinear constitutive and multiscale models, initially developed for conventional fiber composites, are adapted to capture the stochastic response of recycled materials. These models integrate viscoelasticity, viscoplasticity and damage in the polymer matrix and account for anisotropic damage in the strands, addressing the heterogeneity introduced by the recycling process. A modified random sequential adsorption technique replicates the microstructures for nonlinear response modeling. Hypotheses based on microstructural investigations consider processing effects that disrupt the initial chip woven structure and create matrix-rich areas. The model captures anisotropy and variability observed in experimental data, providing a reliable framework for predicting the performance of recycled thermoplastic composites and improving the understanding of the relationship between microstructure and mechanical properties, with a focus on inelastic nonlinear behavior.

## 1. Introduction

Continuous fiber reinforced thermoplastic composites have emerged as a robust solution in the automotive, aerospace, and renewable energy sectors. These materials offer superior mechanical properties, lower manufacturing costs, enhanced impact resistance, and excellent recyclability compared to thermoset counterparts [1–7]. As the range of applications and demand for these composites expand, research into their closed-loop recycling processes has become a pressing concern. Recent environmental directives aim to reduce waste sent to incineration or landfill, stressing the importance of developing sustainable, viable, and economically attractive recycling technologies [8].

Composite recycling technologies are classified into mechanical, thermal, and chemical methods, chosen based on the material type and intended product application [5,6]. Thermal recycling of carbon fibers involves pyrolysis, breaking down polymer matrices in an inert atmosphere above 350 °C, ideally suited for thermoset composites to recover fibers alone. However, Mazzocchetti et al. noted recycled carbon fiber composites generally have lower mechanical properties compared

to virgin carbon fiber composites, with inconsistent crystallinity and impurities observed on the fibers [8,9]. Chemical recycling via solvolysis breaks down polymer matrices using reactive solvents like water or organic solvents at varying temperatures and pressures, achieving clean fiber extraction with minimal mechanical property loss and maintaining fiber length. However, this method is costly due to specialized high-temperature, high-pressure reactors needed, thus lagging behind thermal and mechanical methods in industrial implementation [10].

In contrast, thermoplastic matrices, owing to their meltability compared to non-meltable thermoset polymers, facilitate easier handling in recycling processes, making them particularly suitable for mechanical recycling by shredding or cutting, followed by thermocompression or injection molding processes [1,7,11–14]. Colucci et al. [12] studied mechanical recycling effects on polymeric composite behavior and microstructure via injection molding post artificial aging. Their results indicated minimal impact of recycling and aging on final properties and microstructure. However, the injection molding process, which is shear-intensive, can lead to severe fiber attrition, resulting in short fiber

<footnote><sup>a</sup> Arts et Métiers Institute of Technology, CNRS, Université de Lorraine, LEM3-UMR7239, 4 Rue Augustin Fresnel, Metz, 57078, France</footnote>

<footnote><sup>b</sup> Cetim, 21 Rue de Chemnitz, Mulhouse, 68200, France</footnote>

<footnote>* Corresponding author.
E-mail address: fodil.meraghni@ensam.eu (F. Meraghni).</footnote>
composites with low impact strength [13,14]. Processing parameters such as mold temperature and extruder speed play a crucial role in determining the mechanical properties of recycled composites.

An alternative mechanical recycling route for thermoplastic composites involves direct processing without prior material homogenization. This method includes chopping the material into centimeter-sized chips and using compression molding to form simple or complex parts. Tapper et al. [15] assessed closed-loop recycling of discontinuous carbon fiber polyamide 6 composites, noting a significant decrease in tensile stiffness and strength over two recycling cycles due to fiber misalignment and breakage under high compaction pressures. Other studies have explored similar techniques using long chopped unidirectional carbon fiber thermoplastic tapes [7,16–18]. In contrast, Kiss et al. [4] proposed a variation of this approach by co-molding shredded material as a core layer between continuous fiber skins. These sandwich-like panels demonstrated significant potential for upcycling, restoring original composite laminate properties in flexural and impact loading. Control of processing parameters such as mold temperature and extruder speed remains critical in determining the mechanical properties of these recycled composites.

The literature on modeling nonlinear mechanical behavior in recycled composite materials is scarce. Modeling the mechanical properties of these materials requires a deep understanding of the complex microstructures induced by the recycling process. The challenge lies in accurately identifying the microstructural parameters that have the most influence on the performance of recycled composites, while effectively representing them in the model. In contrast, research on prepreg platelet molded composites and randomly oriented strand composites provides insights applicable to recycled thermoplastic composites. Kravchenko et al. conducted computational analyses on prepreg platelet molded composites (PPMCs), focusing on stochastic meso-scale morphology and its impact on tensile properties under uniaxial tension [19–21]. Their work utilized progressive failure analysis to simulate distributions of tensile properties, emphasizing how factors such as platelet geometry and system morphology contribute to variations in effective tensile properties. The study validated predictions using experimental data, highlighting the influence of flow-induced fiber orientation and compaction deformations on composite properties. Harper et al. [22] developed a geometrical modeling scheme to create realistic random fiber architectures for discontinuous fiber composites, facilitating subsequent modeling of their mechanical properties. Tang et al. [23] investigated the fatigue behavior of chopped carbon fiber chip-reinforced composites through experimental and numerical methods.

The structures of recycled composites are well-suited for homogenization and multi-scale approaches to predict their overall mechanical responses. Among the available theories, mean-field [24–27] and full-field approaches are particularly effective [28–30]. Several studies have used a combination of the two approaches for hierarchical composite materials, particularly for thermoplastic woven composites [31–33]. In recent years, Fast Fourier Transform (FFT)-based methods have been developed as an alternative to the Finite Element Method (FEM) for solving microstructural boundary value problems. Based on the foundational work of Moulinec and Suquet [34], these methods use a regular grid discretization to solve equilibrium equations in heterogeneous media [35,36].

The authors' previous work focused on detailed microstructural investigations of recycled thermoplastic composites to understand the relationship between microstructure and mechanical properties [37, 38]. The recycled material was initially modeled using a meso-macro scheme to represent nonlinear flexural behavior [37], and subsequently with a multi-scale predictive model that captured the observed anisotropy in the material [38]. This research laid the groundwork for detailing the microstructure generation process, capturing anisotropy through Orientation Distribution Function (ODF), and incorporating fiber structure degradation after recycling into the model. The current

work focuses on recycling glass fiber polyamide 6 composites. This matrix is highly sensitive to relative humidity, affecting the glass transition temperature, rheology, and damage mechanisms [39–42]. These polymers exhibit inelastic nonlinear behavior driven by viscoelastic and viscoplastic deformation mechanisms. Numerous constitutive models have been presented in the literature for these polymers [43–46].

The novelty of this work lies in the following key contributions:

* Advanced nonlinear constitutive and multiscale models are adapted to address the unique, heterogeneous microstructure of recycled composites. These models, initially designed for conventional fiber composites (woven and short fiber), incorporate viscoelasticity and viscoplasticity with damage in the polymer matrix, as well as anisotropic damage in the strands. To address the challenges of the microstructure complexity of a recycled material, a modified random sequential adsorption technique has been integrated into the multiscale framework.

* Simplifying hypotheses are formulated based on extensive microstructural investigations and material response examinations. These hypotheses account for the processing effects that disrupt the initial chip woven structure and create matrix-rich areas. By incorporating these experimental findings into the nonlinear homogenization framework, the model effectively captures the inherent heterogeneity and variability observed in experimental data.

The paper is structured as follows: Section 2 provides an overview of the examined material, the recycling process, and the resulting recycled composite microstructure. Section 3 outlines the theoretical background of periodic homogenization and the microstructure generation process for nonlinear analyses. Section 4 presents the constitutive models used to describe the matrix and strand behaviors. Finally, Section 5 discusses the experimental and numerical results concerning the nonlinear macroscopic response of recycled thermoplastic composites.

## 2. Recycled glass fiber reinforced polyamide 6 composites

### 2.1. *Recycling process*

The recycled composites studied in this work are manufactured and supplied by Cetim, a Research and Technology Organization (RTO) that contributed to the development of mechanical recycling of polymer composites. Their modular production line, “Thermosaïc®” recycles end-of-life parts and production offcuts into consolidated plates through compression molding. The process begins by shredding the material into small chips using a mill. The chips are then distributed evenly for proper consolidation.

The chips undergo two compression stages: the first stage preheats and pre-compacts the strands with a load not exceeding 50 kN and a temperature lower than the resin's fusion point. The second stage involves a higher temperature and load, not exceeding 1300 kN, to ensure proper matrix fusion without undesired squeeze flow. The formed plate is then cooled and ejected at the end of the production line. This process is summarized in Fig. 1.

To limit the number of factors and variables, this study focuses on recycled composites produced from uniformly cut chips sized 20 x 20 x 1 mm³. The original virgin material is a PA-GF woven composite with a fiber volume fraction of 47%, referenced as Tepex® Dynalite 102-RG600(x)/PA-GF/47%, and supplied by Bond-Laminates GmbH (Am Patbergschen Dorn 11-D-59929 Brilon, Germany). The resulting recycled plate is 5 mm thick. For the present work, the studied mechanically recycled composite will be referred to as RTPC (Recycled ThermoPlastic Composite).

The mechanical performance of these recycled composites is closely tied to the processing parameters used in their production. Selecting the optimal parameters for the best mechanical performance is challenging and can be expensive without numerical tools for modeling the process
<figure><img src="image_3.png" /><figcaption>Fig. 1. Overview of the Thermosaic® production line detailing the mechanical recycling process for thermoplastic composites [38].</figcaption></figure>

and mechanical performance. Specifically, modeling the mechanical behavior requires a deep understanding of the material microstructure and its variability with respect to the mechanical response. A comprehensive microstructural investigation, including an analysis of this variability, was conducted in previous work [38]. The main findings will be summarized in Section 2.2.

## 2.2. RTPC microstructure description and investigation

Optical microscopy observations reveals that the composite microstructure can be analyzed at two different scales: a microscopic scale and a mesoscopic scale. The main reason for this scale separation is the irregular distribution of the matrix in the material, as shown in Fig. 2. This heterogeneity is characterized by the presence of regions that are rich in matrix and others that are rich in fibers. At the mesoscopic scale, the fibers form a wave pattern with mostly in-plane orientation. It is possible to distinguish two different strands based on changes in orientation or the presence of matrix-rich regions as shown in Fig. 2. At the microscopic level, the microstructural morphology within the strands is more homogeneous, with multiple fibers having a similar orientation.

These microstructural observations provide crucial quantitative measurements for numerical simulations, offering detailed insights about the material's microstructure, particularly the volume fractions at relevant scales. At the microscale, the average fiber ratio inside the strand, $v_{f/s}$, is 75.36%. At the mesoscale, the strand ratio within the composite is computed assuming a global fiber volume fraction of 47%, obtained by matrix burn-off, resulting in a value of 62.36%. This data serves as input for the microstructure generation process preceding multiscale modeling.

X-ray tomography investigation of the studied material reveals that the woven structure is maintained at the plate center across its thickness, as depicted in Fig. 3(b). Notably, in several instances, the woven structure deviates from the initial 0° and 90° yarn orientations, likely due to material flow during thermocompression. Additionally, the upper and lower surfaces (Figs. 3(a) and 3(c)) exhibit fewer instances of woven structure, with multiple unidirectional fiber bundles apparent instead. This heterogeneity may result from uneven temperature and pressure distributions within the mold during processing.

## 3. Multi-scale modeling of recycled glass fiber reinforced polyamide 6 composites

A multi-scale model is generally based on a Representative Volume Element (RVE) describing the microstructural morphology of the composite over which its mechanical response can be objectively assessed.

To do so, the scale separation assumption states that the characteristic dimensions of the RVE are much smaller than the characteristic dimensions of the macroscopic medium. Basically, this implies that any macroscopic material point is associated with an RVE domain $\mathfrak{B}$ in which any microscopic material point is denoted by $\mathbf{x}$.

## 3.1. Background in full-field multi-scale modeling

As common in multi-scale modeling, the connection between the macroscopic and microscopic scales is established by averaging the stress and strain over the RVE domain $\mathfrak{B}$. Therefore, the macroscopic stress and strain are given by:

$$ \bar{\sigma} = \frac{1}{V} \int_{\mathfrak{B}} \sigma(\mathbf{x}) \, dV, \quad \bar{\epsilon} = \frac{1}{V} \int_{\mathfrak{B}} \epsilon(\mathbf{x}) \, dV, \qquad (1) $$

where $V$ denotes volume of the whole RVE domain $\mathfrak{B}$.

Assuming that the RVE can be represented as a repeating element, it can be idealized as a periodic unit cell in which the displacement field $\mathbf{u}$ takes the following additive form:

$$ \mathbf{u}(\mathbf{x}) = \bar{\epsilon} \cdot \mathbf{x} + \mathbf{u}'(\mathbf{x}) + \mathbf{u}_0, \quad \forall \mathbf{x} \in \mathfrak{B}, \qquad (2) $$

where the first term depicts the affine partition of the displacement field, while the second term, $\mathbf{u}'$ represents a periodic fluctuation that maintains identical values at each pair of opposite points, $\mathbf{x}_+$ and $\mathbf{x}_-$ within the RVE boundary $\partial\mathfrak{B}$:

$$ \mathbf{u}'(\mathbf{x}_+) = \mathbf{u}'(\mathbf{x}_-) \quad \forall \mathbf{x}_+, \mathbf{x}_- \in \partial\mathfrak{B}. \qquad (3) $$

Besides, it is worth noticing that, due to its periodicity, the part of the strain associated with $\mathbf{u}'$ vanishes once averaged, ensuring a full consistency between (1) and (2). Note that, in (2), $\mathbf{u}_0$ denotes an eventual rigid body motion.

Inserting (2) into (3), one can express the Periodic Boundary Conditions (PBCs) in terms of displacement $\mathbf{u}$, while directly involving the macroscopic strain $\bar{\epsilon}$. Those apply to the unit cell's boundaries $\partial\mathfrak{B}$ such that:

$$ \mathbf{u}(\mathbf{x}_+) - \mathbf{u}(\mathbf{x}_-) = \bar{\epsilon} \cdot (\mathbf{x}_+ - \mathbf{x}_-) \quad \forall \mathbf{x}_+, \mathbf{x}_- \in \partial\mathfrak{B} \qquad (4) $$

Based on this framework, the macroscopic stress-strain response of the composite results from the solution of the unit cell problem. In this work, the latter is solved through the finite element method, using periodically meshed unit cells as in what follows.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Optical light microscopy image of the RTPC microstructure, showing observed heterogeneities that lead to an analysis by scale separation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Comparison of the fiber structures at different relative out-of-plane coordinates $\delta^*$ for an RTPC specimen using X-ray micro-computed tomography. Micrographs show the upper surface (a), center (b), and lower surface (c) of the plate, highlighting the variability and heterogeneity in composite microstructure across the thickness.</figcaption></figure>

## 3.2. Microstructure generation for recycled glass fiber reinforced polyamide 6 composites

### 3.2.1. Microstructure generation process

The primary challenge in Recycled ThermoPlastic Composite RVE generation for this research is achieving the high strand volume fraction necessary to accurately represent the real microstructures of these materials. To attain this value, the RVE is constructed in multiple layers rather than by placing inclusions directly in the 3-D RVE. The microstructure generation algorithm is based on the Random Sequential Absorption (RSA) approach [47–50]. Strands are sequentially placed in each layer of the RVE according to a user-defined Orientation Distribution Function (ODF) until the required strand ratio is reached. The implementation, done in MATLAB, utilizes the **polyshape** and **intersect** functions to detect intersections between the generated strands.

A Python script is used to generate a finite element model using data from the Matlab generation code. The Matlab code supplies the coordinates of the strand vertices, the orientations of the strands, and the strand volume fraction during the generation process. The Python script then uses this data to generate a Finite Element (FE) mesh with the help of ABAQUS. This process includes positioning the strands with the vertex data, eliminating the unit cell edges, adding the matrix material, and assigning material orientations to each strand.

The generation process begins by randomly selecting strand center coordinates within the designated domain. To ensure periodicity within the RVE, any strand intersecting a boundary is duplicated on the opposite side. For strands positioned in a corner, duplication occurs on all other corners. Next, intersections between strands are checked. If an intersection is found, the iteration stops, and a counter increments to track the number of attempts. If the maximum number of attempts is
<figure><img src="image_6.png" /><figcaption>Fig. 4. (a) Generated unit cell showing the RVE. (b) RVE with the matrix hidden to better illustrate the spatial distribution of the strands through the thickness.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Illustration of the modified strand deposition process. (a) Strands are deposited in the first square until the maximum number of trials is achieved. (b) The process is repeated sequentially for each subsequent square, ensuring minimum spacing requirements are met.</figcaption></figure>

reached, the current RVE is discarded, and a new one is generated. This loop continues until the desired strand ratio is achieved. A flowchart detailing the generation process is presented in [38]. An example of a generated unit cell is shown in Fig. 4. Fig. 4(a) displays the RVE, while Fig. 4(b) shows the RVE with the matrix hidden.

For this work, a RVE unit cell is comprised of multiple layers. Each RVE, measuring 100×100×5 mm³, consists of four layers. Strands of 20×20 mm² are deposited until saturation ($v_{s/c}$ of 56%), with a condition of spacing between them to prevent high stress concentration. A MATLAB function computes the minimum euclidean distance between polygons, enforcing a minimum of 4 mm separation between strands [51]. However, this additional requirement significantly increased runtime. To mitigate this, the generation process was adjusted: the 100 × 100 mm area was divided into four squares. Strands are deposited in each square sequentially until saturation, ensuring compliance with spacing criteria. This process is illustrated in Fig. 5, where the generation process begins by depositing strands in the first square until a maximum amount of trials is achieved (Fig. 5(a)), before moving to the next square (Fig. 5(b)).

### 3.2.2. Discussion on the mechanical response of recycled composites and model fine-tuning

Previously performed tensile tests over different specimen orientations have shown that significant discrepancy in mechanical response was observed between the 0° (manufacturing direction) and the 90° orientations [38]. This can be attributed to the complex material flow inside the mold during thermocompression, that has likely oriented more fibers closer to the 0° direction, making it stiffer. The above leads to the conclusion that the RTPC material is not macroscopically transversely isotropic as was expected, but rather exhibits a certain

<figure><img src="image_8.png" /><figcaption>Fig. 6. Adopted orientation distribution function (ODF) for the strand generation.</figcaption></figure>

anisotropy. To capture the material's anisotropy, a specific Orientation Distribution Function (ODF) was applied during microstructure generation, as illustrated in Fig. 6. Given the complexity and heterogeneity of the material's microstructure, experimentally quantifying an ODF is challenging. Consequently, a numerical approach was employed, involving the generation of multiple RVEs with various ODFs to closely replicate the experimental observations [38].

Initially, the strands were modeled as a unidirectional composite. This simplification was based on various microstructural observations that revealed the material's complexity after processing: some chips
<figure><img src="image_9.png" /><figcaption>Fig. 7. Rheological scheme of the constitutive model for the thermoplastic matrix.</figcaption></figure>

retained a woven structure, while others displayed a unidirectional architecture. A mean-field Mori–Tanaka homogenization scheme calculates the strand material properties [52–54]. The computed elastic constants are then used for mesoscale full-field homogenization. However, this model’s implementation results in an overestimation of the mechanical response of RTPC material when compared to experimental data. This discrepancy may be due to the fibers sustaining damage during processing, which changes their orientations and reduces their properties. To correct this, a calibration process based on experimental results is implemented. The revised model introduces a correction coefficient, which accounts for the microstructural complexity of the strands and the degradation the chips undergo during processing. After calibration, the recovered Young’s and shear moduli values from micro-homogenization are reduced by 35%. Accounting for fiber degradation, Young’s moduli for randomly oriented strand RVEs approach experimental values. For more details about model fine-tuning to capture microstructural effects on material response, readers are referred to [38].

## 4. Constitutive models of the constituents

In the considered RVEs, two specific constitutive models are employed to describe the microscopic response of the matrix and strand phases within the unit cell. They are briefly presented in this section.

### 4.1. Constitutive model for the thermoplastic matrix

Thermoplastic matrices like polyamide 6 are known to exhibit a complex dissipative behavior combining both solids and fluids properties along with progressive degradation [55]. To describe this type of behavior, the phenomenological multi-mechanism constitutive model developed by [46,56] is utilized to represent the matrix phase in the present multi-scale modeling approach. Note that a brief overview of this constitutive model is given in this section for the sake of conciseness. However, readers are referred to above-mentioned references for a comprehensive description.

This constitutive model is based on a thermodynamically-consistent formulation integrating both viscoelastic and viscoplastic deformations as well as ductile damage, as depicted by the rheological scheme in Fig. 7. The latter consists in one single spring subjected to an elastic strain $\epsilon_e$, N viscoelastic Kelvin–Voigt branches, each one of them being subjected to a viscoelastic strain $\epsilon_{v_i}$, and a viscoplastic branch subjected to a plastic strain $\epsilon_p$ along with as hardening measure. All of these branches act under the effect of the effective stress $\tilde{\sigma}$ to represent damage which is regarded as an isotropic phenomenon through the damage variable $D$.

The constitutive equations of this model, which are presented briefly below, consider the total strain $\epsilon$ as the only external state variable, whereas $\epsilon_{v_i}$, $r$, $\epsilon_p$ and $D$ act as internal state variables.

$$ \sigma = (1 - D)C_e : \left( \epsilon - \sum_{i=1}^{N} \epsilon_{v_i} - \epsilon_p \right) \quad (5) $$

$$ -\sigma_{v_i} = (1 - D)C_{v_i} : \epsilon_{v_i} - \sigma \quad (6) $$

$$ R = Kr^n \quad (7) $$

$$ \dot{\epsilon}_{v_i} = \frac{V_{v_i}^{-1}}{1 - D} : \sigma_{v_i} \quad (8) $$

$$ \dot{\epsilon}_p = \frac{3 \mathbf{Dev}(\sigma)}{2 \mathbf{eq}(\sigma)} \frac{\dot{r}}{1 - D} \quad (9) $$

$$ \dot{D} = \Lambda(D) \left(\frac{Y}{S}\right)^{\beta} \frac{\dot{r}}{1 - D} \quad (10) $$

$$ f = \frac{\mathrm{eq}(\sigma)}{1 - D} - R(r) - R_0 \quad (11) $$

$$ \text{Active } (\dot{r}) > 0 \text{ if } f > 0, \langle f \rangle_+ = H \dot{r}^m \quad (12) $$

In these constitutive equations, $C_e$, $C_{v_i}$ and $V_{v_i}$ are the elastic stiffness tensor, viscoelastic stiffness tensor and a viscosity tensor respectively. Those are classically formulated for bulk isotropic materials and are defined by the Young’s modulus $E_e$ and $E_{v_i}$ and the viscosity $\eta_{v_i}$ respectively, as well as the Poisson’s ratio $\nu$. The value of $\nu$ presented in Table 1 holds for the elastic response as for all Kelvin–Voigt branches. Moreover, $R_0$ represents the initial yield threshold. $K$ and $n$ are hardening parameters. This model assumes isotropic hardening, which simplifies the representation of the material’s response. While kinematic hardening can be significant in cyclic loading scenarios, particularly for polymers [57–61], its experimental verification is challenging due to complications such as buckling in compression tests. For this reason, the current work focuses on isotropic hardening parameters, as they provide a more straightforward approach to capturing the primary viscoplastic behavior. However, the influence of tension-compression asymmetry and the hydrostatic sensitivity of the viscoplastic response can significantly influence the polymer’s behavior. These factors have been addressed in other studies [60]. $H$ and $m$ are viscous parameters. $S$ and $\beta$ are damage parameters. It should be mentioned that the damage-related parameter $\beta$, typically positive for metallic materials with slow initial and rapid later-stage damage propagation, is set to negative values to capture the opposite trend observed in thermoplastic polyamide, where damage initiates quickly but grows slowly. A damage saturation function $\Lambda(D)$ is employed in the matrix phase to ensure that the damage coefficient $D$ does not exceed its limit. The function is defined as shown in Eq. (13), where the initial and ultimate damage saturation limits, $D_{ini} = 0.45$ and $D_{ult} = 0.6$, are determined from a tensile test on a polyamide specimen. This approach prevents numerical and convergence issues in simulations.

$$ \Lambda(D) = \begin{cases} 1, & D \le D_{\text{ini}}, \\ \sqrt{1 - \left(\frac{D - D_{\text{ini}}}{D_{\text{ult}} - D_{\text{ini}}}\right)^2}, & D_{\text{ini}} < D < D_{\text{ult}}, \\ 0, & D = D_{\text{ult}}. \end{cases} \quad (13) $$

Furthermore, note that $\mathrm{eq}(\sigma)$ denotes the von Mises stress, while $\mathbf{Dev}(\sigma)$ stands for the deviatoric part of the stress. It is worth noticing that, in this model, viscoplasticity and damage evolve via a common multiplier $r$, which stands for the hardening variable. This is done to
<table><caption>**Table 1**
Parameters for the constitutive model of the thermoplastic matrix, considering 50% of relative humidity and room temperature as environmental conditions.</caption><thead><tr><th>Feature</th><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="4">Viscoelasticity</td><td>Single spring</td><td>E<sub>e</sub></td><td>2731</td><td>MPa</td></tr><tr><td rowspan="3">Kelvin-Voigt branch 1</td><td>E<sub>v1</sub></td><td>8766</td><td>MPa</td></tr><tr><td>η<sub>v1</sub></td><td>1395</td><td>MPa.s</td></tr><tr><td>τ<sub>v1</sub> = η<sub>v1</sub>/E<sub>v1</sub></td><td>0.16</td><td>s</td></tr><tr><td rowspan="3">Kelvin-Voigt branch 2</td><td>E<sub>v2</sub></td><td>13 754</td><td>MPa</td></tr><tr><td>η<sub>v2</sub></td><td>165 601</td><td>MPa.s</td></tr><tr><td>τ<sub>v2</sub> = η<sub>v2</sub>/E<sub>v2</sub></td><td>12.04</td><td>s</td></tr><tr><td rowspan="3">Kelvin-Voigt branch 3</td><td>E<sub>v3</sub></td><td>15 010</td><td>MPa</td></tr><tr><td>η<sub>v3</sub></td><td>457 955</td><td>MPa.s</td></tr><tr><td>τ<sub>v3</sub> = η<sub>v3</sub>/E<sub>v3</sub></td><td>30.51</td><td>s</td></tr><tr><td rowspan="3">Kelvin-Voigt branch 4</td><td>E<sub>v4</sub></td><td>11 634</td><td>MPa</td></tr><tr><td>η<sub>v4</sub></td><td>1 307 516</td><td>MPa.s</td></tr><tr><td>τ<sub>v4</sub> = η<sub>v4</sub>/E<sub>v4</sub></td><td>112.39</td><td>s</td></tr><tr><td colspan="2">Poisson ratio</td><td>ν</td><td>0.3</td><td>-</td></tr><tr><td rowspan="5">Viscoplasticity<br/>coupled to damage</td><td>Yield threshold</td><td>R<sub>0</sub></td><td>4.86</td><td>MPa</td></tr><tr><td rowspan="2">Hardening function</td><td>K</td><td>1304.33</td><td>MPa</td></tr><tr><td>n</td><td>0.674</td><td>-</td></tr><tr><td rowspan="2">Viscous stress function</td><td>H</td><td>47.35</td><td>MPa.s<sup>m</sup></td></tr><tr><td>m</td><td>0.068</td><td>-</td></tr><tr><td rowspan="2">Damage</td><td>S</td><td>21.607</td><td>MPa</td></tr><tr><td>β</td><td>-1.105</td><td>-</td></tr></tbody></table>

transcribe the fact that, in thermoplasics, damage is known to be closely connected to yielding [55], making those two mechanisms directly coupled.

The material parameters listed in Table 1 were identified from the polyamide matrix by following a dedicated procedure [46,56,62], considering 50% of relative humidity and room temperature. The identification process employs the Levenberg-Marquardt algorithm [63] to fit the model parameters to experimental data, with the experimental procedure and detailed methodology thoroughly reported in [46]. However, although these environmental conditions are retained in this work, note that other sets of parameters can be found in [41,60] for different environmental conditions in terms of relative humidity, using the Nelder-Mead approach.

## 4.2. Constitutive model for the strands

As discussed earlier in Section 3.2.2, the strands are regarded as equivalent to a unidirectional composite. In such a material, the presence of long unidirectional fibers makes the longitudinal response, i.e., tension 11, linear elastic up to brittle failure due to fibers breakage. However, the transverse response, i.e., tension 22 and/or in-plane shear 12, exhibits a much more progressive degradation. The latter is caused by the development of micro-cracks initiating by debonding at the fiber/matrix interfaces and then propagating by coalescence. To describe this phenomenology, the constitutive model specifically developed by [31,56] for these type of materials is utilized to represent the strand phase in the present multi-scale modeling approach. As for the previous section, a brief overview of this constitutive model is given in here. However, reader are referred to the above-mentioned references for a comprehensive description.

This constitutive model is based on a thermodynamically-consistent formulation, while some features are based on mean-field micromechanics. Indeed, this model introduces an anisotropic damage via a micromechanical description of a RVE containing a micro-crack density $\gamma_c$, as illustrated in Fig. 8. As a known phenomenon in this type of material, damage is accompanied by permanent strains caused by the non-closure effect of the micro-cracks. This aspect is then described through an inelastic strain $\epsilon_s$.

The constitutive equations of this model consider the total strain $\epsilon$ as the only external state variable, whereas $\gamma_c$ and $\epsilon_s$ act as internal state

variables. In these constitutive equations, the initial fourth order stiffness tensor $C_0$, defined as transversely isotropic, is gradually reduced by a tensor $D(\gamma_c)$. The latter is defined as a function of the micro-crack density $\gamma_c$ using mean-field micromechanics [64]. Under this configuration, the transverse shear and in-plane shear stiffnesses are the most affected by the presence of micro-cracks, while the longitudinal stiffness is nearly unaffected.

Below is a summary of the key constitutive equations:

$$ \sigma = [C_0 - D(\gamma_c)] : (\epsilon - \epsilon_s) \quad (14) $$

$$ \dot{\epsilon}_s = \frac{\partial H_s(\sigma)}{\partial \sigma} \dot{\gamma}_c \quad (15) $$

The local mean stress, $\sigma_0$, in the virgin part of the material is expressed as:

$$ \sigma_0 = C_0 : \epsilon_0 \quad \text{with} \quad \epsilon_0 = A_0(\gamma_c) : (\epsilon - \epsilon_s) \quad (16) $$

where $A_0$ stands for the fourth order strain concentration tensor [31].

The criterion governing damage initiation and development is formulated as:

$$ H_c(\sigma_0) = \sqrt{\left(\frac{\sigma_{022}}{R_{22}}\right)^2 + \left(\frac{\sigma_{012}}{R_{12}}\right)^2} \quad (17) $$

The criterion governing the direction of the inelastic strain is defined as:

$$ H_s(\sigma) = \sqrt{(a_{22} \sigma_{22})^2 + (a_{12} \sigma_{12})^2} \quad (18) $$

and micro-crack density, $\gamma_c$, given as:

$$ \gamma_c = \gamma_c^\infty \left[ 1 - \exp\left( - \left[ \frac{\langle \sup(H_c) - 1 \rangle_+}{S} \right]^\beta \right) \right] \quad (19) $$

Since conducting mechanical tests on isolated strands is nearly impossible, it remains very challenging to identify the parameters of the strands constitutive model. Nonetheless, a set of parameters were previously identified by [31,65] for yarns within a polyamide based woven composite, which can be reasonably considered as close to the strands material. Starting from the values reported in these works as initial guess, the parameters of the model were adjusted to match the macroscopic responses of the recycled glass fiber reinforced polyamide 6 composite. The obtained value are listed in Table 2.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Principle of the constitutive model for the strands: (i) Initially transversely isotropic medium. (ii) Embedding a micro-crack density $\gamma_c$ into this medium. (iii) Evaluation of the stiffness reduction $D(\gamma_c)$ induced by $\gamma_c$ via mean-field micromechanics.</figcaption></figure>

<table><caption>**Table 2**
Parameters for the constitutive model of the strands.</caption><thead><tr><th>Feature</th><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="4">Transversely isotropic stiffness tensor<br>(non-null components)</td><td>C<sub>0,111</sub></td><td>45,365</td><td>MPa</td></tr><tr><td>C<sub>0,112</sub> = C<sub>0,113</sub></td><td>6745</td><td>MPa</td></tr><tr><td>C<sub>0,222</sub> = C<sub>0,223</sub></td><td>22,101</td><td>MPa</td></tr><tr><td>C<sub>0,223</sub></td><td>7742</td><td>MPa</td></tr><tr><td>C<sub>0,122</sub> = C<sub>0,133</sub></td><td>8193</td><td>MPa</td></tr><tr><td>C<sub>0,223</sub> = 1/2(C<sub>0,222</sub> - C<sub>0,223</sub>)</td><td>7179</td><td>MPa</td></tr><tr><td>Pure transverse tension threshold</td><td>R<sub>22</sub></td><td>40</td><td>MPa</td></tr><tr><td>Pure in-plane shear threshold</td><td>R<sub>12</sub></td><td>20</td><td>MPa</td></tr><tr><td>Weibull length parameter</td><td>S</td><td>12.3</td><td>-</td></tr><tr><td>Weibull exponent parameter</td><td>β</td><td>2.75</td><td>-</td></tr><tr><td>Micro-cracks saturation</td><td>γ<sub>c</sub><sup>∞</sup></td><td>0.025</td><td>-</td></tr><tr><td>Transverse tension anelasticity parameter</td><td>a<sub>22</sub></td><td>3.60</td><td>-</td></tr><tr><td>In-plane shear anelasticity parameter</td><td>a<sub>12</sub></td><td>2.15</td><td>-</td></tr></tbody></table>

## 5. Macroscopic response of recycled thermoplastic composites

The multi-scale model for recycled thermoplastic is obtained from the generated RVEs in which the previously discussed constitutive models are introduced to describe the behavior of the matrix and the strands. Note that these constitutive models both make use of local damage formulations, which cannot be safely employed beyond the strain softening onset due to well-known strain localization issues [66]. In this study, stopping the analyses when such critical point is attained would be very restrictive as it may occur at the very first loading stages at a few locations with high stress concentrations. Therefore, to overcome this, damage in the matrix is virtually prevented from exceeding the value of 0.6 at which the strain softening usually occurs with the identified parameters in Table 1. In the strands, the damage saturation parameter $\gamma_c^\infty$ is set to a value at which there is no strain softening. Besides, all the simulations presented here were performed at moderate macroscopic strain levels so that these critical damage levels are only reached in very limited parts of the RVE having a minor influence on the macroscopic response of the whole composite.

In this section, the multi-scale model is employed to reproduce uni-axial tests performed on three types of specimens, each one of them having the reference direction oriented at 0°, 45° and 90°, as shown in Fig. 9. As detailed in Section 2.2, the examined recycled composite exhibits a macroscopically heterogeneous structure, posing a challenge in determining the optimal specimen size to ensure representative macroscopic mechanical response. Previous work addressed this by varying the widths of tensile specimens [38]. Findings indicate that the influence of specimen width between 25 mm and 40 mm is minimal compared to factors like specimen orientation. Specimens were cut off from 5 mm thick plates using a water jet hyperbaric machining process. To maintain consistent testing conditions, all specimens were

conditioned to 50% of relative humidity<fnref n="1" /> before being tested at room temperature, specifically at 20 ± 1 °C. It is worth mentioning that the viscoelastic behavior of the material can be heavily influenced by temperature [67]; however, this potential effect was not considered within the scope of the current study. To comprehensively understand the mechanical performance of the recycled composite, several tests are conducted. Firstly, cyclic tests are performed to assess the material's response to repeated loading and unloading cycles. These tests provide insights into the material's damping properties, energy dissipation during hysteresis, and the accumulation of damage over multiple cycles. Secondly, creep and strain recovery tests are conducted to evaluate the material's viscoelastic time-dependent deformation behavior. These tests aim at characterizing how the material responds to sustained loading and its ability to recover deformation over time.

### 5.1. Monotonic loading: effects of strand positions and orientations

To model the response of the recycled composite under monotonic loading, the matrix is represented using a viscoelastic viscoplastic with ductile damage constitutive model (discussed in Section 4.1). Damage initiated by micro-cracks within the strands is addressed using an anisotropic damage constitutive model detailed in Section 4.2. While the previous version of the model exhibits good repeatability in predicting effective properties, nonlinear behavior modeling may introduce variability in results. Therefore, for this work, 10 RVEs are generated,

<footnote n="1">¹ The conditioning was done by initially drying the specimens in a chamber at 80 °C. The specimens were next placed in a hygrometric chamber at an elevated temperature to accelerate the water uptake before lowering the temperature to near-room levels to ensure a uniform water distribution throughout the specimen thickness.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Recycled composite tensile tests specimens cut in 0°, 45° and 90° orientations.</figcaption></figure>

with random strand positions and according to the orientation distribution function (ODF) shown in Fig. 6. It is important to note that the current generation process selects orientation angles randomly within predefined 20° intervals to maintain fidelity to the overall ODF. For instance, during the initial strand deposition, a random angle between 0° and 20° is chosen until the specified number of strands for that interval is deposited, before moving on to the next interval.

The analyses are strain controlled, with the RVEs loaded to 1.2% strain in both 0° and 90° orientations, matching the strain at failure for most tested specimens. The results are plotted as stress-strain curves in Fig. 10(a) and compared to previously performed experimental tensile tests on RTPC specimens [38]. A mean response curve is plotted for numerical data, along with a 95% confidence interval represented by a gray shaded area to assess variability. The model accurately captures the experimental trends, especially the anisotropy between the 0° and 90° orientations through the assigned ODF. The 0° orientation results fall within the range of experimental curves, while the 90° orientation results show a slightly more rigid response. The model shows acceptable levels of variability and dispersion, with numerical variability being less than the experimental. The 90° orientation results show little variability compared to the more significant variability in the 0° orientation. This could be due to the ODF selecting angles randomly within predefined 20° intervals.

To better understand the source of this variability, 10 additional RVEs are generated using a new orientation protocol. Instead of randomly selecting angles, a fixed ODF with known angle values that best fit the mean experimental data is used. This reduces uncontrolled parameters to only the random spatial position of the strands. The results are presented in Fig. 10(b) as stress-strain curves, showing significantly reduced variability, especially for the 0° orientation results. These adjustments lead to the understanding that variability is heavily influenced by the generation ODF. This conclusion can be exploited in model development by taking into account experimental variability and adjusting the orientation protocol to match it.

## 5.2. Loading and unloading cycles at increasing strain levels

Cyclic loading-unloading tests were conducted at four increasing strain levels: $\epsilon = 0.3\%$, $\epsilon = 0.6\%$, $\epsilon = 0.9\%$, and $\epsilon = 1.2\%$. Loading and unloading stages were controlled with a strain rate of approximately $\dot{\epsilon} = 8 \times 10^{-3} \text{s}^{-1}$. It should be noted that strain cannot be perfectly regulated by the testing machine. Hence, the strain rate is not perfectly uniform during the test but is generally close to the indicated value. This strain rate was chosen to stay in accordance with the values used to calibrate the constitutive model, namely, the characteristic times of Kelvin-Voigt branches. During unloading cycles, the load was reduced until reaching a threshold value close to $F = 0 \text{ N}$ (here $F = 50 \text{ N}$) to prevent specimen compression. Fig. 11(a) shows the applied strain over time that was recorded by the extensometer for the specimen cut along the 0° direction, while Fig. 11(b) illustrates the stress response highlighting the stress-controlled unloading steps.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Stress-strain curves for RVEs with (a) random and (b) fixed orientation distributions compared to experimental data, showing significant improvement in variability with fixed orientation.</figcaption></figure>

Fig. 12 shows the strain-stress responses for 0°, 45°, and 90° specimen orientations, illustrating the anisotropy observed in prior monotonic tests [38]. A significant stiffness discrepancy between the 0° and 90° orientations is evident in Fig. 12(a), consistent with monotonic tensile test results. The 0° specimens exhibit the highest stiffness, the 90° specimens the lowest, and the 45° specimens show intermediate stiffness (Figs. 12(b) and 12(c)). This anisotropy supports assumptions of preferential fiber orientation induced by the recycling process, aligning with typical behavior in short fiber composites where fibers predominantly align in a direction that maximizes stiffness, while the transverse direction exhibits minimal stiffness [68].

Hysteresis loops, indicating energy dissipation during loading-unloading, are most pronounced at 90°, indicating higher internal friction and energy loss. At 45°, hysteresis loops are less prominent, suggesting moderate energy dissipation. At 0°, hysteresis loops are minimal, reflecting minimal energy loss and more elastic behavior.
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Strain recorded by the extensometer for a specimen cut along 0° direction and its stress response upon loading unloading stages with increasing strain levels of ε = 0.3%, ε = 0.6%, ε = 0.9%, and ε = 1.2%.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Strain-stress curves for cyclic loading and unloading tests for various specimen orientations, demonstrating anisotropy consistent with prior monotonic tests. (a) Discrepancy between 0° and 90° orientations, with 0° showing maximum stiffness and 90° showing minimal stiffness. (b) and (c) Intermediate stiffness at 45°, illustrating a gradual transition between 0° and 90°.</figcaption></figure>
The experimental loading paths were reproduced numerically on five RVEs, with loadings at the four increasing strain levels: $\epsilon = 0.3\%$, $\epsilon = 0.6\%$, $\epsilon = 0.9\%$, and $\epsilon = 1.2\%$ and unloadings at zero stress. To remain consistent with the experimental data, the exact loading paths of Fig. 11 were applied to the model to ensure proper strain rates. These RVEs were generated following the previously established ODF (Fig. 6), with orientation angles randomly set within predefined 20° intervals to reproduce the experimental variability in the numerical results over multiple RVEs.

Simulation results are presented in Fig. 13 in the form of stress-strain curves for three material orientations (0°, 45°, and 90°) across five RVEs, compared against experimental data (dashed lines). As expected, material response varies significantly among orientations. The 0° orientation (Fig. 13(a)) exhibits a more rigid response compared to 45° (Fig. 13(c)) and 90° (Fig. 13(e)). For the 90° orientation, numerical results closely match experimental findings with less variability across different RVEs compared to the 0° and 45° orientations, where responses varied more among RVEs. Simulations generally align closely with experiments across all orientations, despite minor discrepancies mainly for the 0° and 45° orientations. These discrepancies are acceptable for these materials, where the recycling process strongly influences the subsequent microstructure, resulting in significant variability in experimental results.

The viscoelastic viscoplastic matrix model effectively captures hysteresis loops, indicating energy dissipation in both matrix and strands, crucial for simulating stiffness reduction from accumulated damage. Specifically, the 90° orientation shows increased energy dissipation per cycle, with notable hysteresis loops in experiments. However, the model's ability to replicate this higher dissipation is limited, partly due to insufficient representation of viscous effects. Adding Kelvin-Voigt branches with varying characteristic times could enhance the model's ability to replicate these effects. Moreover, the current model does not consider thermomechanical couplings, which could incorporate temperature-dependent material property changes, potentially improving overall accuracy.

After analyzing the material's average macroscopic response, the local fields at the microscopic level were evaluated (Fig. 13). The analysis focused on the distribution of damage, assessed through the micro-crack density, denoted as $\gamma_c$. The results showed that the highest micro-crack densities occurred in the 0° loading direction, followed by the 45° and 90° directions.

For the 0° and 90° loading directions, the highest micro-crack densities were observed on strands with local orientations closer to 45°. Specifically, for the 0° loading direction, the maximum damage occurred on a strand with a local material orientation of -37.58°, as shown in Fig. 13(b). Similarly, for the 90° loading direction, the highest damage was found on a strand with a local orientation of 39.49°, as shown in Fig. 13(f). In contrast, for the 45° loading direction, the greatest damage was recorded on a strand with a local orientation closer to 0° (0.022°), as depicted in Fig. 13(d).

From a modeling perspective, although identifying the exact reasons for these observations is challenging due to the interplay of multiple factors, one possible explanation is that damage is primarily initiated by in-plane shear. This is because, in this model configuration, the shear activation threshold ($R_{12} = 20$ MPa) is lower than the tensile activation threshold ($R_{22} = 40$ MPa), as shown in Table 2.

### 5.3. Creep and strain recovery

Creep and strain recovery tests complement cyclic loading-unloading tests by evaluating long-term deformation and strain recovery behaviors of materials. These tests provide insights into both viscoelastic and viscoplastic properties that load-unload tests may not fully capture. The loading and unloading phases were conducted at a strain rate of $\dot{\epsilon} = 2 \times 10^{-3} \text{ s}^{-1}$. During the creep stage, the material was subjected to 50% of its ultimate stress ($\sigma_{ult}$) with respect to each orientation. This

precaution was taken to prevent the specimen from breaking in the middle of the test, particularly considering the significant variability in these values for the studied recycled composite. In the recovery step, the load was reduced to 50 N to prevent compression upon unloading, and this load was maintained for 5 min for both the creep and recovery stages. The applied stresses with respect to time are plotted in Figs. 14(a), 14(b) and 14(c) for the 0°, 45° and 90° orientations respectively.

The results of creep and strain recovery tests are presented in Fig. 15 for the three orientations. Due to the significant variability observed, an average curve is plotted to illustrate the mean response. Stress-strain curves represent the overall behavior of the tested specimens, while strain responses highlight the creep behavior after load stabilization and strain recovery upon load removal. In the 0° orientation (Figs. 15(a) and 15(b)), significant variability is observed among specimens. In two cases, noticeable creep behavior is observed, with strain gradually increasing under a constant load, while the stiffest specimen exhibits less creep. Specimens with greater creep have more open hysteresis loops. The 90° orientation (Figs. 15(e) and 15(f)) exhibits more pronounced creep response than the 0°, with significant strain accumulation and wider hysteresis loops, suggesting greater energy dissipation. The 45° orientation (Figs. 15(c) and 15(d)) shows consistent results for all specimens, with creep behavior similar to the 90° orientation.

After demonstrating the model's ability to capture phenomena in cyclic loading and unloading, including energy dissipation and permanent stiffness reduction, the creep and strain recovery numerical response of the RTPC model is examined. Experimental loading paths in Fig. 14 were reproduced numerically on the same five RVEs. Simulations were stress-controlled by rapidly loading to 50% of the ultimate stress ($\sigma_{ult}$) for each orientation, maintaining the load for 5 min, then rapidly unloading to near zero stress and holding for another 5 min. The focus is on how strain accumulates under constant load, recovers at near zero stress, and how well the modeling results match experimental findings.

Simulation results are presented as stress-strain curves in Fig. 16, showing the overall response for all RVEs across the three material orientations (0°, 45°, and 90°) and comparing them to experimental data (dashed lines) in the Figs. 16(a), 16(c) and 16(e). Strain evolution curves for the three orientations are also plotted to observe strain accumulation during creep and its recovery after unloading (Figs. 16(b), 16(d) and 16(f)). A mean response curve is plotted for numerical data, along with a 95% confidence interval represented by a gray shaded area, to assess variability and compare it across material orientations.

After further investigation of stress-strain curves, the results for the 0° orientation show that the 5 RVEs fall within the experimental curves' upper and lower bounds. Expected variability is observed in the same magnitude for both monotonic and cyclic loading and unloading cases, as detailed in Sections 5.1 and 5.2. For the 45° orientation, numerical results slightly overestimate the experimental response on average, with less variability compared to the 0° orientation. Meanwhile, for the 90° orientation, numerical results are very close to the experimental data, showing minimal variability across RVEs. Moreover, significant energy dissipation and stiffness loss are observed in the experimental results for the 45° and 90° orientations, which the numerical model does not perfectly capture.

Focusing on strain evolution curves, for the 0° orientations, simulation results align with the main experimental trends, with the 95% confidence interval spectrum of numerical results falling within the experimental data range. Strain accumulation stabilizes similarly in simulation results after loading. In contrast, for the 45° orientation, numerical results show lower strain levels compared to experiments, where significant strain accumulation is observed. In the case of the 90° orientation, numerical results generally match the experimental data
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Comparison of numerical and experimental results for cyclic loading and unloading at increasing strain levels for various specimen orientations: (a) 0°, (b) 45° and (c) 90°. Loading and unloading stages were controlled with a strain rate of $\dot{\epsilon} = 8 \times 10^{-3} \text{ s}^{-1}$. Comparison of local damage distribution assessed through micro-crack density $\gamma_c$ for the three loading direction: (d) 0°, (e) 45° and (f) 90°.</figcaption></figure>

range, though strain accumulation was more pronounced in some specimens. Regarding strain recovery, the model consistently shows higher recovery than observed experimentally across all material orientations. As anticipated, variability in numerical results decreases relative to material orientation, indicated by a narrower 95% confidence interval for 45° compared to 0°, and even narrower for 90° compared to both.

Overall, considering the variability in mechanical response due to RTPC's complex microstructure and loading configurations, this

study's modeling efforts yield satisfactory results. Attempts to reproduce experimental variability generally align within acceptable ranges. However, the model struggles to capture significant energy dissipation during cyclic loadings and creep-recovery phases, especially for 45° and 90° orientations, due to insufficient representation of viscous effects. Enhancing the model with Kelvin-Voigt branches of varying characteristic times could address these limitations. Additionally, the current model lacks consideration of thermomechanical couplings, which could
<figure><img src="image_16.png" /><figcaption>Fig. 14. Applied stresses over time during creep and strain recovery tests for different material orientations: (a) 0°, (b) 45°, and (c) 90°.</figcaption></figure>

incorporate temperature-dependent material property changes and potentially enhance accuracy. Notably, the viscoelastic-viscoplastic with ductile damage matrix model used for this work was calibrated based on virgin polymer data. Conducting new experiments with polymer matrices subjected to processing conditions similar to those in recycled composites could yield more appropriate model parameters, thereby significantly improving results.

## 6. Conclusions and further work

This study presents a novel numerical and experimental investigation into the inelastic nonlinear response of a mechanically recycled composite. The material undergoes cyclic loading-unloading and creep-recovery tests to characterize its damping properties, energy dissipation during hysteresis, and time-dependent behavior. Experimental results indicate anisotropy consistent with previous monotonic tensile tests, likely due to preferential fiber orientation induced during the molding process. Pronounced hysteresis loops are more apparent in specimens cut perpendicularly to the production direction. Significant variability is observed in creep-recovery tests, with some specimens exhibiting stiffer responses and others displaying strain accumulation upon load stabilization and less strain recovery after load removal.

The nonlinear behavior is attributed mainly to the polymer matrix's viscoelastic and viscoplastic mechanisms. A viscoelastic viscoplastic with continuum damage constitutive model is considered for the matrix to capture these effects, and an anisotropic damage model is

applied to the strands to account for micro-cracks under loading. The proposed model effectively captures the monotonic loading response, and multiple RVEs have been generated to study the effects of strand positions and orientations. The generation process produces numerical variability less than the experimental dispersion, suggesting variability is influenced by the generation ODF.

For cyclic loading-unloading tests, the model effectively captures the experimentally observed anisotropy and the general trends of hysteresis loops induced by dissipation and stiffness reduction from accumulated damage. For creep-recovery responses, the model's results generally fall within the experimental data range, though high levels of dissipation and stiffness loss in certain orientations are captured less accurately, likely due to insufficient representation of viscous effects. Overall, considering the variability in mechanical response due to the recycled composite's complex microstructure and loading configurations, the study's modeling efforts yield satisfactory results, generally aligning within acceptable ranges of experimental variability.

Further research should focus on modeling the recycling process to complete the process-microstructure-property linkage. This includes understanding the impact of processing parameters on the resulting microstructure after recycling and using the proposed model to predict the mechanical response. To accelerate the development of recycled composites, the proposed nonlinear multiscale model could inform data-driven approaches based on neural networks [69–71]. The authors have developed similar modeling tools for woven composites, offering a solid thermodynamic foundation that enhances prediction reliability
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Creep recovery test results for various specimen orientations, illustrating variability and mean behavior. (a) and (b) 0° orientation: noticeable creep in some specimens, less prominent in the stiffest one. (c) and (d) 45° orientation: showing pronounced creep. (e) and (f) 90° orientation: consistent results similar to 45°, significant strain accumulation and increased energy dissipation.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Stress-strain and strain evolution curves for various specimen orientations, comparing simulation results to experimental data. (a), (c), (e) Stress-strain curves for 0°, 45°, and 90° orientations: simulations fall within experimental bounds, with minimal variability. (b), (d), (f) Strain evolution curves: simulations align with experimental trends, showing stabilized strain accumulation and higher strain recovery. Variability decreases with orientation, indicated by narrower confidence intervals.</figcaption></figure>

and admissibility [65,72]. This approach could be integrated with commercial finite element codes for multi-level simulations of structures. Applying this to recycled composites would involve defining appropriate Internal State Variables (ISVs) to capture the thermodynamical states of the strands with respect to their orientations. Moreover, temperature and strain rate sensitivity significantly affect the behavior of polymers, as shown in previous studies [73–75]. In future work, incorporating these effects within the proposed modeling framework will further enhance the predictive capability of the material response, particularly under varying environmental conditions.

CRediT authorship contribution statement

S.E. Sekkal: Writing – original draft, Validation, Methodology, Investigation, Formal analysis. F. Meraghni: Writing – review & editing, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. G. Chatzigeorgiou: Writing – review & editing, Investigation, Formal analysis, Conceptualization. F. Praud: Writing – review & editing, Investigation, Formal analysis. N. Durand: Resources, Investigation, Data curation.
Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

This work has been conducted in the framework of the project
UpCyclComp-CaMéX-IA consortium (Campus des Métiers et des Quali-
fications d'Excellence) de la région Grand Est funded in the framework
of PIA3 (2021-2028). The authors gratefully acknowledge CETIM for
supporting this work, providing the test materials along with the shared
knowledge and practical procedures. They also acknowledge the sup-
port of the High Performances Computation Center (HPC - Cassiopée)
of Arts et Métiers Institute of Technology, which was made available
to carry out the research and the related numerical results reported in
this paper.

Data availability

Data will be made available on request.

References

[1] Asmatulu E, Twomey J, Overcash M. Recycling of fiber-reinforced composites and direct structural composite recycling concept. J Compos Mater 2014;48:593-608.

[2] Khurshid MF, Hengstermann M, Hasan MMB, Abdkader A, Cherif C. Recent developments in the processing of waste carbon fibre for thermoplastic composites—A review. J Compos Mater 2020;54:1925-44.

[3] Thomason JL, Yang L, Meier R. The properties of glass fibres after conditioning at composite recycling temperatures. Compos Part A: Appl Sci Manuf 2014;61:201-8.

[4] Kiss P, Stadlbauer W, Burgstaller C, Stadler H, Fehringer S, Haeuserer F, Archodoulaki V-M. In-house recycling of carbon- and glass fibre-reinforced thermoplastic composite laminate waste into high-performance sheet materials. Compos Part A: Appl Sci Manuf 2020;139:106110.

[5] Krauklis AE, Karl CW, Gagani AI, Jørgensen JK. Composite material recycling technology—State-of-the-art and sustainable development for the 2020s. J Compos Sci 2021;5:28.

[6] Oliveux G, Dandy LO, Leeke GA. Current status of recycling of fibre reinforced polymers: Review of technologies, reuse and resulting properties. Prog Mater Sci 2015;72:61-99.

[7] Rasheed MIA. Compression molding of chopped woven thermoplastic composite flakes. 2016.

[8] Pakdel E, Kashi S, Varley R, Wang X. Recent progress in recycling carbon fibre reinforced composites and dry carbon fibre wastes. Resour Conserv Recycl 2021;166:105340.

[9] Mazzocchetti L, Benelli T, D'Angelo E, Leonardi C, Zattini G, Giorgini L. Validation of carbon fibers recycling by pyro-gasification: The influence of oxidation conditions to obtain clean fibers and promote fiber/matrix adhesion in epoxy composites. Compos Part A: Appl Sci Manuf 2018;112:504-14.

[10] Khalil Y. Comparative environmental and human health evaluations of thermol-
ysis and solvolysis recycling technologies of carbon fiber reinforced polymer
waste. Waste Manage 2018;76:767-78.

[11] Gonçalves RM, Martinho A, Oliveira JP. Recycling of reinforced glass fibers waste: Current status. Materials 2022;15:1596.

[12] Colucci G, Ostrovskaya O, Frache A, Martorana B, Badini C. The effect of mechanical recycling on the microstructure and properties of PA66 composites reinforced with carbon fibers. J Appl Polym Sci 2015;132(29).

[13] Franzen B, Klason C, Kubat J, Kitano T. Fibre degradation during processing of short fibre reinforced thermoplastics. Composites 1989;20(1):65-76.

[14] Subasinghe ADL, Das R, Bhattacharyya D. Fiber dispersion during compound-
ing/injection molding of PP/kenaf composites: Flammability and mechanical
properties. Mater Des 2015;86:500-7.

[15] Tapper RJ, Longana ML, Hamerton I, Potter KD. A closed-loop recycling process for discontinuous carbon fibre polyamide 6 composites. Compos Part B: Eng 2019;179:107418.

[16] Schinner G, Brandt J, Richter H. Recycling carbon-fiber-reinforced thermoplastic composites. J Thermoplast Compos Mater 1996;9:239-45.

[17] Eguemann N, Giger L, Roux M, Dransfeld C, Thiebaud F, Perreux D. Compression moulding of complex parts for the aerospace with discontinuous novel and recycled thermoplastic composite materials. In: 19th international conference on composite materials. 2013. p. 1-11.

[18] Howell DD, Fukumoto S. Compression molding of long chopped fiber thermoplastic composites. 2014, p. 13-6.

[19] Kravchenko SG, Sommer DE, Denos BR, Avery WB, Pipes RB. Structure-property relationship for a prepreg platelet molded composite with engineered meso-morphology. Compos Struct 2019;210:430-45.

[20] Kravchenko SG, Sommer DE, Denos BR, Favaloro AJ, Tow CM, Avery WB, Pipes RB. Tensile properties of a stochastic prepreg platelet molded composite. Compos Part A: Appl Sci Manuf 2019;124:105507.

[21] Sattar S, Laredo BB, Kravchenko SG, Kravchenko OG. Effect of platelet length and stochastic morphology on flexural behavior of prepreg platelet molded composites. Polym Compos 2023;44(4):2122-37.

[22] Harper LT, Qian CC, Luchoo R, Warrior NA. 3D geometric modelling of discontinuous fibre composites using a force-directed algorithm. J Compos Mater 2017;51:2389-406.

[23] Tang H, Zhou G, Chen Z, Huang L, Avery K, Li Y, Liu H, Guo H, Kang H, Zeng D. Fatigue behavior analysis and multi-scale modelling of chopped carbon fiber chip-reinforced composites under tension-tension loading condition. Compos Struct 2019;215:85-97.

[24] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21(5):571-4.

[25] Hill R. A self-consistent mechanics of composite materials. J Mech Phys Solids 1965;13(4):213-22.

[26] Nemat-Nasser S, Hori M. Micromechanics: overall properties of heterogeneous materials. Elsevier; 2013.

[27] Castañeda PP, Willis JR. The effect of spatial distribution on the effective behavior of composite materials and cracked media. J Mech Phys Solids 1995;43(12):1919-51.

[28] Sanchez-Palencia E. Comportements local et macroscopique d'un type de milieux physiques hétérogènes. Internat J Engrg Sci 1974;12(4):331-51.

[29] Suquet P. Elements of homogenization for inelastic solid mechanics. Homog Tech Compos Media 1987;272:193-278.

[30] Michel J-C, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure: a computational approach. Comput Methods Appl Mech Engrg 1999;172(1-4):109-43.

[31] Praud F, Chatzigeorgiou G, Chemisky Y, Meraghi F. Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites. Compos Struct 2017;182:223-36.

[32] Tikarrouchine E, Benaarbia A, Chatzigeorgiou G, Meraghi F. Non-linear FE2 multiscale simulation of damage, micro and macroscopic strains in polyamide 66-woven composite structures: Analysis and experimental validation. Compos Struct 2021;255:112926.

[33] Tikarrouchine E, Chatzigeorgiou G, Chemisky Y, Meraghi F. Fully coupled thermo-viscoplastic analysis of composite structures by means of multi-scale three-dimensional finite element computations. Int J Solids Struct 2019;164:120-40.

[34] Moulinec H, Suquet P. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Engrg 1998;157(1-2):69-94.

[35] Schneider M. A review of nonlinear FFT-based computational homogenization methods. Acta Mech 2021;232(6):2051-100.

[36] Lucarini S, Upadhyay MV, Segurado J. FFT based approaches in micromechanics: fundamentals, methods and applications. Modelling Simul Mater Sci Eng 2021;30(2):023002.

[37] Nachtane M, Meraghi F, Chatzigeorgiou G, Harper L, Pelascini F. Multiscale viscoplastic modeling of recycled glass fiber-reinforced thermoplastic composites: Experimental and numerical investigations. Compos Part B: Eng 2022;242:110087.

[38] Sekkal SE, Meraghi F, Chatzigeorgiou G, Peltier L, Durand N. Experimental and multi-scale investigation of the mechanical behavior of mechanically recycled glass fiber reinforced thermoplastic composites. Compos Part B: Eng 2023;264:110925.

[39] Arruda EM, Boyce MC, Jayachandran R. Effects of strain rate, temperature and thermomechanical coupling on the finite strain deformation of glassy polymers. Mech Mater 1995;19(2-3):193-212.

[40] Chekkour R, Benaarbia A, Chatzigeorgiou G, Meraghi F, Robert G. Effect of thermo-hydro glycol aging on the damage mechanisms of short glass-fiber reinforced polyamide 66. Compos Part A: Appl Sci Manuf 2023;165:107358.

[41] Satouri S, Chekkour R, Chatzigeorgiou G, Meraghi F, Robert G. Numerical-experimental approach to identify the effect of relative humidity on the material parameters of a rate-dependent damage model for polyamide 66. Mech Mater 2023;184:104735.

[42] Miqoi N, Pomarede P, Meraghi F, Declercq NF, Guillaumat L, Coz GL, Dela-lande S. Detection and evaluation of barely visible impact damage in woven glass fabric reinforced polyamide 6.6/6 composite using ultrasonic imaging, X-ray tomography and optical profilometry. Int J Damage Mech 2021;30:323-48.

[43] Arruda EM, Boyce MC. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. J Mech Phys Solids 1993;41(2):389-412.

[44] Billon N. New constitutive modeling for time-dependent mechanical behavior of polymers close to glass transition: Fundamentals and experimental validation. J Appl Polym Sci 2012;125(6):4390-401.
[45] Krairi A, Doghri I. A thermodynamically-based constitutive model for thermo-plastic polymers coupling viscoelasticity, viscoplasticity and ductile damage. Int J Plast 2014;60:163-81.

[46] Praud F, Chatzigeorgiou G, Bikard J, Meraghni F. Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation. Mech Mater 2017;114:9-29.

[47] Widom B. Random sequential addition of hard spheres to a volume. J Chem Phys 1966;44:3888-94.

[48] Praud F, Schneider K, Chatzigeorgiou G, Meraghni F. Microstructure generation and full-field multi-scale analyses for short fiber reinforced thermoplastics: Application to PA66GF composites. Compos Struct 2024;341:118175.

[49] Schneider K, Klusemann B, Bargmann S. Automatic three-dimensional geometry and mesh generation of periodic representative volume elements for matrix-inclusion composites. Adv Eng Softw 2016;99:177-88.

[50] Schneider K. Computational micromechanics of matrix-inclusion composites [Ph.D. thesis], 2019.

[51] Jacquenot G. Minimum distance between two polygons. 2024, Matlab Central File Exchange. Retrieved June 20, 2024.

[52] Anagnostou D, Chatzigeorgiou G, Chemisky Y, Meraghni F. Hierarchical micromechanical modeling of the viscoelastic behavior coupled to damage in SMC and SMC-hybrid composites. Compos Part B: Eng 2018;151:8-24.

[53] Barral M, Chatzigeorgiou G, Meraghni F, Léon R. Homogenization using modified Mori-Tanaka and TFA framework for elastoplastic-viscoelastic-viscoplastic composites: Theory and numerical validation. Int J Plast 2020;127:102632.

[54] Chen Q, Chatzigeorgiou G, Meraghni F. Extended mean-field homogenization of viscoelastic-viscoplastic polymer composites undergoing hybrid progressive degradation induced by interface debonding and matrix ductile damage. Int J Solids Struct 2021;210:1-17.

[55] Detrez F, Cantournet S, Seguela R. Plasticity/damage coupling in semi-crystalline polymers prior to yielding: Micromechanisms and damage law identification. Polymer 2011;52:1998-2008.

[56] Praud F. Multi-scale modelling of thermoplastic-based woven composites, cyclic and time-dependent behaviour [Ph.D. thesis], Ecole nationale supérieure d'arts et métiers-ENSAM; 2018.

[57] da Costa Mattos H, Reis J, De Medeiros L, Monteiro A, Teixeira S, Chaves E. Analysis of the cyclic tensile behaviour of an elasto-viscoplastic polyamide. Polym Test 2017;58:40-7.

[58] Motta E, Reis J, da Costa Mattos H. Analysis of the cyclic tensile behaviour of an elasto-viscoplastic polyvinylidene fluoride (PVDF). Polym Test 2018;67:503-12.

[59] Motta E, Reis J, da Costa Mattos H. Modelling the cyclic elasto-viscoplastic behaviour of polymers. Polym Test 2019;78:105991.

[60] Satouri S, Chatzigeorgiou G, Meraghni F, Robert G. Viscoelastic-viscoplastic model with ductile damage accounting for tension-compression asymmetry and hydrostatic pressure effect for polyamide 66. Eur J Mech A Solids 2025;110:105491.

[61] Krairi A, Doghri I, Robert G. Multiscale high cycle fatigue models for neat and short fiber reinforced thermoplastic polymers. Int J Fatigue 2016;92:179-92.

[62] Praud F, Chatzigeorgiou G, Meraghni F. Fully integrated multi-scale modelling of damage and time-dependency in thermoplastic-based woven composites. Int J Damage Mech 2021;30:163-95.

[63] Lourakis MI, et al. A brief description of the Levenberg-Marquardt algorithm implemented by levmar. Found Res Technol 2005;4(1):1-6.

[64] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21:571-4.

[65] Idrissi MEF, Praud F, Meraghni F, Chinesta F, Chatzigeorgiou G. Multiscale Thermodynamics-Informed Neural Networks (MuTINN) towards fast and frugal inelastic computation of woven composite structures. J Mech Phys Solids 2024;186:105604.

[66] De Borst R, Sluys L-J, Mühlhaus H-B, Pamin J. Fundamental issues in finite element analyses of localization of deformation. Eng Comput 1993;10(2):99-121.

[67] Kießling R, Ihlemann J. The dependence of the viscoelastic properties of Polyamide 6 on temperature and moisture content: experiment and modeling. Mater Today: Proc 2020;32:83-7.

[68] Mortazavian S, Fatemi A. Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites. Compos Part B: Eng 2015;72:116-29.

[69] Ghavamian F, Simone A. Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network. Comput Methods Appl Mech Engrg 2019;357:112594.

[70] Ghane E, Fagerström M, Mirkhalaf M. Recurrent neural networks and transfer learning for predicting elasto-plasticity in woven composites. Eur J Mech A Solids 2024;105378.

[71] Chen Q, Jia R, Pang S. Deep long short-term memory neural network for accelerated elastoplastic analysis of heterogeneous materials: An integrated data-driven surrogate approach. Compos Struct 2021;264:113688.

[72] Idrissi MEF, Praud F, Champaney V, Chinesta F, Meraghni F. Multiparametric modeling of composite materials based on non-intrusive PGD informed by multiscale analyses: Application for real-time stiffness prediction of woven composites. Compos Struct 2022;302:116228.

[73] Reis J, Pacheco L, da Costa Mattos H. Influence of the temperature and strain rate on the tensile behavior of post-consumer recycled high-density polyethylene. Polym Test 2013;32(8):1576-81.

[74] Pacheco L. Temperature and variable strain rate sensitivity in recycled HDPE [Ph.D. thesis], Ph. D Dissertation, Engineering School, Universidade Federal Fluminense, Brazil; 2014.

[75] Benaarbia A, Chatzigeorgiou G, Kiefer B, Meraghni F. A fully coupled thermo-viscoelastic-viscoplastic-damage framework to study the cyclic variability of the Taylor-Quinney coefficient for semi-crystalline polymers. Int J Mech Sci 2019;163:105128.