# An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite

Stéphane André, Julien Boisse, Camille Noûs

► To cite this version:

Stéphane André, Julien Boisse, Camille Noûs. An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite. Journal of Theoretical, Computational and Applied Mechanics, 2021, pp.1-19. 10.46298/jtcam.6450. hal-02562362v3

HAL Id: hal-02562362

https://hal.science/hal-02562362v3

Submitted on 28 Apr 2021

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
**Identifiers**

DOI 10.46298/jtcam.6450

OAI hal-02562362v3

**History**

Received May 6, 2020

Accepted Feb 12, 2021

Published Apr 16, 2021

**Associate Editor**

Laurence Brassart

**Reviewers**

François Willot

Anonymous

**Open Review**

DOI 10.5281/zenodo.4696158

**Supplementary Material**

Software
[SoftwareHeritage]

**Licence**

CC BY 4.0
©The Authors

# An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite

©Stéphane André<fnref n="1" />, ©Julien Boisse<fnref n="1" />, and ©Camille Noûs<fnref n="2" />

FFT-based solvers are increasingly used by many researcher groups interested in modelling the mechanical behavior associated to a heterogeneous microstructure. A development is reported here that concerns the viscoelastic behavior of composite structures generally studied experimentally through Dynamic Mechanical Analysis (DMA). A parallelized computation code developed with complex-valued quantities provides virtual DMA experiments directly in the frequency domain on a heterogeneous system described by a voxel grid of mechanical properties. The achieved precision and computation times are very good. An effort has been made to illustrate the application of such a virtual DMA tool through two examples from the literature: the modelling of glassy/amorphous systems at a small scale and the modelling of experimental data obtained in temperature sweeping mode by DMA on a particulate composite made of glass beads and a polystyrene matrix, at a larger scale. Both examples show how virtual DMA can contribute to question, analyze, and understand relaxation phenomena on either theoretical or experimental points of view.

**Keywords:** Dynamic Mechanical Analysis, FFT solver, viscoelasticity, homogenization, glassy polymers, particle composites

## 1 Introduction

Dynamic Mechanical Analysis (DMA) is known to be a privileged tool to study materials (especially polymers and rubbers) whose rheological behavior is viscoelastic by nature, i.e. introduces irreversible dissipation of mechanical energy into heat. The technique indeed measures a conservative (storage) or dissipative (loss) modulus or compliance, which are the real and imaginary parts of their complex nature $M^*(\omega) = M'(\omega) + jM''(\omega)$ for instance for the modulus. It relies on harmonic steady-state excitations in strain, applied at pulsation $\omega$ on a material specimen and on the recording of the corresponding output stress signal. This latter is analyzed with respect to the original input both in amplitude decrease (modulus $|M^*|$) and phase lag $\delta$ (reflecting the damping phenomenon associated to some viscous component of the behavior). From that information, complex algebra gives eventually access to $M'(\omega)$ and $M''(\omega)$ with a sweep in frequency allowing for full dynamical characterization of the material. Such a spectroscopic probing of the material exists in many different fields of physics like in thermal science to produce thermal conductivity/diffusivity measurements (Cahill 1990) but the main corpus of publications probably resorts to dielectric properties measurements (Nielsen and Landel 1994). Experimental data on the susceptibility of physical processes as function of frequency is, in general, a central tool to develop physical model of relaxation processes (Havriliak Jr. and Havriliak 1995; Jonscher 1996).

DMA is estimated 100 times more sensitive to glass transition than scanning calorimetry (Menard and Menard 2015) and this is probably the reason why the sweeping temperature mode is generally elected by material scientists as the perfect way of investigating subtle

<footnote n="1">¹ Université de Lorraine, CNRS, LEMTA, F-54000 Nancy, France</footnote>

<footnote n="2">² Cogitamus Laboratory, F-75005 Paris, France</footnote>
microstructural transitions like the motions of polymer chains for example (Seidel 2008). However, for solid-liquid rheological characterization, the forced frequency sweep is more desirable as it allows behavior's law modeling assumptions to be checked and associated material parameters to be identified. The targeted information is obviously the relaxation spectrum which can be precisely investigated with the frequency scanning. Emri and Tschoegl (1995) synthesize the spectrum from temporal data in a collocation-like approach but check the consistency with DMA results in frequency mode. In (Kim and Lee 2009) on the contrary, experimental characterization of Frequency Response Functions (FRF) are used to identify parameters originating from constitutive rheological models of damping materials. One key issue of this study precisely results in the solution offered by virtual DMA to provide quickly these data for any kind of material and to test theoretical concepts underlying such models.

It should be pointed that this technique is before anything else a perfect tool to investigate the behavior of a material at the scale of a Representative Volume Element (RVE). Samples are generally of small (macroscopic) sizes and excitations of low magnitude (small perturbations theoretical framework). Independently of excitation modes and specimen geometry (torsion, flexure, compression), this allows retrieving local material parameters directly from rheological models. These parameters are the exact reflect of statistically averaged microstructural evolutions.

The concern of this paper is classical in scientific calculus: offering a simulation path to replace experiments and exceed their intrinsic limitations or, in other words, to offer a virtual DMA simulator capable of analyzing any kind of heterogeneous or composite material. Such simulation tool has been already striven after in the past through Finite element approach. Brinson and Knauss (1992) for example have modified a Finite Element code to make the computation of real and imaginary complex moduli directly possible, solving the boundary value problem with complex variables.

New computational tools, more efficient, have made recently possible to rejuvenate this idea thanks to two important advances in computational science:

1. Synthesis of virtual composite microstructure using generators working on various mathematical basis (Ghazvinian et al. 2014; Quey 2021; Salnikov et al. 2015) or alternatively synthesis of real microstructure based on tomography imaging followed by appropriate image treatment (Uchic et al. 2007). Virtual DMA performed with FE codes was initially limited to very simple idealized microstructures.

2. Development of the spectral approach (Fast-Fourier Transform operator) to handle the resolution of local equilibrium equations, see (Moulinec and Suquet 1998; Roters et al. 2019) for a wide review on the topic, especially for works developed using the Damask code. Of course, the harmonic steady-state regime can be obtained from time domain simulations based for example on Finite Element approach (Masurel et al. 2015), but such approaches are very much less efficient than the spectral ones directly operating in frequency domain.

The idea is before anything else to make the confrontation of both approaches (simulation and experiment) a source of better knowledge arising from their respective drawbacks.

1. Simulation will always be limited by the idealization of the composite organization: for example, if one can faithfully and easily reproduce particulate composite by considering real size distributions of particulate, real volume fractions..., continuity assumptions are generally considered between particles and the matrix. Overcoming this complexity is gaining attention however. An homogenization problem in thermal science has recently addressed the case of imperfect interfaces within the FFT spectral approach (Monchiet 2018).

2. Experiments based on DMA can carry a bias inherent to the technique: as an example, one can cite the drift in measurement signals with the very long-time durations of these experiments to get low frequency information.

Of course, the second key idea is the homogenization problem of composite materials. It is in that direction that very recent computations of virtual DMA kind have been performed with FFT solvers. In a probably landmark article, Figliuzzi et al. (2016) studied the composite material made of a rubber matrix filled with carbon black fillers. The effective mechanical behavior resulting from various models considered for the multiscale morphological microstructure description has
been successfully computed in 3D with an FFT solver in space and harmonic complex treatment of the dynamic part. The frequency dependence of the effective complex bulk modulus and of the effective shear modulus was directly achieved and shown as a valuable way to assess the performance of analytical effective models. A second study of that type was used by Gallican and Brenner (2019) to investigate the overall properties of composite materials with fractional viscoelastic constituents. In the case of particulate composites with polydisperse spherical elastic inclusions, the authors develop FFT based computations of the overall complex moduli in frequency domain. Such outputs were helpful to assess mean-field approximate models and a generalized fractional effective model derived from exact mathematical (asymptotic) relations constraining the adjustment of an effective relaxation spectrum.

It is worthwhile to mention that whatever the virtual DMA computation tool (FEM or FFT-based solvers), a comfortable aspect is that viscoelastic theories in the linear regime require only infinitesimal deformation analysis for strains and displacements.

The results presented here enter into this movement and favors an account of the benefits that can be expected from virtual DMA to explain physical mechanisms or question theoretical model assumptions as well as experimental data. In Section 2 will be firstly described the extension to the harmonic steady-state regime of spectral FFT solvers used to solve local equilibrium equations generally in the steady state or for temporal responses. Validation of the numerical implementation will consider the academic heterogeneous checkerboard structure with constituents of Standard Linear Solid (SLS) rheological behavior. Additionally, it will illustrate how outputs of virtual DMA can be used to identify material transfer functions, a concept which will allow to compute the response of the material to any kind of solicitation. Section 3 will be devoted to illustrate what a virtual DMA solver can bring to material engineering science in the future. Two test cases will be considered, each of them showing a different aspect of the subject. Example 1, in Section 3.1, will compare the results based on the paper by Masurel et al. (2015) for a multi-material made of a collection of Maxwell units and obtained by time-domain calculations using the FE code Zébulon (Armines & ONERA 2021). Finally, example 2, in Section 3.2, will illustrate the benefit of using virtual DMA when trying to interpret experimental results. It will rely on data relative to strong experimental works made by Alberola and collaborators (Agbossou et al. 1993; Alberola and Mele 1996) on a Particulate Composite of glass beads in a polystyrene matrix. It shows also that harmonic effective behavior can be calculated directly from the knowledge of either the temperature or frequency-discretized moduli of the constituents. This option can be very useful for the experimentalist where DMA data are often produced at a given frequency with a sweep in temperature.

