Multiscale Viscoelastic Analysis of Plain Weave Textile Composites

Modeling of viscoelastic behavior of heterogeneous materials has been mostly limited to a two-phase material system with reinforcements being elastic. The isotropy of the matrix phase then considerably simplified the choice of a suitable constitute model. In some applica- tions, the macroscopic isotropy derived from homogenization then allowed for adopting the same constitutive model on the macro-scale with model parameters obtained computationally through virtual laboratory tests [1]. Modeling of asphalt mixtures within a fully uncoupled multiscale computational framework is just one particular example [2].

This approach is typically precludedin fibrous composites where the aligned fibers, whether periodic [3, 4] or randomly distributed [5, 6, 7] in the transverse cross-section, generate a higher order material symmetry at the macroscopic level. The need for a fully coupled mod- eling approach then immediately arises. When exploiting the finite element method (FEM) on both the micro- and meso-scale, the analysis is typically performed in the framework of FE2 computational strategy [8, 9, 10]. Although often supported by parallel computing [11, 12, 13], such an analysis may prove computationally very expensive. Classical averaging schemes [14] such as the Mori-Tanaka (MT) method [15] then offer a suitable alternative.

However, whenapplied in itsstandard formulation assuming,ingeneral,ann-pointaver-agingwithconcentrationorlocalizationfactorsofeachphaser=1,...,nderivedonthe basis of elasticity, the method provides estimates ofthe macroscopic response, which are too stiff when compared to finite element simulations. This is generally attributed to an insuf- ficient refinement of local fields unable to represent, for example, localization of inelastic strains driven by microstructure even in the absence of softening. Several routs have been explored to address this issue. Among others, Dvorak’s transformation field analysis (TFA) [16, 17] combined for fibrous composites with the PHA model to calculate the transformation influence functions [18] deserves particular attention. A review of other potential approaches including the affine formulation proposed by Masson et al. [19] can be found in [9]. To keep the analysis as efficient as possible, a relatively simple adjustment of the original two-phase formulation of Hashin-Shtrikman variational principles was proposed in [20] by replacing the linear elastic comparison medium with the one depending on a gradually evolving secant shear modulus of the matrix phase supported by the choice of the generalized Leonov consti- tutive model to represent its viscoelastic response. However, in the light of the MT method, the formulation proposed in [2] and further advanced in [21] proved even more elegant and will be employed also in the present study.

The remainder of the paper is organized as follows. The theoretical framework of both the first-order homogenization method and the MT micromechanical model is briefly outlined in Section 2. This is followed by a short description of the generalized Leonov material model adopted in numerical simulations both at the level of yarns and textile ply presented in Sec- tion 4. The essential findings are finally summarized in Section 5.

As already intimated in the introductory part, predicting the macroscopic viscoelastic response of plain weave textile composites relies on a coupled two-scale analysis. Similar to standard FE2 scheme, we consider a statistically equivalent periodic unit cell (SEPUC) at the level of textile plies discretized, herein, by means of constant strain tetrahedral elements. However, unlike FE2 scheme, the increments of local stress averages for given increments of average strain in each element in the yarn are found from the application of the MT averaging scheme. Both methods are shortly described in a sequel.

Consider a representative volume element (RVE) $\Omega$ in the form of SEPUC loaded on its outer boundary, identified by an outward unit normal $n_i$, by either prescribed displacements $u_i=E_{i j} x_j$ or tractions $p_i=\sum_{i j} n_j$, which in turn generate microscopically uniform strains $E_{i j}$ or stresses $\Sigma_{i j}$ in an equivalent homogeneous medium. The local strain field $\varepsilon_{i j}$ then admits the following decomposition written in the vector-matrix notation as

where $\varepsilon^*(x)$ represents the fluctuation strain associated with the fluctuation part of local displacements $u^*$ which, henceforth, is assumed periodic. Writing the increment of local stresses $\Delta \boldsymbol{\sigma}$ as