# 2 A FFT Spectral solver extension to harmonic regime

## 2.1 Standard mathematical procedure and numerical implementation

We consider a heterogeneous multi-material system, with a Representative Volume Element V, subject to oscillatory boundary conditions with angular frequency $\omega$. The heterogeneous harmonic displacement field $\tilde{u}(\mathbf{r}, t)$ in such harmonic regime can be expressed with a complex form:

$$ \tilde{u}(\mathbf{r}, t) = \mathbf{u}^*(\mathbf{r})e^{j\omega t} \qquad (1) $$

where $\mathbf{u}^*(\mathbf{r})$ denotes the complex displacement amplitude, whose real part corresponds to the physical displacement field. Considering the harmonic regime, and in the small perturbation approximation, the spectral method based on the Lippman-Schwinger equation associated to a discrete Green operator is straightforward. The complex amplitudes of the strain $\varepsilon^*(\mathbf{r})$ and stress $\sigma^*(\mathbf{r})$ are second rank tensors, which depend on the frequency $\omega$.

Now considering a heterogeneous medium made of linear viscoelastic phases and thanks to the correspondence principle (Hashin 1970), local constitutive laws expressed below in terms of relaxation functions can be put directly in the complex form (Laplace-Carson transforms in time where the Laplace variable $p$ is set to $j\omega$), with a linear relation:

$$ \sigma^*(\mathbf{r}, j\omega) = \mathbf{H}^*(\mathbf{r}, j\omega) : \varepsilon^*(\mathbf{r}, j\omega) \qquad (2) $$
The complex local stiffness tensor $\mathbf{H}^*(\mathbf{r}, j\omega)$ accounts for the storage and loss moduli.

In the frequency (spectral) domain, the local equations expressing the constitutive law in Equation (2), the mechanical equilibrium, and the overall prescribed strain loading produce the system

$$ \left\{ \begin{array}{ll} \sigma^*(\mathbf{r}, j\omega) = \mathbf{H}^*(\mathbf{r}, j\omega) : \boldsymbol{\varepsilon}^*(\mathbf{r}, j\omega) & \forall \mathbf{r} \in V, \\ \mathbf{div} \, \boldsymbol{\sigma}^*(\mathbf{r}, j\omega) = \mathbf{0} & \omega \in [0; +\infty[, \\ \langle \boldsymbol{\varepsilon}^*(j\omega) \rangle = \bar{\boldsymbol{\varepsilon}}^*. \end{array} \right. \quad (3) $$

It can be solved by making use of a numerical scheme based on the Fourier Transform (Figliuzzi et al. 2016; Moulinec and Suquet 1998) considering the RVE with periodic boundary conditions. The elastic classical scheme introduced by Moulinec and Suquet has to be transposed in the Laplace-Carson domain with complex-valued quantities. The Lippman-Schwinger equation is transformed in the Fourier domain ($\hat{\bullet}$ symbol) and depends now on spatial frequencies $\mathbf{k}$ and dynamical testing frequencies $\omega$:

$$ \hat{\boldsymbol{\varepsilon}}^*(\mathbf{k}, j\omega) = -\hat{\Gamma}_0^*(\mathbf{k}) : \hat{\boldsymbol{\tau}}^*(\mathbf{k}, j\omega), \quad \forall \mathbf{k} \neq \mathbf{0}, \quad \hat{\boldsymbol{\varepsilon}}^*(\mathbf{0}, j\omega) = \bar{\boldsymbol{\varepsilon}}^*, \quad \omega \in [0; +\infty[ \quad (4) $$

The periodic Green's operator in Fourier domain is denoted $\hat{\Gamma}_0^*$ and is associated to the reference medium of isotropic stiffness tensor $\mathbf{H}_0^*$, which has to be real (Figliuzzi et al. 2016); $\hat{\boldsymbol{\tau}}^*(\mathbf{k}, j\omega)$ is the Fourier transform of the polarization field $\boldsymbol{\tau}^*(\mathbf{r}, j\omega) = \boldsymbol{\sigma}^*(\mathbf{r}, j\omega) - \mathbf{H}_0^* : \boldsymbol{\varepsilon}^*(\mathbf{r}, j\omega)$; $\bar{\boldsymbol{\varepsilon}}^*$ denotes the average over the whole RVE (Fourier frequency $\mathbf{k} = \mathbf{0}$) of the local complex strain field $\boldsymbol{\varepsilon}^*(\mathbf{r}, j\omega)$. Equation (4) can be solved in Fourier domain (spatial discretization) and for the harmonic regime described by any angular-frequency set in $[0; +\infty[$ and with fully complex local $\boldsymbol{\varepsilon}^*(\mathbf{r}, j\omega)$ and $\boldsymbol{\sigma}^*(\mathbf{r}, j\omega)$ fields. A numerical scheme is necessary, of fixed-point type, to obtain convergence along an iterative process over the polarization field variable, where iterate $n + 1$ is updated from the previous iterate $n$. High-performance accelerated numerical schemes are available for this purpose (Eyre and Milton 1999; Moulinec and Silva 2014).

A convergence test is computed at each iteration. It consists in checking the deviation from equilibrium $e_{eq}$, from compatibility $e_{comp}$ and from the prescribed loading conditions $e_{load}$. The iterative procedure is stopped when these three criteria are smaller than some prescribed values (most of the times fixed to $10^{-4}$).

The deviation from equilibrium is obtained with the following dimensionless criterion, calculated in the Fourier space:

$$ e_{eq} = \frac{\| \mathbf{div} \, \boldsymbol{\sigma}^* \, \mathrm{d}V \|}{\| \boldsymbol{\sigma}^* \cdot \mathbf{n} \, \mathrm{d}S_n \|} = \frac{\sqrt{\sum_k \hat{\mathbf{f}}_{eq} \cdot \hat{\mathbf{f}}'_{eq}}}{\sqrt{\sum_k \{\hat{\mathbf{f}}_1 \cdot \hat{\mathbf{f}}'_1 + \hat{\mathbf{f}}_2 \cdot \hat{\mathbf{f}}'_2 + \hat{\mathbf{f}}_3 \cdot \hat{\mathbf{f}}'_3\}}} \quad (5) $$

where the prime symbol refers to a complex conjugate and with the forces

$$ \begin{aligned} \hat{\mathbf{f}}_{eq}(\mathbf{k}) &= \hat{\boldsymbol{\sigma}}^*(\mathbf{k}) : j \mathbf{k} \, \mathrm{d}V, & \hat{\mathbf{f}}_1(\mathbf{k}) &= \hat{\boldsymbol{\sigma}}^*(\mathbf{k}) : \mathbf{n}_1 \, \mathrm{d}x_2 \, \mathrm{d}x_3, \\ \hat{\mathbf{f}}_2(\mathbf{k}) &= \hat{\boldsymbol{\sigma}}^*(\mathbf{k}) : \mathbf{n}_2 \, \mathrm{d}x_1 \, \mathrm{d}x_3, & \hat{\mathbf{f}}_3(\mathbf{k}) &= \hat{\boldsymbol{\sigma}}^*(\mathbf{k}) : \mathbf{n}_3 \, \mathrm{d}x_2 \, \mathrm{d}x_1. \end{aligned} \quad (6) $$

The deviation from compatibility is obtained with the following expression calculated in the Fourier space

$$ e_{comp} = \frac{\max_k(\max_{i=1,...,6}(|\hat{e}_i(\mathbf{k})|)) \cdot \max(L_1, L_2, L_3)^2}{\sqrt{\sum_k \hat{\varepsilon}_{mn}^*(\mathbf{k}) \cdot \hat{\varepsilon}_{mn}'(\mathbf{k})}} \quad (7) $$

where $L_i = N_i \, \mathrm{d}x_i$ is the size of the box in directions $i = 1, 2, 3$ and with $\hat{e}_i(\mathbf{k})$ expressed for the six compatibility equations in terms of the tensorial components $\hat{\varepsilon}_{mn}^*$.

In the case of a prescribed macroscopic strain $\bar{\boldsymbol{\varepsilon}}_{mn}^*$, the deviation from the prescribed loading is given by

$$ e_{load} = \frac{\|\langle\boldsymbol{\varepsilon}^*\rangle - \bar{\boldsymbol{\varepsilon}}^*\|}{\|\bar{\boldsymbol{\varepsilon}}^*\|} = \frac{\sqrt{\langle\langle\boldsymbol{\varepsilon}_{mn}^*\rangle - \bar{\boldsymbol{\varepsilon}}_{mn}^*\rangle : \langle\langle\boldsymbol{\varepsilon}_{mn}^*\rangle - \bar{\boldsymbol{\varepsilon}}_{mn}^*\rangle'}}{\sqrt{\bar{\boldsymbol{\varepsilon}}_{mn}^* \bar{\boldsymbol{\varepsilon}}_{mn}^*}} \quad (8) $$
Other loading conditions, like prescribed stress or multi-directional loadings, can be treated in a similar way. In the Eyre and Milton scheme, and considering a multi-material system with *M* linear viscoelastic phases the optimal choice for bulk *K* and shear *μ* moduli for the isotropic reference medium is

$$ 
\begin{aligned}
\mu_0 &= \sqrt{\max_{i,\omega}(\mathfrak{R}[\mu_i(j\omega)]) \cdot \min_{k,\omega}(\mathfrak{R}[\mu_k(j\omega)])} \\
K_0 &= \sqrt{\max_{i,\omega}(\mathfrak{R}[K_i(j\omega)]) \cdot \min_{k,\omega}(\mathfrak{R}[K_k(j\omega)])}
\end{aligned}
\qquad (9) $$

for $i \neq k$, $i, k \in 1, \dots, M$, $\omega > 0$ and $\omega \in [\omega_{min}; \omega_{max}]$. This mathematical procedure has been implemented as a new option in a spectral solver software named CraFT (Moulinec et al. 2004) for Composite response and Fourier Transform developed at LMA laboratory by several scientists since the initial version raised by Moulinec and Suquet. This option takes advantage of all recent developments of the solver in terms of accelerated convergence schemes and is available in parallelized configuration for enhanced computation times. The harmonic version is simply available by selecting an appropriate compiling option which proceeds with a set of functions that have been modified to handle the complex nature of the variables (65 functions over a total of 263 including the behavior law implementation and solver operations). The main change with respect to the temporal version lies in the specification of frequency vectors in the Fourier domain due to the two-sided character of the Discrete Fourier Transform. The code modifications render it even more regular in terms of tuplet management. Indeed, for a real function, say a stress tensor component, the 3D double array of dimension $n_0 \times n_1 \times n_2$ (indices in the three directions determining which voxel is considered) produces a DFT output of $(n_0/2 + 1) \times n_1 \times n_2$ complex numbers<fnref n="1" />. In the harmonic version, we have to handle 3D data arrays of complex numbers which restore a simpler approach of the code. Additionally, the Fourier domain equilibrium criterion in Equation (5) has been implemented, which is formulated in terms of a true normalized force balance (volumetric over surface for each voxel). In comparison with the original CraFT criterion based on the stress tensor divergence presenting a dimensional dependence, this results in a unique criterion specification.

Spectral solvers imposed themselves as an alternative approach for computing mechanical problems on heterogeneous structures with the following recognized advantages: precision and very fast computational times. Note that for stochastic materials, the latter is a crucial one as multiple simulations are required for identical statistical figures to obtain an averaged meaningful response of the heterogeneous RVE. Among the drawbacks, one can cite two of them, being already subjected to a common proposal for reducing their nuisance: (i) spectral solvers can produce high frequency oscillations on the mechanical fields close to interfaces; (ii) because the voxel distribution is fixed once (no remeshing), the precision of the method can deteriorate rapidly in the case of complicated interfaces and/or multiple length scales. For both drawbacks, the composite voxel approach offers substantial improvements: a filtering effect for the first drawback (Gélébart and Ouaki 2015) and more accurate predictions in the latter case (Charière et al. 2020; Mareau and Robert 2017).

## 2.2 Numerical Validation: the checkerboard case

We consider the academic case of the checkerboard composite, classically used for validating the FFT simulations (Lebensohn et al. 2005). The microstructure is made of a periodic 2D, two-phase composite whose unit cell consists of four tiles with the crystallographic orientations of the two pairs of opposite grains being at +90° and 0°. Each material has a rheological behavior described by the three-parameters SLS solid

$$ \frac{\bar{\sigma}(p)}{\bar{\epsilon}(p)}|_{p=j\omega} = G^{(k)}(j\omega) = \frac{G_r^{(k)} + G_g^{(k)} j\omega \tau^{(k)}}{1 + j\omega \tau^{(k)}}, \quad k = 1, 2 \qquad (10) $$

where $G_g$, $G_r$, and $\tau$ denote the glass (instantaneous) modulus, the relaxed modulus and a relaxation time, respectively. The effective behavior of this checkerboard structure with two

<footnote n="1">¹ See the link [DFT-Array].</footnote>
constituents and two orientations submitted to an overall shear excitation has an analytical solution. The Hashin principle states that

$$ \tilde{G}(j\omega) = \sqrt{G^{(1)}(j\omega)G^{(2)}(j\omega)} \qquad (11) $$

which implies

$$ \tilde{G}(j\omega) = \tilde{G}_g \sqrt{\frac{(j\omega + 1/\tau_z^{(1)})(j\omega + 1/\tau_z^{(2)})}{(j\omega + 1/\tau^{(1)})(j\omega + 1/\tau^{(2)})}} \qquad (12) $$

with

$$ \tilde{G}_g = \sqrt{G_g^{(1)} G_g^{(2)}}, \quad \tilde{G}_r = \sqrt{G_r^{(1)} G_r^{(2)}} \quad \text{and} \quad 1/\tau_z^{(k)} = G_r^{(k)}/(G_g^{(k)} \tau^{(k)}) $$

where $\tau_z^{(k)}$ is the relaxation time which appears in the zeros of the rational function in Equation (12).

Real and imaginary parts of the effective complex modulus $\tilde{G}$ can be calculated exactly from Equation (12) and directly compared to the numerical FFT solver outputs in virtual DMA mode. Figure 1 plots the results in the case of a strong contrast in mechanical properties—same values as those reported by Suquet (2012). Table 1 synthesizes the input and output data.

<figure><img src="image_2.png" /><figcaption>**Figure 1** Validation of the CraFT virtual DMA implementation. Comparison to the (a) exact effective real (storage) and (b) imaginary (loss) moduli: $\tilde{G}(j\omega)$ defined in Equation (12) and $Y(j\omega)$, in Equation (13).</figcaption></figure>

<table><caption>**Table 1** Inputs/Outputs of the checkerboard test case. Units: relaxation times in second and moduli in MPa.</caption><thead><tr><th></th><th>Instantaneous (glassy) modulus</th><th>Relaxed modulus</th><th>Relaxation times</th></tr></thead><tbody><tr><td>Constituent 1</td><td>10</td><td>9</td><td>5</td></tr><tr><td>Constituent 2</td><td>100</td><td>50</td><td>0.05</td></tr><tr><td>Effective exact <span style="text-decoration: overline;">G</span></td><td><span style="text-decoration: overline;">G</span><sub>g</sub> = &radic;1000</td><td><span style="text-decoration: overline;">G</span><sub>r</sub> = &radic;450</td><td>&tau;<sup>(1)</sup> = 5</td><td>&tau;<sup>(2)</sup> = 0.05</td><td>&tau;<sub>z</sub><sup>(1)</sup> = 50/9</td><td>&tau;<sub>z</sub><sup>(2)</sup> = 0.1</td></tr><tr><td>Effective identified model</td><td><span style="text-decoration: overline;">G</span><sub>g</sub> = 31.62</td><td><span style="text-decoration: overline;">G</span><sub>r</sub> = 21.21</td><td>&tau;&#770;<sub>p1</sub> = 4.97</td><td>&tau;&#770;<sub>p2</sub> = 0.059</td><td>&tau;&#770;<sub>z0</sub> = 4.42</td><td></td></tr></tbody></table>

It is clear from Figure 1 that the virtual DMA computations are in perfect agreement with the exact solution. One of the virtual DMA advantages can already be illustrated here: inverse parameter identification (i.e. metrology of material parameters) can be applied to virtual DMA outputs, providing an effective model has been selected. From theoretical considerations (Gallican et al. 2017), it is known that the effective material built from two SLS-type constituents does not follow this rheological behavior itself. This is clearly evidenced in Figure 1. The two-order-of-magnitude ratio between $\tau^{(1)}$ and $\tau^{(2)}$ implies a marked spectrum of the dissipative (loss) modulus. An effective model could be proposed through the collocation method applied with a two-term Prony series. In Laplace domain, this leads to a rational admittance $Y(j\omega)$ having two poles (the relaxation times) and one zero:

$$ Y(j\omega) = G_g + (G_r - G_g) \frac{(1 + j\omega\tau_{z0})}{(1 + j\omega\tau_{p1})(1 + j\omega\tau_{p2})} \qquad (13) $$
Table 1 gives the virtual measurements of the material parameters of this model as resulting from a least-square minimization of the residuals between a (here virtual) DMA experiment and the candidate effective model, see Equation (13). The agreement is perfect which means that the mathematical structure of the collocation model (Biot model) catches all the physical content of the signal and would be able to predict the response of the composite material to any input excitation. The moduli $\tilde{G}_g$ and $\tilde{G}_r$ are perfectly recovered and the relaxation times (the two poles of the admittance fraction in this case) are shown to be very close to those of the constitutive materials.

# 3 Virtual DMA Applications and Interest

## 3.1 Application 1: Heterogeneous dynamics in glass/amorphous polymers

### 3.1.1 The heterogeneous problem

The composite material of interest in this section was investigated by Masurel et al. (2015) as an academic approach to understand dynamical heterogeneities which take place in amorphous polymers in the glass transition. We named it a *multi-material composite* as it is built from a grid of elements with individual behavior being of Standard three-parameter Maxwell rheological type, see Figure 2(a). It can be considered as a kind of construction of an *averaged* object representing the entire system at all sub-scales. Observed relaxation processes in disordered systems reveal a kind of universal nature, with now well identified features and/or explicative ideas (cooperative dynamics, cluster formations, energy-criterion arguments...). This multi-material composite aims at catching some of these general features from basic ingredients that are used in a composite view. This is fully motivated by experimental evidences that non-exponential responses of a macroscopic sample results for example from a relaxation-rate distribution of independently randomly sized mesoscopic regions (Jurlewicz and Weron 2002). Here the domains are of identical

<figure><img src="image_3.png" /><figcaption>**Figure 2** Heterogeneous multi-material definition with a total virtual DMA simulation over 25600 voxels. (a) Heterogeneous multi-material of 40 × 40 domains. Each domain *i* has a Standard Linear Solid Maxwell rheological behavior with uniform glassy and rubber moduli over the RVE. Randomly distributed relaxation times $\tau_i$. (b) Zoom over the multi-material: each domain is made of 16 voxels.</figcaption></figure>

sizes but dynamics on each domain are randomly distributed. A few identical elements can be associated to form a domain *i* having a constant relaxation time $\tau_i$ (same mechanical response) and domains are associated randomly to build the polymer *coarse-grained* model. The spectrum of relaxation times has been chosen according to a normal distribution of width *w* (in s) with maximum centered on $\tau_{max} = 1$ s and in logarithmic scale. The glassy and rubber elastic constants ($G_g$ and $G_r$ respectively), being identical for all domain, they applied then to the homogenized composite. This is a strong assumption of the approach which has the advantage to isolate the only contribution of heterodynamics to the homogenized relaxation spectrum (if any).

The objective of the study was to observe the effective behavior of such model in terms of effective relaxation spectrum as well as the local mechanical fields interplay. To this aim, FE simulations were performed using Zébulon code in 2D plane strain condition (frame axis 1, 2),
assuming incompressibility with a bulk modulus set to $K = 10^5 G_g$, periodic boundary conditions, and a pure shearing test (imposed macroscopic shear step strain $\epsilon_{12} = 0.01$).

We must point here that the study produces *time domain* simulation and *therefore does not give a direct way of computing the complex, frequency dependent, modulus $G^*(\omega)$*. Masurel et al. (2015) had to post-process the simulated data in order to identify the frequency complex modulus from the time response of the macroscopic $\langle\sigma_{12}(t)\rangle$.

### 3.1.2 Spectral solver direct solution compared to other theoretical published results

Craft solver in harmonic and parallelized version was used to perform calculations directly in frequency domain (virtual DMA). The exact same inputs have been used as those given by Masurel et al. (2015) and Figure 2(a) shows an example of the random relaxation times assignment over the voxels used in the FFT solver. In Figure 3 are presented the real and imaginary parts of the complex modulus for various cases of the ratio $G_g/G_r$ as published by Masurel et al. (2015) (dotted curves) and the present results obtained with a 3D spectral solver (solid lines). The results

<figure><img src="image_4.png" /><figcaption>**Figure 3** CraFT virtual DMA simulations of heterogeneous dynamics. (a) Storage and (b) dissipative moduli versus frequency for various contrasts $G_g/G_r$ in the glassy and relaxed moduli. Data from (Masurel et al. 2015).</figcaption></figure>

are very close also some differences can be observed for large $G_g/G_r$ ratio which constitutes the most severe case to simulate. It is not useful here to discuss the origin of such discrepancies because they can be numerous. It can just be underlined that:

1. The spectral solver is much more adapted to the multi-material case presented here than FE simulations as the spectral solver directly proceeds on the cubic (or squared) voxels. There is no error introduced here through some meshing of FE type (edge effects, location of quadrature nodes...) as the geometry of the problem is perfectly respected into the FFT scheme. This means that no errors can be made due to the spatial resolution considered. In the FE approach, on the contrary, the authors had to perform a preliminary sensitivity analysis to obtain a representative multi-material which preserves accessible computation times. No such insurance is required in the FFT approach.

2. As already said, the procedure to obtain the complex modulus is direct with the FFT solver and indirect in the FE approach, which must be coupled to the identification of a discretized weighted function.

Some errors are introduced in each of this feature specific to the FE approach and precludes any further investigation at this stage. Note that 3D simulations were probably beyond reach with FE approach whereas they can be done in very reasonable CPU times with the FFT solver. For the same test case, a comparison is also given in Appendix A regarding simulations in time domain (stress relaxations).

### 3.1.3 Convenience of virtual DMA

The direct results of the virtual DMA spectral solver can be very useful in such kind of applications to gain insight about physical theories of disordered structures. The multi-material proposed by Masurel et al. (2015) aimed at offering some alternative modeling of disordered microstructures in order to understand complex cooperative phenomena during a relaxation test (of any kind). The
DMA observables virtually obtained here can now be confronted to physical models established mainly from microscopic and statistical considerations in order to see if some consistency can be achieved from both approaches. Such types of physical models are plethora in the literature and for a pure pedagogical objective, we will here focus on two of them: the HN model (Havriliak Jr. and Havriliak 1995) and the JPC model established by Cavaille et al. (1989).

**HN model** Initially proposed on empirical basis as a fitting model for DMA outputs (real and imaginary components of the material relaxation function in frequency domain) obtained on complex systems, this model further receives physical significance from probabilistic approaches to relaxation phenomenon (Jurlewicz and Weron 2002).

The HN relaxation function has the analytical expression

$$G_{\text{HN}}^{*}(\omega) = \frac{1}{(1 + (j\omega/\omega_{M})^{\nu})^{\gamma}} \qquad (14)$$

where $\gamma$ and $\nu$ are non-integer coefficients.

**JPC model** This model comes from a molecular kinetic theory for the rheology of glass, which is considered valid in its extension to amorphous polymers and is twofold.

1. Physical mechanisms associated to uncorrelated rotational-transitional motions of molecules within a frozen bath of randomly distributed high energy sites. Assumed of a fast character, they generate an anelastic deformation which characterizes the so-called $\beta$ relaxation.

2. Physical mechanisms associated to cooperative, but hierarchically constrained, motions of cluster of molecules (referred to as shear microdomains in the paper for additional reasons). Assumed of slow character, they generate an irrecoverable deformation which characterizes the so-called $\alpha$ relaxation.

These mechanisms are fairly well recognized in the glass material scientific community, and have obviously been considered by Masurel et al. (2015) to propose such kind of heterodynamics modelling. It is then naturally considered here as capable to treat relaxation correctly regardless of the precise nature and scale of local interactions. As it offers a direct formulation of real and imaginary moduli as function of frequency, we can illustrate the interest of virtual DMA computations to confirm or validate assumptions introduced in physical approaches. The rheological equations read

$$G'_{\text{JPC}}(\omega) = \frac{J'(\omega)}{J'(\omega)^2 + J''(\omega)^2} + G_R \quad \text{and} \quad G''_{\text{JPC}}(\omega) = \frac{J''(\omega)}{J'(\omega)^2 + J''(\omega)^2} \qquad (15)$$

where the terms $J'(\omega)$ and $J''(\omega)$ are the sums of two contributions ($\beta$ and $\alpha$ relaxations)

$$J'(\omega) = J'_{\beta} + \frac{1}{G_g} \left[ 1 + (\omega\tau_m)^{-\kappa} \cos\left(\frac{\kappa\pi}{2}\right) + C(\omega\tau_m)^{-\xi} \cos\left(\frac{\xi\pi}{2}\right) \right] \qquad (16a)$$

$$J''(\omega) = J''_{\beta} + \frac{1}{G_g} \left[ (\omega\tau_m)^{-\kappa} \sin\left(\frac{\kappa\pi}{2}\right) + C(\omega\tau_m)^{-\xi} \sin\left(\frac{\xi\pi}{2}\right) \right] \qquad (16b)$$

where $J'_{\beta}$ and $J''_{\beta}$ are provided by Cavaille et al. (1989, Equation (27)) along with all additional necessary information. Formally, this model reaches an amount of ten parameters which would necessitate a very detailed sensitivity analysis study if any strong and metrological application to real data is required. Here we will perform an inverse parameter estimation without any prior sensitivity analysis considering that

1. some of them are known ($G_g$, $G_r$ which can be identified from low and high frequency asymptotes of the real and imaginary moduli but are known here as the inputs of the virtual DMA computations, $T_g$ the glass transition temperature of the polymer),

2. the others respect realistic physical numerical values,

3. no noise was added to virtual DMA data.

Figure 4 compares both HN and JPC models in their adjustment performance to the virtual DMA results obtained on the multidomain material and for the case $G_g/G_r = 10^3$. Adjustment of
<figure><img src="image_5.png" /><figcaption>**Figure 4** Virtual DMA outputs on the heterogeneous dynamics system:
(a) Cole-Cole plot
(b) storage modulus
(c) dissipative modulus</figcaption></figure>

the models to the virtual DMA data is obtained through a Least-Square minimization of the residuals norm, considering either a simplex or Levenberg-Marquardt algorithm in order to check the good behavior of the parameter estimation process (sensitivity to initial values, stability of the results). Input parameters for virtual DMA simulations of the elastic unrelaxed $G_g$ and relaxed $G_R$ moduli (row (1) in Table 2) are considered as known parameters in the Parameter Estimation Problem (PEP) for both the HN and JPC models. As a result, all supplementary identified parameters are involved in the sole modeling approach of the relaxation(s) kinetics (rows (2) and (3)). Note that in the JPC model, the parameter $\tau_{max}$ refers to the kinetics of the $\beta$ relaxation, the characteristic time associated to the $\alpha$ relaxation being explicitly absent from the parameters of the model but implicitly present and related to it through an assumption that parameter $\tau_{\alpha} = \tau_{max}\beta(A_{\alpha}G_g/\kappa)^{-1/\kappa}$—see (Cavaille et al. 1989, Equation (48)).

Figure 4(b) and Figure 4(c) plot the real and imaginary components, respectively, of the complex equivalent modulus at convergence of the PEP as a function of frequency. Figure 4(a) shows the corresponding Cole-Cole representation. It is clear that the HN model is unable to describe the outputs of virtual DMA simulations. The Cole-Cole plot especially brings out this deficiency. On the contrary, the JPC model captures the heterodynamics approach much more precisely. Parameters identified for the JPC model (row (3) in Table 2) can be seen to be very close or at least in the same order of magnitude as those considered for the simulation of the generic amorphous polymer in (Cavaille et al. 1989). This finally confirms that the ingredients used to simulate virtual DMA on a glassy system model are in phase with the ideas of physicists on the underlying mechanisms and demonstrate from our point of view all the benefits that can be drawn from spectral solvers used in virtual DMA mode.

## 3.2 Application 2: Virtual versus experimental DMA on polystyrene/glass beads composites

This second example concerns now a material which is a particulate composite made of glass beads (GB) considered to be elastic and embedded in a viscoelastic matrix of Polystyrene (PS). Experimental DMA results have been obtained with a temperature sweep at a given frequency. A temperature variation directly affects the relaxation time spectrum of the matrix and the overall composite behavior. Around the glass transition, microstructural configuration changes
<table><caption>**Table 2** Input parameters for the virtual DMA heterodynamics approach and the parameter estimated for the HN and JPC models. (4) Values considered in (Cavaille et al. 1989) for a generic amorphous polymer. Units: G<sub>g</sub> and G<sub>r</sub> in GPa, A<sub>α</sub><sup>-1</sup> and A<sub>β</sub><sup>-1</sup> in GPa<sup>-1</sup>, U and ΔU in kJ/mol; w: width of the normal distribution of times τ<sub>i</sub>.</caption><thead><tr><th colspan="5">Virtual DMA data</th></tr><tr><th>Input data</th><th>G<sub>g</sub></th><th>G<sub>r</sub></th><th>τ<sub>max</sub></th><th>w</th></tr></thead><tbody><tr><td>(1)</td><td>1</td><td>10<sup>-3</sup></td><td>1</td><td>2</td></tr><tr><th colspan="5">HN model</th></tr><tr><th>Known parameters</th><th>G<sub>g</sub></th><th>G<sub>r</sub></th><th></th><th></th></tr><tr><td></td><td>1</td><td>10<sup>-3</sup></td><td></td><td></td></tr><tr><th>Identified parameters</th><th></th><th></th><th>τ<sub>max</sub></th><th>ν</th></tr><tr><td>(2)</td><td></td><td></td><td>0.244</td><td>0.532</td></tr><tr><th colspan="5">JPC model</th></tr><tr><th>Known parameters</th><th>G<sub>g</sub></th><th>G<sub>r</sub></th><th>T<sub>g</sub></th><th></th></tr><tr><td></td><td>1</td><td>10<sup>-3</sup></td><td>355</td><td></td></tr><tr><th>Identified parameters</th><th>τ<sub>max</sub></th><th>κ</th><th>χ</th><th>A<sub>α</sub><sup>-1</sup></th></tr><tr><td>(3)</td><td>0.156</td><td>0.37</td><td>0.75</td><td>0.29</td></tr><tr><td>(4) (Cavaille et al. 1989)</td><td>8</td><td>0.3</td><td>0.95</td><td>0.3</td></tr><tr><td></td><td></td><td></td><td></td><td>0.1</td></tr><tr><td></td><td></td><td></td><td></td><td>U</td></tr><tr><td></td><td></td><td></td><td></td><td>ΔU</td></tr></tbody></table>

occur and probing kinetic response times makes DMA a very sensitive and adapted tool. All experimental data collected for this part were published in a series of three papers by Alberola and collaborators (Agbossou et al. 1993; Alberola and Bergeret 1994; Alberola and Mele 1996) and reported in terms of storage modulus and loss tangent variables. As a consequence, this example will first show how the virtual DMA code can work with tabulated pairs (real and imaginary parts) informing about the material behavior at each temperature. Secondly it will again illustrate the efficiency of a spectral code to produce either 2D or 3D DMA simulations for any type of composites from the knowledge of the individual constitutive materials. Finally, it will show that virtual DMA, as a computational tool applying on composite microstructures which mimic the topological reality in the best way, can either lead to refute previous theories or effective models or to cast doubt on the experimental data (which can advantageously lead to reconsider the metrology associated to its apparatus or the possible source of bias).

### 3.2.1 Data brought by experimental DMA

The DMA experiments on PS/GB particulate composites were conducted on a Metravib equipment. The aim was to check the modeling abilities of homogenization approaches of Generalized Self-Consistent (GSC) type. These later were based on available approaches in the elastic framework and extended to viscoelastic behavior through the correspondence principle.

In (Agbossou et al. 1993, Figures 1 and 7), DMA data are given for neat PS in terms of storage modulus, tan δ and complex Poisson coefficient respectively. These data have been digitized to provide (E*, ν*) or, after conversion, (K*, μ*) for implementation in CraFT solver within a tabulated format. Errors have certainly been introduced in this process. The main outputs of virtual DMA reported here will show that they have no impact on the results, those being rather fully determined by the extrapolation made by Agbossou et al. (1993) to produce experimental Poisson coefficient values.

These data are available in the 25-185 °C temperature range, 5-100 Hz frequency range, for composites with f<sub>v</sub> = 6 %, 15 %, 21 %, 35 %, 50 % volume fraction of glass beads, in two different size distributions: d<sub>1</sub> = 1-45 µm and d<sub>2</sub> = 70-110 µm. Unfortunately, no precise information is known regarding the real topology of the microstructures obtained from these distributions. At large volume fractions, aggregates or clusters of particles can have been formed, with unknown packing arrangement, in possible anisotropic configurations, etc. (Alberola and Mele 1996). The authors also mentioned an existing porosity of about 8 vol % in the specimen. This of course precludes perfectly realistic simulations, yet easily achieved with spectral solvers. Virtual DMA simulations will then be performed assuming a given topology and as a result some discrepancies can be expected. Regarding the constituent mechanical properties, ultrasonic measurements have
provided elastic constants of the polymer matrix with assumed high accuracy (E = 3.69 GPa, G = 1.38 GPa, ν = 0.33) and at room temperature the DMA results are in accordance with them. Note that DMA technique provides only indirectly the intrinsic moduli through overall rigidity measurements, which depend on the type of mechanical configuration (tension, flexion) and on length dimensions. In the absence of any details about that point, we observe that the authors report a good consistency between DMA data and US measurements. The glass transition temperature of PS is about 100 °C. Regarding GB, elastic constants are given as E = 73 GPa and ν = 0.2.

Virtual DMA isochronal scan simulations are performed directly from the knowledge of the pair (E', tan δ) versus temperature at 5 Hz for the PS matrix, as measured via DMA by Agbossou et al. (1993, Figures 1(a) and 1(b)), and the knowledge of a given distribution size ($d_1$ or $d_2$) and given volume fraction of glass beads $f_v$. Note that the log(E'') curve was not made available by the authors for pure PS but available for all types of composites. As a result, and because tan δ tends asymptotically to 0 at low temperatures (At 5 Hz and far from the transition region, a pure elastic behavior is expected.), the data collected from the figures are obviously affected by artifacts which will impact the calculated imaginary moduli of the composites.

### 3.2.2 Conclusions drawn from Alberola and collaborators

In the above-mentioned papers, experimental data deliver the following observations.

* As the volume fraction of particles increases (from 0 to 50% in their study),

(i) the magnitude of mechanical relaxation characterized by the maximum value of the tan δ peak decreases;

(ii) a shift of the tan δ peak towards higher temperatures is reported. This effect is not obvious from the data (Agbossou et al. 1993, Figure 3(c)) and could result from a small experimental bias (with regard to the temperature sweep and the thermalization of the specimen for example);

(iii) both glassy and relaxed moduli at low and high temperatures respectively show an increase.

* At constant volume fraction,

(iv) the increase in average size of the particles ($d_1 \rightarrow d_2$) results in an increase in the tan δ peak maximum (dissipation), or the reinforcement effect increases with decreasing the size of the particles. We have ($E_{vR}$)$_{d_1}$ > ($E_{vR}$)$_{d_2}$ and (tan δ)$_{d_1}$ < (tan δ)$_{d_2}$. Size effects depend obviously on the nature of the matrix/fillers couple and the review paper by Fu et al. (2008) reports for example that the modulus of a composite with high volume fraction of GB in an epoxy resin is nearly insensitive to size effects.

A GSC mean field theoretical approach (Kerner's and Christensen and Lo's models for the bulk and shear moduli respectively) was used in complex variables to explain these data. Only the observed fact (iii) was reproduced by the model, based on a *constant and real valued* Poisson ratio (whatever its value between the low temperature $v = 0.33$ and high temperature $v \rightarrow 0.5$). Especially, behaviors (i) and (ii) were not captured by the models: drop in maximum and shift in temperature of tan δ with an increase in volume fraction of GB (Agbossou et al. 1993, Figures 5 and 6). This led the authors to invoke a complex-valued Poisson coefficient and a strategy to determine it on the considered temperature range. This point will deserve a large discussion in the next section.

### 3.2.3 Virtual DMA results and implications regarding a comprehensive analysis of the role played by the Poisson ratio of the matrix

DMA simulations of the particulate composite are based on a synthetic material created with the Random Sequential Algorithm. For large volume fractions (up to 50%), this latter is combined with a metropolis algorithm to reduce computation times. Note that without additional reported information, the two size distributions of the particles are considered as Gaussian. The chosen $R_{max}/R_{min}$ ratio and volume fraction $V_f$ allow distributing a number of spheres N of mean radius $r_{mean}$ in a square box of unit dimensions.

**Influence of the volume fraction of the particles** Because this point was the major fact discussed by Agbossou et al. (1993), we first present some results obtained with the virtual
DMA solver along with the available experimental data. Virtual DMA simulations have been
performed for three temperature-independent *matrix* Poisson coefficients: $v_1$ = 0.33, $v_2$ = 0.49
and $v_3$ = 0.4999. Note that the values 0.49 and 0.4999 are selected to represent the incompressible
behavior expected in rubbery state (at high temperature) but without any possible precise
experimental confirmation. Figure 5(a) and Figure 5(b) compare experimental and simulated data
for the storage modulus and the loss angle tangent in the case of a size distribution $d_1$. The

<figure><img src="image_6.png" /><figcaption>**Figure 5** Comparison of virtual DAM simulations and DMA experiments for two volume fractions. (a) Storage modulus and (b) loss tangent versus temperature. Distribution $d_1$.</figcaption></figure>

following observations can be made (same as for distribution $d_2$).

1. For a volume fraction of 6% of particles, virtual DMA solutions are very close to the experimental data for the pure PS matrix. The latter data are used as input for the matrix properties so that this result is expected. Any mean-field approach would lead to the same results, as a matter of fact for small concentrations. Note that experimental data (full circles) nevertheless depart substantially from the matrix data. This naturally results from a mismatch between the intrinsic ideal character of theoretical models (whatever they are) and the experimental conditions actually followed. On one hand, the deficiencies in the model could be relative to a lack of knowledge on the specimen (porosity, aggregates, interphases, statistical distribution of sizes unknown...). On the other hand, one could incriminate biases in the experimentation. Regarding DMA one has to keep in mind that intrinsic measurements, like the real and imaginary moduli provided here, are obtained through an overall stiffness measurement. They do already depend from an identification model and its assumptions (stability in size, loading conditions through grips...). Because this discrepancy arises only in the transition range between 110 °C-130 °C, a bias due to the temperature regime could also be suspected, especially in the temperature sweep mode.

2. For a volume fraction of 50% of particles, virtual DMA solutions produce a tan $\delta$ curve very close to the experimental ones with a strong peak drop at the transition temperature—compare with (Agbossou et al. 1993, Figures 5 and 6). The computations show also a very high sensitivity to the value of the Poisson ratio: considering $v_2$ or $v_3$ values for the incompressible state highly affects the solution. The same sensitivity was shown previously when computing bounds of the homogenized bulk modulus, see (Gibiansky and Lakes 1993, Figure 4). These two results clearly escape to the mean-field approach ('GSC model' curve reported in Figure 5(b), quasi-insensitive to the Poisson ratio). In virtual DMA solutions with constant real Poisson ratio, experimental data are well described by the value $v_2$ = 0.49 for both the storage modulus and loss angle. It is worthwhile to note that this sensitivity logically begins only as the material enters into the vitreous transition *i.e.* above temperatures of 110 °C, when the mechanical properties of the PS matrix drastically decreases towards the rubber state. Thus, the value of 0.49 is very sounded. The temperature dependence of the Poisson ratio is taken into account in the next subsection and confirms that point.

3. The increase in volume fraction of glass beads makes both glassy and rubber moduli greater at low and high temperature respectively. This fact was the only one observed similarly with mean-field GSC approaches ('GSC model' curve in Figure 5(a)). However, DMA simulations produce a
slightly higher reinforcement at low temperature compared to the experiment. The porosity of
8 % considered by the authors as a result of the composite elaboration process can explain this
fact. Note that according to theoretical simulations by Remillat (2007, Figure 2), our ratio of
elastic moduli particle/matrix is 20 and we would expect at this volume fraction of particles a
2.5 factor gain in stiffness. Virtual DMA gives an homogenized elasticity of about 9.8 GPa at
60 °C for a 50 % volume fraction which falls exactly as 2.5 × 3.7 GPa. This again could confirm a
pretty good consistency of the virtual DMA predictions and consequently that the discrepancies
with experimental data may primarily be due to a specimen microstructure mismatch between
idealized and real composites. Regarding the behavior of the storage modulus at high temperature,
it can be clearly seen in Figure 5(a) that the GSC model underestimates its value when compared
to virtual DMA.

Virtual DMA thus directly confirms observations (i) and (iii), without any requirement of a
Poisson coefficient being a complex quantity, the central argument of the paper by Agbossou
et al. (1993), to correct the model-experiment inconsistency at high volume fraction of particles.
However, their two-phase model (Christensen and Lo 1979) clearly overestimates the maximum
of tan δ peak—Agbossou et al. 1993, Figure 8(c) and GSC model in Figure 5(b). Later in (Alberola
and Mele 1996), observation (i) was reproduced theoretically using a four-phase model (Hervé
and Zaoui 1993), invoking possible aggregation of particles at high volume fractions (50 %) but
considering a Poisson coefficient as a real quantity. Again, it is clear that virtual DMA strongly
questions this assumption as it is not required here in our simulations. The better agreement
obtained by the authors between experimental data and their model probably just results from
an increase in the degrees-of-freedom of the model, going from the two-phase to four-phase
self-consistent approach. Note also that the complex-valued Poisson coefficient extrapolated
by Agbossou et al. (1993) for the whole considered temperature range, when converted into
bulk and shear moduli (for introduction in the constitutive law of the Craft solver), produces a
negative real bulk modulus which is not physical.

Observation (ii) is not confirmed by the virtual simulations. The peak location of the
composite is totally determined by the tan δ peak of the matrix (input data). It is an obvious
manifestation of the temperature determined physical transition of the matrix only and to
the known fact that a high damping of the compliant phase gives rise to high damping in the
composite although the moduli of this phase are much smaller in absolute value than the moduli
of the stiff elastic phase (Gibiansky and Lakes 1997). The loss factor curve of the mixture is ruled
by that of the matrix, regardless of the volume fraction, if the inclusion is either very stiff or very
soft. In that view, it extends the result established by Hashin (1970), according to which the
effective and matrix loss tangents are the same, but only very far from the transition temperature,
that is far from the tan δ peak, and for very low tan δ (weak viscoelastic effect). Our conclusion is
that the experimental-theoretical mismatch stems from the sole thermal aspects of the problem:
DMA experiments in temperature sweep mode assumes a quasi-stationary state; mechanical
simulations assume perfect thermalization of both phases.

**Influence of the size distribution of the particles** Simulations shown in Figure 6 were
performed for the two size distributions of the particles *d*₁ and *d*₂. Regarding experimental
observation (iv), virtual DMA simulations contradict this point. In the case of constant Poisson ratio
for the matrix equal to 0.49, we observe from the simulations that an increase in average size of the
particles at the same constant volume fraction (f₀ = 50 %) increases the reinforcement: (E*ν*<sub>T</sub>)*d*₁*<sub>1</sub>* <
(E*ν*<sub>T</sub>)*d*₂* and (tan δ)*d*₁* > (tan δ)*d*₂*. This is the exact contrary of what is reported experimentally.
Note that this question could not have been addressed by Alberola and collaborators within
the Self-Consistent models they used (taking into account a volume fraction parameter but
no parameter relative to any particle size). For the two specimen at 50 % volume fraction and
distributions *d*₁ or *d*₂, the conclusion is then that the model can not explain the data when an
equal matrix Poisson ratio behavior is considered.

With some reasons given previously, an experimental bias could be invoked to explain this
conflict. But because it is clear from Figure 6 that the sensitivity to the Poisson ratio is of crucial
concern, the Poisson ratio of the matrix is now considered to vary with temperature (noted vT in
that case) according to a S-shape function, starting at vT = v1 at room temperature, and with an
<figure><img src="image_7.png" /><figcaption>**Figure 6** Storage modulus versus temperature. Distributions $d_1$ and $d_2$ for volume fraction $f_v = 50\%$.</figcaption></figure>

asymptotic value $v_T$ at high temperature (180 °C) that is supposed to tend towards 1/2 because the matrix tends to an incompressible behavior. This fact was demonstrated experimentally long times ago in the works of Kono (1960) and Waterman (1963) with $v_T$ of PolyStyrene considered to vary from 0.35 → 0.5 or 0.3 → 0.45 respectively for the two authors. If we apply the procedure reported in (Agbossou et al. 1993) to reconstruct $v_T$ from experimental DMA data, with the constraint $v_T → v_2 = 0.49$ and $v_T$ always considered as a real-valued quantity, it can be seen in the plots of Figure 6 that the extreme sensitivity to the Poisson ratio could explain such conflict: for a same volume fraction of GB, a Poisson ratio considered as varying with temperature and $v_T → v_2 = 0.49$ for distribution $d_1$ leads to an overall curve of rigidity now surpassing the curve obtained for distribution $d_2$ with a Poisson ratio considered as constant and equal to $v_2 = 0.49$. DMA experimental data have been obtained up to 190 °C, a temperature still very far from the fusion temperature where a drastic drop in modulus occurs. Hence the rubbery state can be more or less established. Additionally, DMA experiments conducted under temperature sweep condition impose a given temperature rate. Specimen with different size distribution for the glass beads can easily produce slightly different thermal conductivity, a parameter which can change substantially the thermal regime in sweeping mode and hence change the behavior with respect to the incompressibility limit. The presence of voids in the highly concentrate case (already discussed) can also affect this parameter. Anyway, if the reasons explaining this possible conflict cannot be firmly established, on the other hand, this issue shows how virtual DMA simulation can contribute to open avenues for reflection. In the present case, it suggests strong investigations relative to the precise measurement of Poisson ratio of the matrix with temperature, especially when it tends to a rubbery state.

Finally, as shown in the first example, outputs of virtual DMA in terms of mechanical fields can be valuable to get deeper analysis of phase transition. We select in Figure 7 the volumetric strain fields of the matrix, calculated for the same parameters as for Figure 6. Two maps are shown, obtained for a temperature below (60 °C) and above (160 °C) the glass transition temperature $T_g$ and for an input average stress of 1 Mpa. These fields are normalized with respect to the longitudinal modulus as it changes drastically when crossing the transition zone. When compared to the map at 60 °C, it can be seen that the volumetric strain obtained above $T_g$ tends to a more homogeneous zero-value everywhere in the matrix, except in the very limited area confined between two beads. Note that in this simulation, the Poisson ratio tends only to 0.49. Figure 7(a) shows the temperature-dependent overall volumetric strain averaged over a selection of ten different configurations along with the volumetric strains in both the GB and PS matrix phases. The transition regime marks the decrease of both volumetric strains. In the GB, it tends to zero within numerical precision when raising the temperature. In the matrix, it decreases strongly of nearly one order of magnitude at the transition. It slightly increases beyond probably because the matrix sustains the full imposed mechanical excitation and, as shown in the strain maps, strong volume expansion may exist locally because Poisson's ratio is not exactly equal to 0.5.
<figure><img src="image_8.png" /><figcaption>**Figure 7** Volumetric strain. Same input data as Figure 6:
(a) Evolution with temperature averaged over ten configurations
(b) Map at T = 60 °C
(c) Map at T = 160 °C</figcaption></figure>

## 4 Conclusion

The study presented here was intentionally based upon data provided in previously published papers, either theoretical or purely experimental, to show the potentialities of virtual DMA. As a simulation tool based on an FFT solver, it turns out to be a very competitive alternative path to classical FEM in terms of computation times. For experimentalists, it offers a direct access to the frequency-dependent observables of this widely used technique: storage and dissipative moduli, loss tangent. For a given heterogeneous RVE, it provides at the same time the full-field strain or stress maps, in either 2D or 3D, which can be valued to push forward our understanding of its mechanical behavior. Due to a favorable precision/computation times ratio, it is expected that in a near future, such codes will be able to perform inverse parameter identification on the collected experimental data. Hence it will help to question the validity of homogenized models or the experimental technique itself with regards to the way data are produced. But one of the most interesting perspective that virtual DMA offers is maybe the possibility to check the consistency of an idealized microstructure with respect to viscoelasticity. Based on MEB, X-ray tomography, AFM or any other imaging technique of the microstructural organization of some material, the complementary observation of the macroscopic behavior of this microstructure in a DMA experiment, both real and virtual, will allow fixing this microstructure idealization as consistent either in terms of morphological parameters (volume fraction of constituents, aspect ratio of the phase...) and mechanical parameters (moduli, Poisson ratio, relaxation time spectra...). Considering such a way to firmly assess the microstructure model, routes will then be open to consider the variety of expected behaviors at a larger deformation range (plasticity, hardening, damage...).

**Acknowledgements** The authors wish to thank Hervé Moulinec (Laboratoire de Mécanique et Acoustique, Marseille, UMR CNRS 7031) for a fruitful introduction to CraFT and guidance for use.

## A Simulation of local fields in the transient state case of stress relaxations. Comparison between CraFT spectral solver and FEM

The test case of Section 3.1 is considered (Masurel et al. 2015). The local behavior of the spectral solver approach can also be compared to the FEM ones by comparing maps of the microdomain's
stress relaxations (simulations performed now in time domain for both numerical approaches). The exact same behavior is obtained in terms of the various stages leading to the establishment of a network of paths supporting the global stress, as illustrated in Figure A.1.

<figure><img src="image_9.png" /><figcaption>**Figure A.1** Calculated reduced shear modulus $\sigma_{12}^i(t)/\langle\sigma_{12}^i(t)\rangle$ of the macroscopic relaxation at different times: CraFT [left] and (Masurel et al. 2015, Fig. 11) [right] (Reprinted (adapted) with permission. Copyright (2015) American Chemical Society). Calculations performed with $G_g/G_r = 10^3$, $\mu = \log(\tau_{max}) = 0$, $w = 2$. $t^* = (a) 3 \times 10^{-4}$, (b) $10^{-2}$, (c) $10^{-1}$, (d) 1, (e) 10, (f) 40, (g) 400, (h) $10^4$, (i) $10^6$.</figcaption></figure>

CraFT (and FFT solvers more generally) appears then as a numerical way of obtaining both harmonic (complex) and temporal simulations in the same tool and in a single simulation run. DMA apparatus offer generally the additional possibility to perform creep or relaxation tests. Based on such a set of composite experimental data, it is then possible to chain both computations in complex and time variables in order to perform inverse identification of material parameters. For such purpose, reduced computational times is of utmost importance. The well-known performance of FFT solvers in this domain is increased in this perspective, and would be undoubtedly a nice option for material scientists.

## References

Agbossou, A., A. Bergeret, K. Benzarti, and N. Alberola (1993). Modelling of the viscoelastic behaviour of amorphous thermoplastic/glass beads composites based on the evaluation of the complex Poisson's ratio of the polymer matrix. *Journal of Materials Science* 28(7):1963–1972. [DOI].

Alberola, N. and A. Bergeret (1994). Physical modeling of the interphase in amorphous thermo-plastic/glass bead composites. *Polymer Composites* 15(6):442–452. [DOI].

Alberola, N. D. and P. Mele (1996). Viscoelasticity of polymers filled by rigid or soft particles: Theory and experiment. *Polymer Composites* 17(5):751–759. [DOI].

Armines & ONERA (2021). *Zébulon*. Version 9.1.1. [URL].

Brinson, L. C. and W. G. Knauss (1992). Finite Element Analysis of Multiphase Viscoelastic Solids. *Journal of Applied Mechanics* 59(4):730–737. [DOI].

Cahill, D. G. (1990). Thermal conductivity measurement from 30 to 750 K: the $3\omega$ method. *Review of Scientific Instruments* 61(2):802–808. [DOI].

Cavaille, J. Y., J. Perez, and G. P. Johari (1989). Molecular theory for the rheology of glasses and polymers. *Physical Review B* 39(4):2411–2422. [DOI].

Charière, R., A. Marano, and L. Gélébart (2020). Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites. *International Journal of Solids and Structures* 182-183:1–14. [DOI], [OAI].

Christensen, R. and K. Lo (1979). Solutions for effective shear properties in three phase sphere and cylinder models. *Journal of the Mechanics and Physics of Solids* 27(4):315–330. [DOI].

Emri, I. and N. Tschoegl (1995). Determination of mechanical spectra from experimental responses. *Time Dependent Problems in Mechanics* 32(6):817–826. [DOI].
Eyre, D. J. and G. W. Milton (1999). A fast numerical scheme for computing the response of composites using grid refinement. *The European Physical Journal Applied Physics* 6(1):41-47. [DOI].

Figliuzzi, B., D. Jeulin, M. Faessel, F. Willot, M. Koishi, and N. Kowatari (2016). Modelling the Microstructure and the Viscoelastic Behaviour of Carbon Black Filled Rubber Materials from 3D Simulations. *Technische Mechanik - European Journal of Engineering Mechanics* 32(1-2):22-46. [DOI].

Fu, S.-Y., X.-Q. Feng, B. Lauke, and Y.-W. Mai (2008). Effects of particle size, particle/matrix interface adhesion and particle loading on mechanical properties of particulate-polymer composites. *Composites Part B: Engineering* 39(6):933-961. [DOI].

Gallican, V. and R. Brenner (2019). Homogenization estimates for the effective response of fractional viscoelastic particulate composites. *Continuum Mechanics and Thermodynamics* 31(3):823-840. [DOI], [OAI].

Gallican, V., R. Brenner, and P. Suquet (2017). Exact asymptotic relations for the effective response of linear viscoelastic heterogeneous media. *Comptes Rendus Mécanique* 345(11):742-751. [DOI].

Ghazvinian, E., M. Diederichs, and R. Quey (2014). 3D random Voronoi grain-based models for simulation of brittle rock damage and fabric-guided micro-fracturing. *Journal of Rock Mechanics and Geotechnical Engineering* 6(6):506-521. [DOI].

Gibiansky, L. and R. Lakes (1993). Bounds on the complex bulk modulus of a two-phase viscoelastic composite with arbitrary volume fractions of the components. *Mechanics of Materials* 16(3):317-331. [DOI].

Gibiansky, L. and R. Lakes (1997). Bounds on the complex bulk and shear moduli of a two-dimensional two-phase viscoelastic composite. *Mechanics of Materials* 25(2):79-95. [DOI].

Gélébart, L. and F. Ouaki (2015). Filtering material properties to improve FFT-based methods for numerical homogenization. *Journal of Computational Physics* 294:90-95. [DOI], [ARXIV].

Hashin, Z. (1970). Complex moduli of viscoelastic composites—I. General theory and application to particulate composites. *International Journal of Solids and Structures* 6(5):539-552. [DOI].

Havriliak Jr., S. and S. J. Havriliak (1995). Comparison of dielectric theories that explicitly include viscoelastic parameters. *Journal of Polymer Science Part B: Polymer Physics* 33(16):2245-2252. [DOI].

Hervé, E. and A. Zaoui (1993). *n*-Layered inclusion-based micromechanical modelling. *International Journal of Engineering Science* 31(1):1-10. [DOI].

Jonscher, A. K. (1996). *Universal relaxation law*. English. Chelsea Dielectrics Press.

Jurlewicz, A. and K. Weron (2002). Relaxation of dynamically correlated clusters. *Journal of Non-Crystalline Solids* 305(1):112-121. [DOI], [CITESEERX].

Kim, S.-Y. and D.-H. Lee (2009). Identification of fractional-derivative-model parameters of viscoelastic materials from measured FRFs. *Journal of Sound and Vibration* 324(3):570-586. [DOI].

Kono, R. (1960). The Dynamic Bulk Viscosity of Polystyrene and Polymethyl Methacrylate. *Journal of the Physical Society of Japan* 15(4):718-725. [DOI].

Lebensohn, R., O. Castelnau, R. Brenner, and P. Gilormini (2005). Study of the antiplane deformation of linear 2-D polycrystals with different microstructures. *International Journal of Solids and Structures* 42(20):5441-5459. [DOI].

Mareau, C. and C. Robert (2017). Different composite voxel methods for the numerical homogenization of heterogeneous inelastic materials with FFT-based techniques. *Mechanics of Materials* 105:157-165. [DOI], [OAI].

Masurel, R. J., S. Cantournet, A. Dequidt, D. R. Long, H. Montes, and F. Lequeux (2015). Role of Dynamical Heterogeneities on the Viscoelastic Spectrum of Polymers: A Stochastic Continuum Mechanics Model. *Macromolecules* 48(18):6690-6702. [DOI].

Menard, K. and N. Menard (2015). Dynamic Mechanical Analysis in the Analysis of Polymers and Rubbers. *Encyclopedia of Polymer Science and Technology*. John Wiley & Sons. [DOI].

Monchiet, V. (2018). FFT based iterative schemes for composites conductors with non-overlapping fibers and Kapitza interface resistance. *International Journal of Solids and Structures* 135:14-25. [DOI].
Moulinec, H. and F. Silva (2014). Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. *International Journal for Numerical Methods in Engineering* 97(13):960-985. [DOI], [OAI].

Moulinec, H. and P. Suquet (1998). A numerical method for computing the overall response of nonlinear composites with complex microstructure. *Computer Methods in Applied Mechanics and Engineering* 157(1):69-94. [DOI], [ARXIV].

Moulinec, H., P. Suquet, A. Labé, F. Silva, D. Garajeu, and G. Boittin (2004). *CraFT: Composite response and Fourier Transforms*. Version 1.2. [URL].

Nielsen, L. E. and R. F. Landel (1994). *Mechanical properties of polymers and composites*. 2nd ed. M. Dekker.

Quey, R. (2021). *Neper: Polycrystal Generation and Meshing*. Version 4.1.2. [URL].

Remillat, C. (2007). Damping mechanism of polymers filled with elastic particles. *Mechanics of Materials* 39(6):525-537. [DOI].

Roters, F., M. Diehl, P. Shanthraj, P. Eisenlohr, C. Reuber, S. Wong, T. Maiti, A. Ebrahimi, T. Hochrainer, H.-O. Fabritius, S. Nikolov, M. Friák, N. Fujita, N. Grilli, K. Janssens, N. Jia, P. Kok, D. Ma, F. Meier, E. Werner, M. Stricker, D. Weygand, and D. Raabe (2019). DAMASK – The Düsseldorf Advanced Material Simulation Kit for modeling multi-physics crystal plasticity, thermal, and damage phenomena from the single crystal up to the component scale. *Computational Materials Science* 158:420-478. [DOI].

Salnikov, V., D. Choï, and P. Karamian-Surville (2015). On efficient and reliable stochastic generation of RVEs for analysis of composites within the framework of homogenization. *Computational Mechanics* 55(1):127-144. [DOI], [ARXIV].

Molecular Organization and Dynamics (2008). *Characterization and Analysis of Polymers*. Ed. by A. Seidel. Wiley, pp 649-675.

Suquet, P. (2012). Four exact relations for the effective relaxation function of linear viscoelastic composites. *Recent Advances in Micromechanics of Materials* 340(4):387-399. [DOI], [OAI].

Uchic, M. D., L. Holzer, B. J. Inkson, E. L. Principe, and P. Munroe (2007). Three-Dimensional Microstructural Characterization Using Focused Ion Beam Tomography. *MRS Bulletin* 32(5):408-416. [DOI].

Waterman, H. A. (1963). Determination of the complex moduli of viscoelastic materials with the ultrasonic pulse method (Part I). *Kolloid-Zeitschrift und Zeitschrift für Polymere* 192(1):1-8. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license, you may need to obtain permission from the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** Stéphane André: conceptualization, numerical investigation, writing original draft. Julien Boisse: code development, numerical investigation, review. Camille Noûs: collective co-author embodying the grounding role of the collective in the creation of scientific knowledge, proceeding from a collegial organization, from the legacy and criticism of previous work, and from peer-review.

**Supplementary Material** The source code used in the present article can be found at the permalink: swh:1:dir:0135 7308677370b558b1a9ecb6dodca9bb8ea18;origin=https://gitlab.com/lemta-rheosol/craft-virtual-dma;visit=swh:1 snp:ef1a939275f05b667e189afbeed5fd59cca51c9d;anchor=swh:1:rev:ad74f6a7f73c7906f9b36ee28dd006231f42552e.

**Acknowledgements** The authors gratefully acknowledge financial support from the CARNOT label for its specific financial support #ICEEL-2018, the EMPP (Energy, Mechanics, Process, Products) Research Axis of Lorraine University, the Fédération Jacques Villermaux (FR 2863 CNRS). Fundings were brought also through the research program for young scientists of the "Grand-Est" Region of France.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.