allows us to expand the Hill lemma $\left\langle\delta \boldsymbol{\varepsilon}^{\mathrm{I}} \Delta \sigma\right\rangle=\delta \boldsymbol{E}^{\mathrm{I}} \Delta \Sigma(\langle\cdot\rangle$ stands for volume averaging) into the following discretized system of governing equations

where

In equations (2) and (3), $L$ represents an instantaneous stiffness matrix corresponding either to the polymer matrix or to the homogenized stiffness of the yarn provided by the MT scheme, $\Delta \mu$ is the eigenstrain vector here representing the creep strain developed correspondingly in the polymer matrix or the yarn, $\mathbf{B}$ is the standard geometrical matrix, and $\Delta r$ stores the increments of nodal displacements of $\boldsymbol{u}^*$. Further details are available in [14,22].

Suppose that a two-phase composite consisting of elastic aligned fibers embedded into a viscoelastic matrix is loaded by an increment of the mesoscopic strain $\Delta \bar{\varepsilon}$ (this term can be imagined as a strain increment in a specific element of the finite element mesh at the ply level discussed previously in Section 2.1). The local strain increments in the fiber (f) and matrix (m) phase then follow from the TFA analysis as

where $\overparen{\mathbf{A}}_r$ and $\overparen{\mathbf{D}}_{r m}$ are the instantaneous mechanical strain localization factors and strain transformation influence functions, respectively, and are function of instantaneous properties of the matrix phase. In light of the MT method, they are provided by

where the partial strain concentration factor $\widehat{\mathbf{T}}_f$ depends on the shape and orientation of the fiber and instantaneous properties of the matrix [14].

The matrix-phase constitutive equation takes a usual form

where $\widehat{\mathbf{L}}_m$ represents a potential dependence on the current viscoelastic modulus. The local stresses in the fiber phase are provided by

where $\Delta s_f, \Delta \tilde{\sigma}_{f, m}$ represent the deviatoric and mean components of the fiber stress increment $\Delta \sigma_f$. The tensorial notation in Eqs. (12) and (13) is adopted just for the sake of convenience. The increments of mesoscopic stresses $\Delta \bar{\sigma}$ and eigenstrains $\Delta \mu$ follow from

where $c_r, r=f, m$, is the volume fraction of a given phase and $\widehat{\mathbf{B}}_m$ is the so-called stress concentration factor of the matrix phase. The second term in Eq. (14) is called the Levin formula [14].

The damage like parameter $\omega$ was introduced in Eq. (12) 1 in analogy with [2,21] to reduce the stresses carried by the fiber phase. One may associate that with a debonding like failure, which in turn simulates the formation of shear bands in real microstructures once loading the composite beyond the elastic limit. Its evolution as a function of the equivalent deviatoric stress in the matrix phase $\tau_{m, e q}=\sqrt{\frac{1}{2} s_{m, i j} s_{m, i j}}$ is proposed on the basis of the adopted generalized Leonov nonlinear viscoelastic model in the form, note Eq. (17) ${ }_2$,

where $M, N, T$, are the model parameters, $t$ stands for the current time instant, and $\tau_0$ is the parameter of the generalized Leonov model to be introduced later in Section 3.

The viscoelastic behavior of the matrix phase is assumed to be well represented by the generalized Leonov model which, similarly to von Mises plasticity, limits the nonlinear viscoelastic response to the deviatoric stress components while the bulk response remains elastic so that

Equation (16) resembles the Maxwell chain rheological model with $M$ Maxwell units, where $G_\mu$ is the shear modulus of the elastic spring of the $\mu$-th unit and $e_{i j}^{p, \mu}$ is the deviatoric creep strain developed in the $\mu$-th dashpot in accord with the Eyring flow model

where $\eta_0^\mu$ is the zero shear viscosity, $a_\sigma$ is the stress-dependent shift factor to address the dependence of the response on the current stress stress level and thus also on the applied strain rate, and $\tau_0$ is the model parameter determined experimentally [ $14,21,23$ ].

In the present study, we adopt the data pertinent to 285/500 aero Havel epoxy resin and derived in [21] from a series of creep tests performed at different stress levels and a series of tensile tests carried out at different strain rates. The latter ones allowed us to determine $\tau_0$ while the former ones served to construct the compliance master curve. It has been shown that this experimentally obtained master curve can be well approximated by 10 units of the Dirichlet series. The resulting compliances $J_\mu$ for a-priory selected retardation times $\tau_\mu$ are listed in Table 1. The Laplace transform was then employed to get the corresponding pairs for relaxation times $\theta_\mu$ and stiffnesses $E_\mu$ to describe the relaxation function more suitable for numerical simulations. The shear moduli $G_\mu$ in Eq. (16) ${ }_1$ were found from $E_\mu$ for the selected

Poisson ratio in the matrix $v_m=0.39$. The same Poisson ratio was used to derive the bulk modulus $K=\frac{\sum_\mu E_\mu}{3\left(1-2 v_m\right)}$ in Eq. (16) 2 .

The elastic properties of the examined fibers are stored in Table 2 together with the corresponding volume fractions extracted from the images of real microstructures.

The computational part begins with the calibration of parameters M, N, T of the damage model in Eq. (15). Once these are known, the MT method can be exploited as a stress updater at the level of yarn in the multiscale analysis of plain weave textile composite, which is the principle objective of this contribution. These computational steps are described next.

|The model parameters $M, N, T$ are found by comparing the mesoscopic in-plane shear response of the fiber tow predicted by the detailed finite element simulations and the MT method. In the present study, we adopt a certain simplification and approximate the actual microstructure by the periodic hexagonal array (PHA) model in Fig. 2. Therefore, the two fiber systems differ only geometrically by the volume of the fiber phase. Note that the local $z$-axis is aligned with the fiber direction.

To appreciate the need for a suitable modification of the MT method, we loaded the composite by the prescribed mesoscopic shear strain rate $2 \dot{\bar{\varepsilon}}_{x y}=\dot{\bar{\gamma}}_{x y}=10^{-4} \mathrm{~s}^{-1}$. The corresponding mesoscopic strain-stress diagrams for the two composite systems are plotted in Fig. 1. We clearly see that the original MT formulation with $\omega=0$ significantly overestimates the predictions obtained from PHA simulations. On the contrary, introducing the gradually evolving damage parameter $\omega$ with properly adjusted parameters $N, M, T$ provides almost a perfect match.

An obvious step forward is to check whether fitting the model parameters to shear loading only provides satisfactory results for other types of loading conditions. To that end, we loaded the system in tension by the prescribed tensile strain rate $\dot{\bar{\varepsilon}}_{x x}=10^{-4} \mathrm{~s}^{-1}$ in the direction normal to the fiber. The results appear in Fig. 2(c). While a certain improvement is observed for the basalt fiber-based system, the response obtained for carbon fibers is slightly more compliant than that of PHA model. This may suggest a more general calibration step combining the data from several loading directions and strain rates. While this issue is being currently examined, we accept the values of $N, M, T$ obtained purely from shear loading as satisfactory for further study aimed at a macroscopic behavior of textiles.

It is worth mentioning that all simulations assumed strain control loading conditions. This means that we prescribed all strain components of the mesoscopic strain rate vector $\overline{\bar{\varepsilon}}$, i.e. $\dot{\bar{\varepsilon}}^{\top}=\left\{\dot{\bar{\varepsilon}}_{x x}, 0,0,0\right\}$ in case of transverse tension along the $x$-axis and plane-strain state of stress. Similar loading scenarios are also considered in the next section devoted to the ply level where we substitute the mesoscopic strains $\bar{\varepsilon}$ by their macroscopic counterparts $E$.

Although typically supplied in terms of laminates, we limit our attention in this preliminary study to a single ply composite. The computational model for carbon-fiber-based composite is shown in Fig. 3 identifying basic geometrical data (Fig. 3(a)) to construct the most simple RVE, a periodic unit cell (PUC) of a plain weave textile ply (Fig. 3(b)). Note that the two types of composite systems differ again in their geometrical details only. These are available in Table 3.

More complex geometrical models taking into account various types of imperfections including voids are also available [14, 24], but this goes beyond the present scope.

The multiscale modeling strategy adopted in the present study is evident in Fig. 4. The finite element mesh of basalt-fiber and carbon-fiber based composite consisted of 51087 and 84448 constant strain tetrahedral elements, respectively. Similarly to $\mathrm{FE}^2$ scheme, each element of the yarn was considered as a two-phase unidirectional fibrous composite loaded by the increment of a mesoscopic strain $\Delta \bar{\varepsilon}$ as seen in Fig. 4. But instead of solving this sub-problem by FEM exploiting another PUC, for example, the one in Fig. 2(a), we adopt the MT averaging scheme to provide the piecewise uniform averages of local fields and their corresponding mesoscopic counterparts.

For illustration, we compare the macroscopic response of both material systems to in-plane and out-of-plane shear and in-plane tension. To that end, we loaded the unit cells in turn by the macroscopic in-plane $2 \dot{E}_{x y}=10^{-4} \mathrm{~s}^{-1}$ and out-of-plane $2 \dot{E}_{x z}=10^{-4} \mathrm{~s}^{-1}$ shear strain rates and by the macroscopic tensile strain rate $\dot{E}_{x x}=10^{-4} \mathrm{~s}^{-1}$. The results are plotted in Fig. 5.

To confirm the fact that viscoelastic response is driven by the deviatoric components of the stress field, we plot in Fig. 6(a) the evolution of the mean $\Sigma_m=\frac{1}{3} \Sigma_{i i}$ and equivalent deviatoric $J=\sqrt{\frac{1}{2} S_{i j} S_{i j}}, \quad S_{i j}=\Sigma_{i j}-\Sigma_m \sigma_{i j}$ stresses clearly identifying the elastic bulk response. A rather small deviatoric stresses also explain a more or less elastic response of the textile ply when loaded in tension, recall Fig. 5(b).

When using commercial codes, the direct derivation of shear response in Fig. 5(a) by applying the macroscopic shear strain might not be possible. In that case, the in-plane shear stresses might be derived by loading the composite by combined tension and compression, for example, by setting $\dot{E}_{x x}=-\dot{E}_{y y}=2 \dot{E}_{x y} / 2$. Given the loading condition to construct $2 E_{x y} \times \Sigma_{x y}$ curve in Fig. 5(a), we thus set $\dot{E}_{x x}=-\dot{E}_{y y}=5 \times 10^5$. The results are plotted in Fig. 6(b). An almost perfect agreement with $2 E_{x y}$ loading has been achieved. This is attributed here to a perfect symmetry in the weft and warp directions resulting in zero out-of-plane stresses.

However, a word of caution is needed. When considering the same exercise with PHA model on meso-scale, it is seen in Fig. 7 that the relation $\bar{\sigma}_{x y}=\frac{1}{2}\left(\sigma_{x x}-\sigma_{y y}\right)$ breaks at the onset of nonlinear response. This is because of evolution of axial normal stress $\bar{\sigma}_{z z}$ due to viscoelastic strain in the matrix and the mesoscopic constraint $\bar{\varepsilon}_{z z}=0$. The transverse isotropy valid for elasticity is then lost.

The paper described an application of the MT averaging scheme in multiscale analysis of the viscoelastic response of textile composites. It has been shown that such a relatively complex task can be solved very efficiently, which would not be the case when employing the standard FE2 computational approach, unless a massive parallel computation is used. Nevertheless, there are still some open questions including more complex optimization procedure to tune the parameters of the proposed damage model and to validate the applicability of the MT method experimentally. Both issues are currently under investigation.

The financial support provided by the GAČR grant No. 19-15666S and by the Czech Technical University in Prague within SGS project with the application registered under the No. OHK1-026/21 gratefully acknowledged.