T T Le

J Guilleminot
email: johann.guilleminot@u-pem.fr

C Soize

Stochastic continuum modeling of random interphases from atomistic simulations. Application to a polymer nanocomposite

Keywords: Interphase, maximum entropy, Molecular Dynamics, nanocomposite, probabilistic model, random field, inverse problem

Stochastic continuum modeling of random interphases from atomistic simulations. Application to a polymer nanocomposite

1

Introduction

The prediction of the enhanced physical properties exhibited by heterogeneous materials reinforced with nanoscopic fillers is a great challenge of nanoscale science. For such systems, small-scale interactions (between the fillers and the matrix phase, or between adjacent fillers) are no longer negligible and must be accounted for in a multiscale setting. One difficulty pertaining to modeling issues is that the final properties strongly depend on the coupling between various parameters, such as size effects or filler aggregation [START_REF] Fu | Effects of particle size, particle/matrix interface adhesion and particle loading on mechanical properties of particulate-polymer composites[END_REF][START_REF] Móczó | Polymer micro and nanocomposites: Structure, interactions, properties[END_REF]. Perhaps one the most discussed specificity of nanoreinforced composites is the existence of an interphase region surrounding the inclusions. Such a region is commonly identified with the region confined to inclusion neighborhood where the polymer chains exhibit conformational and geometrical properties that are substantially different from those in the bulk polymer phase. The experimental characterization of such properties, by using high-resolution nuclear magnetic resonance (NMR) techniques for instance, gave rise to an extensive literature: see e.g. [START_REF] Harton | Immobilized polymer layers on spherical nanoparticles[END_REF][START_REF] Holt | Dynamics at the Polymer/Nanoparticle Interface in Poly(2-vinylpyridine)/Silica Nanocomposites[END_REF][START_REF] Papon | Low-Field NMR Investigations of Nanocomposites: Polymer Dynamics and Network Effects[END_REF][START_REF] Fragiadakis | Dynamics near the filler surface in natural rubber-silica nanocomposites[END_REF][START_REF] Petrović | Effect of silica nanoparticles on morphology of segmented polyurethanes[END_REF][START_REF] Berriot | Filler-elastomer interaction in model filled rubbers, a 1 H NMR study[END_REF][START_REF] Berriot | Reinforcement of model filled elastomers: characterization of the crosslinking density at the filler-elastomer interface by 1 H NMR measurements[END_REF] for systems reinforced with Silica nanoinclusions, and [START_REF] Leu | NMR Characterization of Elastomer-Carbon Black Interactions[END_REF] for similar results with carbon black fillers. All these results show a local disturbance in the mobility of the polymer chains, as well as density fluctuations vanishing away from the inclusions. These observations are well correlated by atomistic simulations, which have additionaly highlighted a preferred orientation of the polymer chain segments tangentially to the particle surface; see [START_REF] Barbier | Interface between End-Functionalized PEO Oligomers and a Silica Nanoparticle Studied by Molecular Dynamics Simulations[END_REF][START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF][START_REF] Brown | A molecular dynamics study of a model nanoparticle embedded in a polymer matrix[END_REF][START_REF] Starr | Molecular Dynamics Simulation of a Polymer Melt with a Nanoscopic Particle[END_REF][START_REF] Odegard | Modeling of the mechanical properties of nanoparticle/polymer composites[END_REF][START_REF] Ndoro | Interface of grafted and ungrafted silica nanoparticles with a polystyrene matrix: Atomistic molecular dynamics simulations[END_REF][START_REF] Ndoro | Interface and interphase dynamics of polystyrene chains near grafted and ungrafted silica nanoparticles[END_REF][START_REF] Ghanbari | Interphase structure in silica-polystyrene nanocomposites: A coarse-grained molecular dynamics study[END_REF][START_REF] Chen | Structure and Dynamics Properties at Interphase Region in the Composite of Polystyrene and Cross-Linked Polystyrene Soft Nanoparticle[END_REF] for Molecular Dynamics (MD) simulations, as well as [START_REF] Vacatello | Monte carlo simulations of polymer melts filled with solid nanoparticles[END_REF][START_REF] Vacatello | Monte carlo simulations of the interface between polymer melts and solids. effects of chain stiffness[END_REF][START_REF] Vacatello | Molecular arrangements in polymer-based nanocomposites[END_REF][START_REF] Vacatello | Predicting the molecular arrangements in polymer-based nanocomposites[END_REF][START_REF] Vogiatzis | Monte carlo simulations of a coarse grained model for an athermal all-polystyrene nanocomposite system[END_REF] for the use of a Monte Carlo solver.

From a mechanical point of view, two types of approach were developed so as to estimate the effective properties of nanocomposites while accounting for size effects. A first class of methods relies on the integration of so-called interface models, such as spring-layer or membrane-type imperfect interface ones, in various extensions of classical (e.g. Eshelby-type) micromechanical frameworks. In this case, surface and interfacial energy concepts are invoked and plugged into asymptotic formulations, hence allowing for an explicit dependance of effective moduli on the nanoscale characteristic length (the radius of a spherical filler, for instance); see [START_REF] Duan | A unified scheme for prediction of effective moduli of multiphase composites with interface effects. part i: Theoretical framework[END_REF][START_REF] Duan | A unified scheme for prediction of effective moduli of multiphase composites with interface effects. part ii: Application and scaling laws[END_REF][START_REF] Le Quang | Variational principles and bounds for elastic inhomogeneous materials with coherent imperfect interfaces[END_REF][START_REF] Le Quang | Size-dependent effective thermoelastic properties of nanocomposites with spherically anisotropic phases[END_REF] and the references therein to name a few. A major issue related to such approaches is the definition of the interface stiffness, which may not be necessarily positive-definite (see [START_REF] Shenoy | Atomistic calculations of elastic properties of metallic fcc crystal surfaces[END_REF] for a derivation based on MD simulations, for instance). In the second class of models, the interphase is considered as an additional phase with a finite volume and standard (in the thermodynamical sense) mechanical properties. Such formulations can be subsequently used within (semi-)analytical homogenization schemes, such as the generalized self consistent scheme [START_REF] Odegard | Modeling of the mechanical properties of nanoparticle/polymer composites[END_REF][START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF][START_REF] Marcadon | Confrontation between molecular dynamics and micromechanical approaches to investigate particle size effects on the mechanical behaviour of polymer nanocomposites[END_REF][START_REF] Paliwal | Estimation of anisotropic elastic properties of nanocomposites using atomistic-continuum interphase model[END_REF], or coupled with computational approaches [START_REF] Peng | Modeling of nano-reinforced polymer composites: Microstructure effect on Young's modulus[END_REF]. It should be noted that these approaches are sequential, in the sense that they rely on an a priori characterization of interphase properties. Energy-based coupling techniques were alternatively developed in order to couple, in a concurrent manner, an atomistic description of the interphase with a finite element formulation (see [START_REF] Pfaller | An Arlequinbased method to couple molecular dynamics and finite element simulations of amorphous polymers and nanocomposites[END_REF] and the references therein, for instance).

Despite the intrinsic stochastic nature of the physical phenomena occurring while processing nanocomposites, most of the multiscale models developed so far are deterministic and presume that the interphase properties are known. More recently, a few studies involving some uncertainty quantification aspects were proposed. In a series of papers, Vu Bac and coworkers have investigated sensitivity analysis of multiscale predictions by using forward propagation of parametric uncertainties related either to input parameters of MD simulations or to other parameters (such as the volume fraction of fillers) in a continuum mechanics formulation [START_REF] Vu-Bac | Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs)[END_REF][START_REF] Vu-Bac | Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations[END_REF][START_REF] Vu-Bac | Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters[END_REF][START_REF] Vu-Bac | A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites[END_REF] -see also [START_REF] Rizzi | Uncertainty quantification in md simulations. part i: Forward propagation[END_REF][START_REF] Rizzi | Uncertainty quantification in md simulations. part ii: Bayesian inference of force-field parameters[END_REF] for the modeling of MD finite sampling noise and uncertainties in force field parameters, as well as [START_REF] Salloum | A stochastic multiscale coupling scheme to account for sampling noise in atomistic-to-continuum simulations[END_REF] for a study where uncertainties due to MD finite sampling are exchanged back-and-forth between an atomistic model and a continuum one. In [START_REF] Chang | Multiscale homogenization model for thermoelastic behavior of epoxy-based composites with polydisperse SiC nanoparticles[END_REF], an inverse identification procedure for some parameters modeled as random variables (such as particle radius or distribution, or matrix mechanical properties) is proposed. Finally, uncertainty propagation from nanoscale to structural analysis has been performed in [START_REF] Ghasemi | Uncertainties propagation in metamodel-based probabilistic optimization of CNT/polymer composite structure using stochastic multi-scale modeling[END_REF]. This work is devoted to the complementary and fundamental aspect of model construction and inverse calibration based on atomistic simulations. Here, and contrary to previous works developed elsewhere, the interphase elastic properties are modeled by a non-Gaussian tensor-valued random field. MD simulations are used so as to infer some basic properties of the latter, such the material symmetry locally exhibited or the underlying correlation structure, and to construct a virtual database devoted to the model calibration. This identification step is subsequently addressed by solving a statistical inverse problem stating the equivalence of the apparent properties obtained from atomistic computations and those estimated from stochastic homogenization in a continuum mechanics formulation.

This paper is organized as follows. Section 2 is devoted to the definition and virtual mechanical testing of the nanosystem under consideration. In particular, the atomistic models are introduced, together with the assembly procedure for the nanocomposite. The stochastic model is then introduced in Section 3. The methodology for inverse calibration is finally presented in Section 4.

Modeling and virtual testing of model nanocomposites

This section is devoted to the modeling and mechanical testing of a prototypical nanocomposite made up of a polymer matrix filled with a silica nanoinclusion. The overall methodology associated with the MD simulations and related issues is mostly borrowed from [START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF]. The atomistic descriptions of the polymer and silica particle are first introduced in Sec. 2.1 and Sec. 2.2, respectively. The assembly procedure of the nanostructure is then summarized in Sec. 2.3. The framework for virtual testing and results are finally presented in Sec. 2.4. All MD simulations are carried out in the LAMMPS package [START_REF] Plimpton | Fast parallel algorithms for short-range molecular-dynamics[END_REF], with a time step set to 2 fs. Atomistic images are produced by using the OVITO program [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito-the open visualization tool[END_REF].

Let r i denote the position of atom i, 1 i N , and let r = (r 1 , . . . , r N ). Let r ij = r jr i and r ij = r ij . The notation i -j is used hereafter to denote the bond between atoms i and j.

Molecular model of polymer matrix

In this study, a prototypical generic model of an amorphous linear polymer is considered. More specifically, the system under consideration is made up of n pc polymer chains, each of which containing 1, 000 united-atom CH 2 sites.

Force field description

The functional form of the force field in the polymer system classically reads as

V(r) = V B (r) + V UB (r) , (1)
where

V B (r) = {i,j} ∈ B b V b (r ij ) + {i,j,k} ∈ B θ V θ (θ ijk ) + {i,j,k, } ∈ B φ V φ (φ ijk ) (2)
is a term modeling interactions between m-uplets of bonded atoms (2 m 4), and

V UB (r) = N i=1 j>i V LJ (r ij ) (3)
models non-bonded interactions. In Eq. ( 2), θ ijk denotes the angle (r ji , r jk ) between bonds i-j and j -k, and φ ijk is the angle between the planes spanned by vector pairs (r ij , r kj ) and (r jk , r k ). The sets B b , B θ and B φ are the sets of pairs, triplets and quadruplets of bonded atoms along a polymer chain, all chains being considered. Therefore, card(B b ) = n pc (n P -1), card(B θ ) = n pc (n P -2) and card(B φ ) = n pc (n P -3), with n P the degree of polymerization. The harmonic potentials V b and V θ are defined as follows:

V b (r) = 1 2 K b (r -r 0 ) 2 , (4)
V θ (θ) = 1 2 K θ (cos(θ) -cos(θ 0 )) 2 , (5)
where r 0 and θ 0 are the equilibrium bond length and bond angle, and K b and K θ are model parameters. The torsional potential is given by

V φ (φ) = 5 i=1 A m cos m-1 (φ) , (6)
with {A m } 5 i=1 a set of additional model parameters. Nonbonded pair interactions are governed by a truncated Lennard-Jones 12-6 potential:

V LJ (r) = 4ε σ r 12 - σ r 6 , r r c (7)
= 0 , r > r c [START_REF] Berriot | Filler-elastomer interaction in model filled rubbers, a 1 H NMR study[END_REF] where the well depth ε and collision diameter σ are parameters depending on the polymer system, and r c denotes the cutoff radius (here, r c = 1.4 nm). The values of the force field parameters for the polymer matrix are listed in Tab. 1 (see [START_REF] Mayo | DREIDING: a generic force field for molecular simulations[END_REF][START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF]).

Potential

Initial configuration and relaxation steps

A set of 20 initial configurations is generated for three different systems, corresponding to n pc ∈ {10, 80, 320}, by using the Self Avoiding Random Walk (SARW) algorithm [START_REF] Binder | Monte Carlo and Molecular Dynamics Simulations in Polymer Science[END_REF]. In order to obtain fairly good initial configurational properties, MD relaxation steps are alternately performed in the NVT and NPT ensembles, making use of Nosé-Hoover schemes. More precisely, and following [START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF], simulations are started at 500 K, hence ensuring an amorphous structure through the induced high chain mobility. The latter is illustrated in Fig. 1, for n pc = 10. The systems are subsequently cooled down at 100 K, with a cooling rate of 1 K/ps. A NPT treatment is finally applied at 100 K and at a target isotropic pressure of 0 bar. It should be noticed that the simulations are performed below the glass transition temperature (which is estimated to 264 K for n pc = 10, following the methodology defined in [START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF] with a cooling rate of 0.1 K/ps), and that the above procedure yields samples exhibiting conformational characteristics in accordance with the literature.

Molecular model of Silica

In order to create amorphous SiO 2 nanoinclusions exhibiting an isotropic mechanical behavior, a set a α-quartz silica crystalline structures is first generated. These configurations are then relaxed by using MD simulations with Van der Waals and Coulomb (non-bonded) interactions. The first type of interactions is modeled by a Buckingham potential:

V BCK (r) = B 1 exp - r B 2 - B 3 r 6 (9)
where the values of parameters {B i } 3 i=1 are outlined below in Tab. 2 for the three types of bonds encountered in the system (note that the Van der Waals interaction between Si atoms is not taken into account). The Coulomb interaction is modeled by the following potential:

Bond B 1 (eV) B 2 ( Å) B 3 (eV /
V(r ij ) = q i q j 4π 0 r ij , (10)
where q i and q j are the electric charges of atoms i and j (q Si = 2.4e and q O = -1.2e, with e the charge of a proton), and 0 is the vacuum permittivity. Amorphous structures are then obtained through a heat up and cool down treatment, down to the equilibrium temperature of 100 K. In order to speed up computations, the above non-bonded interactions are subsequently replaced, at equilibrium, by fictitious bonded ones (see [START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF] for a discussion). More precisely, bond stretching and angle bending are modeled by using the harmonic potentiels defined in Eqs [START_REF] Holt | Dynamics at the Polymer/Nanoparticle Interface in Poly(2-vinylpyridine)/Silica Nanocomposites[END_REF][START_REF] Papon | Low-Field NMR Investigations of Nanocomposites: Polymer Dynamics and Network Effects[END_REF], and the values of associated parameters are listed in Tab the Si-O bond length and O-Si-O/Si-O-Si valence angles (at equilibrium) obtained from the atomistic simulations with either non-bonded or bonded interactions is shown in Fig. 2. It is seen that the two probability density functions are reasonably close to each other (note that the same conclusion can be drawn for the valence angles). The difference in variance is more significant for the Si-O bond length, for which the mean values for the non-bonded and bonded cases turn out to be pretty similar (with values equal to 1.6152 Å and 1.6247 Å, respectively). For each value of n pc , spherical nanoparticules of different radii are extracted from MD boxes following the procedure detailed in [START_REF] Brown | Effect of filler particle size on the properties of model nanocomposites[END_REF], so that the volume fraction of silica remains constant and equal to 4.8% regardless of the polymer system. Let R p denote the approximate particle radius. For n pc = 10 (resp. n pc = 80 and n pc = 320), it is found that R p ≈ 1.5 nm (resp. R p ≈ 3 nm and R p ≈ 4.8 nm). It is worth noticing that because of the atomistic model, each inclusion is not perfectly spherical and exhibits some roughness, especially for low target diameters. This characteristic, which is illustrated in Fig. 3 (where atomistic models and associated continuum approximations are represented for two different diameters), is an additional source of randomness within the atomistic description.

Molecular model of nanocomposite

Once the initial, equilibrated configurations of polymer and silica are obtained, nanocomposite samples are created as follows [START_REF] Brown | A molecular dynamics study of a model nanoparticle embedded in a polymer matrix[END_REF]. First, and for a given configuration of each subsystem, the MD box containing the polymer is swelled by a volume corresponding the one occupied by the particle. Second, a soft repulsive potential is gradually introduced, under NVT conditions, at the center of the box. Finally, the silica particle is inserted once the target hole diameter is reached, hence creating an initial configuration for the nanocomposite. This generation step is then followed by a MD relaxation procedure (similar to the one carried out on the polymer system, with a temperature increase up to 500 K, followed by a final decrease down to 100 K at a cooling rate of 1 K/ps) where only non-bonded Lennard-Jones interactions between the CH 2 sites and oxygen atoms are considered. The force field parameters for CH 2 -O interactions are estimated following the Lorentz-Berthelot mixing rules [START_REF] Allen | Computer Simulation of Liquids[END_REF]: σ = 0.31 nm, ε = 0.1288 kcal/mol. A final relaxation step under NPT conditions is performed during 2 ns at 100 K and 0 bar. A snapshot of one nanocomposite sample is shown in Fig. 4.

Morphological and mechanical characterizations 2.4.1. Radial distribution of density

Here, we investigate the radial evolution of polymer density in the nanocomposite system. To this aim, instantaneous estimations of the local density are obtained on spherical layers of constant volume V ρ and averaged over a sufficiently large time interval [0, ∆]. Each layer is specifically defined by an interior radius R i ρ -the exterior radius R e ρ being determined from the latter and

V ρ . Let R m ρ = (R i ρ + R e ρ )
/2 be the mean radius of a given layer. The instantaneous local density in the layer defined by mean radius R m ρ is defined as

ρ n (R m ρ , t) = M p N (R m ρ , t) N A V ρ , (11)
where M p is the molar mass of the polymer, N (R m ρ , t) is the number of CH 2 sites in the layer under consideration at time t and N A is the Avogadro constant. The radial density is then given by region hereafter. Note that such a definition for the interphase region is not unique and basically depends on the monitored quantities of interest (see [START_REF] Voyiatzis | How thick is the polymer interphase in nanocomposites? probing it by local stress anisotropy and gas solubility[END_REF] for numerical investigations -on a more detailed atomistic model and for another nanocomposite system -and insightful comments on this issue). Whereas this dependence may impact some quantitative results (such as the estimation of the correlation lengths) obtained in the remainder of this paper, it does affect neither the proposed probabilistic model nor the overall methodology related to the identification through a statistical inverse problem.

ρ n (R m ρ ) = 1 ∆ ∆ 0 ρ n (R m ρ , t) dt . (12)

Virtual mechanical testing in MD simulations

In order to obtain independent realizations of the apparent stiffness properties of the nanocomposite with n pc = 10 (see [START_REF] Shinoda | Rapid estimation of elastic constants by molecular dynamics simulation under constant stress[END_REF][START_REF] Marcadon | Confrontation between molecular dynamics and micromechanical approaches to investigate particle size effects on the mechanical behaviour of polymer nanocomposites[END_REF]), each configuration is submitted to a mechanical loading under NPT conditions, with a target pressure tensor of the following form:

P αβ (t) ij = p αβ (t) 2 (δ iα δ jβ + δ iβ δ jα ) , (13)
where α and β are fixed integers in {1, 2, 3}, and t → p αβ (t) is a stepwise function corresponding to a quasi-static discretization of a linear loading function (with a loading rate equal to 1 bar per 2 ps) between the initial and final pressures -hence limiting viscosity effects in the polymer system. For a given initial configuration, average deformations of the MD box are then extracted all along the loading and used to compute, by means of an ergodic estimator, the associated realization of the compliance tensor. Because of the roughness of small-diameter inclusions, six different loading conditions are sequentially applied in order to recover the whole slightly anisotropic tensor. These conditions correspond to the choices 1 α 3 and α β in Eq. ( 13). The deformations of the MD box for α = β = 1 and {α = 1, β = 2} are shown in Figs. 6 and7, respectively. The stress-strain curves for both the pure polymer and the nanocomposite system are shown in Fig. 8 (here, ergodic estimators for the components of the stress and strain tensors are computed at each load increment, hence resulting in smoothed stress-strain curves), for the above loading cases. As expected, it is seen that the nanocomposite system exhibits a stiffer mechanical behavior than the polymer matrix.

Stochastic modeling of the elasticity random field in the interphase region

In this section, we address the probabilistic modeling of the elasticity random field in the interphase region, which is denoted by D I . Let (O, e r , e θ , e ϕ ) be the system of spherical coordinates, and denote by x s = (r, θ, ϕ) any generic point of R 3 expressed in this coordinate system. Let x = (x 1 , x 2 , x 3 ) denote any point of R 3 given in the cartesian coordinate system, with x 1 = r sin(ϕ) cos(θ), x 2 = r sin(ϕ) sin(θ) and x 3 = r cos(ϕ). From the results of the MD simulations, it is postulated that:

• the interphase region can be represented as a spherical shell of constant thickness e I = 2 nm:

D I := {x s = (r, θ, ϕ)|r ∈ [R p , R p + e I ], θ ∈ [0, 2π], ϕ ∈ [0, π]} , (14)
with R p the radius of the particle under consideration (recall that R p is fixed for a given value of n pc )

• the material in the interphase exhibits a transversely isotropic local behavior defined by the unit normal vector n(x s ) = e r (x s ) at point x s ;

• the stiffness tensor random field exhibits a mean value that is independent of x s and a correlation structure with a spherical symmetry (in the spherical coordinate system);

• the probabilistic properties of the elasticity random field are similar in the directions defined by e θ and e ϕ .

It is further assumed that the mean model is independent of x s in spherical coordinates. Whereas such a modeling assumption could be readily relaxed by adding, for instance, a dependence along the radial direction, it is worth noticing that the model thus obtained would require additional data that are not considered here for model calibration.

Overview of the methodology

Let {[C s (x s )], x s ∈ D I } be the random field modeling the stiffness properties in the stochastic interphase. Let Ela T I (n) be the set of transversely isotropic tensors defined by the unit vector n.

The following algebraic decomposition is first considered:

∀x s ∈ D I , [C s (x s )] = η[I 6 ] + [M s (x s )] , (15)
where η 1 is an arbitrary small parameter, [I 6 ] is the 6-by-6 identity matrix and {[M s (x s )], x s ∈ D I } is a random field with values in Ela T I (n), with n = (1, 0, 0) in spherical coordinates. The first term in the right-hand side of Eq. ( 47) is introduced in order to preserve the uniform ellipticity of the stochastic boundary value problem. Next, for all x s fixed in D I , the random matrix [M s (x s )] is normalized according to

[M s (x s )] = [M s ] 1/2 [N s (x s )] [M s ] 1/2 , ∀x s ∈ D I , (16)
where

M s = E{[C s (x s )]} -η[I 6
] (note that the mean function depends on x in the cartesian system), and

{[N s (x s )], x s ∈ D I } is a random field such that E{[N s (x s )]} = [I 6 ], ∀x s ∈ D I .
It can then be shown (see [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF]) that there exists a unique random field

{[G s (x s )], x s ∈ D I } such that ∀x s ∈ D I , [N s (x s )] = exp{[G s (x s )]} (17)
with

[G s (x s )] = 5 i=1 G s i (x s )[E i ] , (18)
where {G s i (x s ), x s ∈ D I }, i ∈ {1, . . . , 5}, are non-Gaussian dependent R-valued random fields and {[E i ]} 5 i=1 denotes the matrix form of the so-called Walpole basis [START_REF] Walpole | Fourth-rank tensors of the thirty-two crystal classes: Multiplication tables[END_REF] of Ela T I (n), with n = (1, 0, 0). Let {G s (x s ), x s ∈ D I } be the R 5 -valued random field such that G s (x s ) = (G s 1 (x s ), . . . , G s 5 (x s )). Following [START_REF] Soize | Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF][START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF], the probabilistic model for random field {G(x s ), x s ∈ D I } is constructed by following a two-step strategy:

• first, the first-order marginal probability density function of the random field is constructed in the context of Information Theory [START_REF] Shannon | A mathematical theory of communication[END_REF] and more specifically, by invoking the principle of maximum entropy [START_REF] Shannon | A mathematical theory of communication[END_REF][START_REF] Jaynes | Information theory and statistical mechanics[END_REF] (with information implied by mathematical properties or inferred from atomistic simulations);

• second, the non-Gaussian vector-valued random field {G s (x s ), x s ∈ D I } is defined through a nonlinear mapping acting on an underlying centered second-order R 5 -valued Gaussian random field, the parametrized correlation structure of which is transported to {G s (x s ), x s ∈ D I }.

Construction of the family of first-order marginal probability distributions

The maximum entropy principle is formulated by considering a constraint on the mean value (raising from Eq. ( 16)),

E exp 5 i=1 G s i (x s )[E i ] = [I 6 ] , (19)
as well as a constraint related to the finiteness of second-order moments for random matrix [N s (x s )] [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF] (see [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]):

5 i=1 E {G s i (x s )} tr([E i ]) = ν , |ν| < +∞ , ( 20
)
where ν is independent of x s following the assumptions stated at the beginning of section 3. It can be shown that an approximate solution writes

p G s (x s ) (g) = c exp (-Φ(g)) , ( 21
)
where c is the normalization constant and the potential function Φ is defined as [START_REF] Staber | Approximate solutions of lagrange multipliers for informationtheoretic random field models[END_REF] Φ

(g) = α × δ -2 [N ] × tr exp 5 i=1 g i [E i ] - 5 i=1 g i [E i ] , ∀g ∈ R 5 . (22)
In Eq. ( 22), α is a parameter which depends on the symmetry under consideration (α ≈ 0.82 for the set Ela T I (n) of transversely isotropic tensors, regardless of the unit normal vector n), and δ [N ] is a scalar parameter independent of x s controlling the level of statistical fluctuations exhibited by [N s (x s )].

3.3. Definition and generator of random field {G s (x s ), x s ∈ D I } 3.3.1. Definition of a x s -dependent family of Wiener processes Let {Ξ(x s ) = (ξ 1 (x s ), . . . , ξ 5 (x s )) , x s ∈ D I } be the centered Gaussian R 5 -valued random field with independent components that are defined as follows. The diagonal entries of the continuous correlation function (x s , y s ) → [R Ξ (x s , y s )] with values in the set of diagonal real matrices are such that

[R Ξ (x s , y s )] ii = E{ξ i (y s )ξ i (x s )} =: ρ (x s , y s ) (23) with [R Ξ (x s , x s )] ii = 1 , 1 i 5 ( 24
)
for all (x s , y s ) in D I × D I . It follows that the correlation structure is completely defined by the normalized correlation function ρ. A separable structure is further assumed and written, for all x s = (r, θ, ϕ) and y s = (r , θ , ϕ ) in D I , as:

∀(x s , y s ) ∈ D I × D I , ρ (x s , y s ) := ρ r (τ r ) × ρ θ (r, r , τ θ ) × ρ ϕ (r, r , τ ϕ ) , (25)
where τ r := r -r , τ θ := θ -θ and τ ϕ := ϕ -ϕ are the radial and angular lags. Here, the above correlation functions are parametrized as:

                   ρ r (τ r ) := 2L r πτ r 2 sin 2 πτ r 2L r , ρ θ (r, r , τ θ ) := exp - 1 2 r + r L θ 2 sin 2 τ θ 2 , ρ ϕ (r, r , τ ϕ ) := exp - 1 2 r + r L ϕ 2 sin 2 τ ϕ 2 , (26)
where L r is the spatial correlation length along the radial direction and (L θ , L ϕ ) is a set of parameters controlling the decorrelation rates along the directions defined by e θ and e ϕ . For given values of r and r in [0, R p ], the functions ρ θ and ρ ϕ are 2π-periodic and satisfy:

• ρ θ (r, r , -τ θ ) = ρ θ (r, r , τ θ ) for all τ θ in [-2π, 2π]; • ρ ϕ (r, r , -τ ϕ ) = ρ ϕ (r, r , τ ϕ ) for all τ ϕ in [-π, π].
Note that the correlation functions ρ r , ρ θ and ρ ϕ could be parametrized by using other algebraic forms satisfying the conditions rising from the spherical description of the correlation structure.

Similarly, the angular correlation functions may be represented by means of expansions in Fourier series. However, the above forms interestingly depend on a small number of parameters and are thus well suited for an inverse identification from limited data. Let W = {W (t, x s ), t 0, x s ∈ D I } be the second-order centered Gaussian R 5 -valued random field defined as follows:

• ∀x s ∈ D I , W (0, x s ) = 0 almost surely;

• the generalized derivative D t W of W (with respect to t) is the normalized Gaussian cylindrical white noise B whose covariance function [C B ] writes, for all ∀(x s , y s ) ∈ D I × D I and for all τ ∈ R:

[C B (x s , y s , t + τ, t)] := δ 0 (τ )[R Ξ (x s , y s )] , (27)
with δ 0 the generalized Delta function at 0. For any x s fixed in D I , the stochastic process {W (t, x), t 0} is then a normalized R 5 -valued Wiener process (note that variable t does not have to be confused with time variable and in particular, with the time variable involved in the atomistic simulations).

Random generator

Let x s be fixed in D I . Let {(U (r, x s ), V (r, x s )), t 0} be the R 5 × R 5 -valued Markov process satisfying the following Itô stochastic differential equation [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF]:

dU (t, x s ) = V (t, x s ) dt dV (t, x s ) = -∇ u Φ(U (t, x s )) - 2 V (t, x s ) dt + √ η dW (t, x s ) , (28)
for all t ∈ R + , where Φ is the potential function defined by Eq. ( 22), ∈]0, +∞[ is a tunable parameter and {W (t, x s ), t 0} is the Wiener process defined above. Initial conditions are written as U (0, x s ) = U 0 (x s ) and V (0, x s ) = V 0 (x s ) almost surely, where the probability distribution of random variable (U 0 (x s ), V 0 (x s )) is assumed to given. It can then be shown that [START_REF] Soize | Construction of probability distributions in high dimension usign the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF] lim

t → +∞ U (t, x s ) law = G s (x s ) , (29)
hence providing a sampling algorithm for the non-Gaussian random field {G s (x s ), x s ∈ D I }. In this work, the stochastic differential equation defined by Eq. ( 28) is discretized by means of a Störmer-Verlet algorithm [START_REF] Verlet | Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules[END_REF][START_REF] Vogelaere | Methods of integration which preserve the contact transformation property of the Hamiltonian equations[END_REF][START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF][START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF] (see e.g. [START_REF] Kloeden | Numerical Solution of Stochastic Differentials Equations[END_REF][START_REF] Milstein | Stochastic Numerics for Mathematical Physics[END_REF][START_REF] Talay | Simulation and Numerical Analysis of Stochastic Differential Systems: a Review[END_REF] for thorough reviews on discretization schemes).

Inverse identification of the random field model from atomistic computations

Overall methodology

Let {[ C MD (ω i )]} 20 i=1 be the set of realizations for the apparent stiffness tensor estimated from MD simulations, which is denoted by [ C MD ] (see section 2.4.2). From the previous section, one has

∀x s ∈ D I , [C s (x s )] = η[I 6 ] + [M s ] 1/2 exp{[G s (x s )]} [M s ] 1/2 , ( 30
)
so that the probabilistic model for random field {C s (x s ), x s ∈ D I } depends on two vectors of hyperparameters, namely:

• a vector h 1 of hyperparameters gathering the coordinates of mean value [M s ] onto matrix basis

{[E i ]} 5 i=1 , with h 1 = (m 1 , . . . , m 5 ) ∈ S 1 ⊂ R 5 ;
• a vector h 2 = (δ [N ] , L r , L a ) ∈ S 2 of hyperparameters related to the level of fluctuations and to the correlation functions, where it is assumed (based from the MD results) that

L θ = L ϕ =: L a .
In this work, the optimal values of the above vectors, denoted by h opt 1 and h opt 2 respectively, are identified by stating the equivalence between the apparent properties obtained from MD simulations (with n pc = 10) and those determined by using standard computational homogenization in a continuum mechanics (CM) formulation. It should be noted that realizations of the local stiffness tensor field cannot be extracted from MD simulations and readily passed into the CM formulation (hence, a calibration based on a statistical inverse problem), for at least two reasons. First, the definition and CM-based interpretation of local properties extracted from MD results (see [START_REF] Papakonstantopoulos | Calculation of local mechanical properties of filled polymers[END_REF][START_REF] Riccardi | Molecular dynamics method to locally resolve poisson's ratio: Mechanical description of the solid-soft-matter interphase[END_REF]) is still an open question whenever confined matter is involved [START_REF] Marcadon | Confrontation between molecular dynamics and micromechanical approaches to investigate particle size effects on the mechanical behaviour of polymer nanocomposites[END_REF]. Second, the loading conditions applied in the atomistic virtual testing (where the local properties are triclinic with respect to the MD reference frame) are very different from those considered in the continuum setting, so that imposing a match, in some probabilistic sense, between the two local fields may not yield reliable results.

Let The properties of each deterministic isotropic phase were estimated from MD simulations and are as follows (k and µ denote the bulk and shear moduli, respectively):

• for the Silica inclusion: k Silica = 36.8, µ Silica = 30.5 (both in GPa).

• for the polymer bulk: k pbk = 5.19, µ pbk = 0.93 (both in GPa).

In the interphase region, the elastic properties are modeled by the random field model (detailed in section 3) that is described, from now on, by using the cartesian coordinate system -hence, the omission of superscript "s".

Homogenization problem for the continuum mechanics modeling

In order to formulate the statistical inverse problem, let us now introduce the upscaling procedure in the continuum mechanics framework. This procedure is based on averaging procedures that are widely used in homogenization theories. However, and as opposed to classical assumptions stated in such theories, it is worthwhile to notice that D is not assumed to be a representative volume element, so that upscaled properties are coined as apparent ones [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF] and exhibit non-negligible statistical fluctuations.

For selected values of h 1 and h 2 , let us consider the following stochastic boundary value problem (SBVP) [START_REF] Bornert | Homogenization in Mechanics of Materials[END_REF]:

-div [σ(x)] = 0 , ∀x ∈ D , (31)
[σ(x)] n(x) = [Σ] n(x) , ∀x ∈ ∂D , (32)
where x → [σ(x)] is the local stress field, n(x) is the unit vector normal to ∂D at point x, and [Σ] is the macroscopic stress tensor. The constitutive law writes:

[σ(x)] = C(x) : [ (x)] , ∀x ∈ D , (33)
where C(x) is the fourth-order tensor representation of elasticity matrix [C(x)]. Note that C(x) is random whenever x ∈ D I , and deterministic elsewhere (here, and for the sake of notational convenience, the same notation is used regardless of the deterministic or random nature of the elasticity tensor at a given point). It can then be shown that there exists a so-called stress concentration tensor field x → B(x) with values in the set of fourth-order tensors (with minor symmetries) such that

[σ(x)] = B(x) : [Σ] , ∀x ∈ D (34) and 1 |D|
D B(x) dx = I , (35)
with I the fourth-order symmetric identity tensor. By construction, one has

B ijk (x) = σ k ij (x) , (36)
where x → [σ k (x)] is the stress field solving the aforementioned stochastic boundary value problem with the macroscopic stress tensor Σ k given by:

Σ k ij = 1 2 (δ ik δ j + δ i δ jk ) . (37)
It follows that the random field x → B(x) can be recovered by solving the above stochastic boundary value problem for six different values of the macroscopic stress tensor Σ k , with k and in {1, 2, 3}. The apparent compliance tensor S defined with respect to the statically uniform boundary conditions is then defined as:

S = 1 |D| D S(x) : B(x) dx , (38)
where x → S(x) is the local fourth-order compliance tensor field. Let C = S -1 be the associated apparent stiffness tensor, and denote by [ C CM (h 1 , h 2 )] its second-order tensor representation (where the dependence on h 1 and h 2 is made explicit for subsequent use). The above apparent tensors depend on the applied boundary conditions (as discussed in the literature [START_REF] Ostoja-Starzewski | Microstructural Randomness and Scaling in Mechanics of Materials[END_REF]), which are chosen in accordance with the boundary conditions that are considered in the MD simulations (it is worth recalling here that the virtual testing is carried out by prescribing a target macroscopic pressure tensor in the NPT ensemble). Similarly, [ C CM (h 1 )] denotes the apparent tensor (under the same boundary conditions) obtained by considering an interphase with deterministic properties defined by vector h 1 in spherical coordinates.

In this work, the mechanical problem defined by Eqs. ( 31) and ( 34) is solved, for any realization of the random field in the interphase region, by the finite element method. Domain D is discretized with 4-node tetrahedron elements (with 1 Gauss integration point) -see Fig. 9. The final computational model is made up of 190, 310 elements (corresponding to 102, 561 degrees of freedom). It is worth noticing that such a mesh density provides at least four integration points per correlation length (in mean), regardless of the direction or the configuration tested while solving the statistical inverse problem, hence ensuring a good sampling of the correlation structure.

Remark on random field sampling by the ergodic method. Assume that each ISDE is solved by using a discretization depending on step ∆t for parameter t, and denote by t * := (k * -1)∆t the first instance of t within the stationary regime. For 1 i j 5, a decorrelation index τ ε ij (x s ) associated with stochastic processes {U i (t, x s ), t t * } and {U j (t, x s ), t t * } (where x s is any point in D I ) is then introduced and such that

E{U j ((k + τ ij )∆t, x s )U i (k∆t, x s )} ε (39)
for some arbitrary decorrelation threshold 0 < ε 1 and k k * . A global decorrelation index τ can be defined as

τ = sup x s ∈ D I max 1 i j 5 τ ε ij (x s ) . (40)
In practice, and for a sufficiently small value of ε 1, this result means that {U ((k+τ ε )∆t, x s ), x s ∈ D I } and {U (k∆t, x s ), x s ∈ D I }, k k * , are weakly correlated samples of {G s (x s ), x s ∈ D I }. This sampling strategy is especially relevant whenever the rate of convergence toward the stationary solution at all points in D I is low and requires a number of iterations much larger than τ ε , so that sampling by the Monte Carlo method would require a larger CPU time. The value of parameter ε must be carefully selected in order to balance sampling quality with computational cost.

Strategy for model calibration

The stochastic model is calibrated by solving a statistical inverse problem that is formulated in a two-step manner. First, the vector h opt 1 of optimal hyperparameters for the mean model is defined as

h opt 1 = arg min h 1 ∈ S 1 L 1 (h 1 ) , (41)
where L 1 is the cost function defined as

L 1 (h 1 ) := P iso (E{[ C MD ]}) -P iso ([ C CM (h 1 )]) F P iso (E{[ C MD ]}) F , (42)
with • F is the Frobenius norm and S 1 the admissible set such that the mean model is ensured to be positive-definite. In Eq. ( 42), P iso is the projection operator, defined in the sense of Euclidean norm, on the set of isotropic second-order tensors (see [START_REF] Moakher | The closest elastic tensor of arbitrary symmetry to an elasticity tensor of lower symmetry[END_REF] and the references therein). Note that this operator can be shown to commute with the operator of mathematical expectation. Second, the vector h opt 2 of optimal hyperparameters related to both the level of statistical fluctuations and the correlation structure is defined by the maximum likelihood method:

h opt 2 = arg max h 2 ∈ S 2 L 2 (h 2 ) , (43)
where L 2 the likelihood function given by

L 2 (h 2 ) = 20 i=1 p [ C CM (h opt 1 ,h 2 )] [ C MD (ω i )] (44)
and p [ C CM (h opt 1 ,h 2 )] is the probability density function of [ C CM (h opt 1 , h 2 )] estimated by the nonparametric statistics based on the use of the kernel estimation method.

Numerical results

Calibration of the mean model

The optimization defined by Eqs. [START_REF] Chang | Multiscale homogenization model for thermoelastic behavior of epoxy-based composites with polydisperse SiC nanoparticles[END_REF][START_REF] Ghasemi | Uncertainties propagation in metamodel-based probabilistic optimization of CNT/polymer composite structure using stochastic multi-scale modeling[END_REF] is solved by using the Matlab function fmincon (with default settings) under inequality constraints arising from the definition of S 1 . The convergence of the algorithm is shown in Fig. 10, where it is seen a good convergence is obtained after 40 iterations.

        (45)
in spherical coordinates (in GPa). A strengthening effect for the shear modulus is thus observed for the mean model in the interphase region (as compared to the bulk properties), which is in accordance with the loss of polymer-chains mobility in the neighborhood of the inclusion, as well as with numerical results obtained elsewhere (see e.g. [START_REF] Rahimi | Mechanical behavior and interphase structure in a silica-polystyrene nanocomposite under uniaxial deformation[END_REF]).

Calibration of the level of fluctuations and correlation structure

Here, the ISDEs are defined with = 9 and solved by using a step ∆t equal to 0.01 (note that these values were obtained from specific convergence analysis). For these parameters, the convergence toward the stationary solution is reached after 20, 000 steps. Furthermore, the parameter ε involved in the ergodic sampling is set to 0.1, which corresponds to a decorrelation index τ ε equal to 200. Below, all estimators are computed on a set of 300 realizations obtained from the ergodic sampling introduced at the end of section 4.1.1. The optimization problem defined by Eqs. [START_REF] Plimpton | Fast parallel algorithms for short-range molecular-dynamics[END_REF][START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito-the open visualization tool[END_REF] is solved by a regular sampling of the hyperparameters over their admissible sets. The correlation length in the radial direction is normalized by the interphase thickness e I , whereas the parameters related to angular decorrelation are normalized by the term 2πR ext , with R ext := R p + e I . The plot of the cost function L 2 is shown in Fig. 11. The optimal values are found to be δ opt

[N ] = 0.3, L opt r ≈ 0.9e I (= 1.8 nm) and L opt a ≈ 3.85 rad×nm, hence showing that the fluctuations in the interphase region are neither negligible nor prone to being smoothed out by an upscaling procedure (given that the condition L opt r e I is not satisfied for instance). The optimal value of parameter L a further shows that for a given radial position, the elasticity field is substantially correlated over a hemisphere (depending on the distance to particle surface). This result turns out to be consistent with MD results, as qualitatively illustrated in Fig. 12. As an illustration, a realization of the calibrated random field {C 11 (x), x ∈ ∂D ext I } (where ∂D ext I denotes the exterior boundary surface of D I ) is shown in Fig. 13. Finally, and from the perspective of validation, it is worth investigating whether the continuum can faithfully reproduce, in some sense, the data generated by the MD simulations. In order to get a physical insight on this issue, the realizations obtained from the MD simulations are plotted together with the estimates of (joint) probability density functions for some arbitrarily chosen components of the apparent tensor (defined by the calibrated random field model and through the SBVP stated in section 4.1.1) in Figs. 14 and15.

It is seen that for the selected components, all the realizations obtained by the MD simulations are associated with non-zero probability levels (note that this result holds for all remaining visualizable combinations of components), hence showing in part the consistency of both the probabilistic model and the identification strategy.

Calibration with MD finite sampling noise

In this section, the effect of MD finite sampling noise on calibration results is investigated. From a physical standpoint, it can be reasonably assumed that such a noise does not impact the correlation structure of the underlying elasticity random field in the interphase region, but may yield an overestimation of δ [N ] . Therefore, the sampling noise is taken into account below by modeling the spatially-constant but random properties of the bulk polymer and Silica inclusion by random matrices, and by considering a triclinic noise independent of x s in the interphase elasticity random field written in spherical coordinates. Let {[C s W N (x s )], x s ∈ D I } be the aforementioned random field in the interphase region. The same probabilistic model is used for all the random matrices and corresponds to the MaxEnt-based model proposed in [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]. For a given phase, the model depends on a prescribed mean value, as well as on a parameter controlling the level of statistical fluctuations. It is deduced from the atomistic simulations that the mean values for the polymer and Silica bulks are isotropic and defined by the properties listed at the end of section 4.1, whereas the mean for the random matrix modeling the noise in the interphase region is set to [I 6 ] (so that triclinic fluctuations are introduced). It is further assumed that the random matrices depend on the same dispersion parameter, denoted by δ noise . Optimal value for parameters δ [N ] and δ noise are then defined as those maximizing the likelihood function

L 3 δ [N ] , δ noise = 20 i=1 p [ C CM (δ [N ] ,δ noise )] [ C MD (ω i )] (46)
over [0, 7/11[×[0, 7/11[. In Eq. ( 46), [ C CM (δ [N ] , δ noise )] is the random variable modeling the apparent stiffness tensor obtained by homogenizing the random microstructure defined as follows:

• the random field {[C s W N (x s )], x s ∈ D I } in the interphase region is written as

∀x s ∈ D I , [C s W N (x s )] = η[I 6 ] + [M s (x s )] 1/2 [Z] [M s (x s )] 1/2 , (47)
where {[M s (x s )] , x s ∈ D I } is defined by the vector h opt 1 and parameters (L opt r , L opt a ) calibrated in section 4.4, as well as by δ [N ] (to be updated); [Z] is the M + 6 (R)-valued random matrix depending on δ noise (to be calibrated) and such that E([Z]) = [I 6 ];

• the random properties in the bulk polymer and Silica particle defined by the mean values computed in section 4.1 and dispersion parameter δ noise .

The plot of the likelihood function δ [N ] , δ noise → L 3 δ [N ] , δ noise is shown in Fig. 16. Optimal values are found to be δ opt [N ] = 0.2 and δ opt noise = 0.1, which shows that neglecting finite sampling noise for the MD system under consideration yields an overestimation of δ [N ] .

Conclusion

This study was devoted to the construction and identification, based on results arising from MD simulations, of a random field model for the elastic properties of the interphase region surrounding nanoscopic fillers. Atomistic simulations were performed on a prototypical nanocomposite and allowed some fundamental features of the random field to be inferred. A stochastic representation was then proposed, making use of a generalized formulation derived elsewhere. The stochastic field thus defined exhibits transverse isotropy in spherical coordinates (in accordance with the MD-predicted geometrical configuration of the polymer chain segments in the neighborhood of the nanoinclusion), and interestingly depends on a low-dimensional hyperparameter. An optimal value of this hyperparameter was then sought by solving a statistical inverse problem based on an upscaling procedure. It is shown that the elasticity random field exhibits nonnegligible fluctuations, and that the estimates of the spatial correlation lengths are actually consistent with characteristic lengths of the atomistic model. These results were finally refined by integrating fluctuations due to MD finite sampling noise, which allows for a better estimation of the parameter related to the level of fluctuations.

Figure 1 :

 1 Figure 1: Snapshot of the polymer chains (n pc = 10) during the MD relaxation step, starting from an initial configuration obtained by the SARW algorithm (from left to right).

Figure 2 :

 2 Figure 2: Probability density functions of Si-O bond length (left) and O-Si-O (black)/Si-O-Si (red) valence angles (right). Solid line: atomistic model with non-bonded interactions. Dashed line: surrogate atomistic model with bonded interactions.

Figure 3 :

 3 Figure 3: Visualization of nanoparticles: atomistic representations (left) and associated polyhedral representations (right). The smallest (resp. largest) particle corresponds to a diameter equal to 3 nm (resp. 6 nm).

Figure 4 :

 4 Figure 4: Snapshot of a nanocomposite sample at equilibrium (n pc = 10). The CH 2 sites appear as small particles; the Si and O atoms are the blue and red big particles, respectively.

Figure 5 :

 5 Figure 5: Radial density distribution of mass for three nanocomposite system, obtained from the nanoparticle surface and normalized by mass density of pure polymer ρ p corresponding.

Figure 6 :

 6 Figure 6: Deformation of the MD box for the loading t → [P 11 (t)] (tensile test): initial configuration (left) and configuration deformed at 10% (right).

Figure 7 :

 7 Figure 7: Deformation of the MD box for the loading t → [P 12 (t)] (shear loading): initial configuration (left) and configuration deformed at 5% (right).

Figure 8 :

 8 Figure 8: Stress-strain curves simulated for the pure polymer and the nanocomposite: tensile test (left) and shear loading (right).

  us now introduce the continuum model. Let D be the three-dimensional open bounded domain in R 3 such that D = (] -L D /2, L D /2[) 3 , with L D = 6.68 (nm), and denote by ∂D its boundary. Below, B O (r) denotes the ball of radius r (in nm), centered at the origin O of the cartesian reference frame (O, e 1 , e 2 , e 3 ). The Silica nanoparticle is then modeled by the ball D Silica := B O (1.5), whereas the interphase region is represented by domain D I := D ∩ (B O (3.5) \ B O (1.5)). It is seen that the diameter of the ball B O (3.5) is slightly larger than the characteristic length of D, but this configuration has a negligible effect in the computational results provided below. The domain occupied by the pure polymer bulk is thus D P := D \ (D Silica ∪ D I ).

Figure 9 :

 9 Figure 9: 3D views of the meshed continuum model (the inclusion appears in magenta, the interphase in white and the polymer matrix in turquoise).

Figure 10 :

 10 Figure 10: Convergence of the solver for the optimization problem defined by Eqs. (41-42).The solution [M s ] = 5 i=1 m opt i [E i ] (recall that h opt 1 = (m opt 1 , . . . , m opt 5 )) in matrix form is found to be

Figure 11 :

 11 Figure 11: 4D plot of the likelihood function h 2 → L 2 (h 2 ) for different combinations of δ [N ] , L r and L a . The values of the cost function are proportional to marker size, as well as to color intensity (infinite values of L 2 are arbitrarily replaced by the value -100 for visualization purposes).

Figure 12 :

 12 Figure 12: Visualization of an instantaneous configuration of polymer chains trapped into spherical shells of thickness 5 Å within the interphase region. The shell on left (resp. right) panel is the closest (resp. furthest) one to the surface of the nanoparticle.

Figure 13 :

 13 Figure 13: Plot of a realization for elasticity random field {C 11 (x), x ∈ ∂D ext I } (in GPa) on the exterior surface of the interphase region (with calibrated hyperparameters).

Figure 14 :

 14 Figure 14: Plot of the estimated probability density function for the continuum model (solid line) and associated MD-based realizations (red circles). Left: component C 11 . Right: component C 12 .

Figure 16 :

 16 Figure 16: Plot of the likelihood function δ [N ] , δ noise → L 3 δ [N ] , δ noise . Former results (with no sampling noise) correspond to the curve δ [N ] → L 3 δ [N ] , 0 .

Table 1 :

 1 Force field parameters used for the atomistic modeling of the polymer phase.

		Parameter	Value	Unit
	V b	K b	700	kcal/mol/ Å2
		b 0	0.15	nm
	V θ	K θ	124.2830	kcal/mol
		θ 0	112.813	deg
	V φ	A 1	2.1109	kcal/mol
		A 2	4.3229	kcal/mol
		A 3	1.1665	kcal/mol
		A 4	-7.6004	kcal/mol
		A 5	0	kcal/mol
	V LJ	σ	0.43	nm
		ε	0.1133	kcal/mol

Table 2 :

 2 Force field parameters used for the atomistic modeling of silica inclusion.

	Å6 )

Table 3 :

 3 Force field parameters used for the atomistic modeling of the silica particle.

	3. A comparison of the probability density functions of

Note that in practice, convergence analysis must be performed with respect to both V ρ and ∆. The graphs of (R m ρ -R p ) → ρ n (R m ρ )/ρ p is shown in Fig.5for different system sizes, where ρ p is the density of pure polymer estimated from MD simulations. Decreasing density oscillations are clearly observed in the neighborhood of the particle, no matter the value of n pc . These fluctuations occur for a distance to the particle surface less than 2 nm, and the density of pure polymer is recovered beyond this limit. It is worth pointing out that the main oscillations are still exhibited at higher temperatures, and that the noise in density estimation increases far from the silica particle. In this region, a tangential reorientation of the polymer segments and a decrease in chain mobility are further observed, hence making the polymer structure much less amorphous. The zone of matrix disturbance thus defined, with a constant thickness equal to 2 nm, is identified with the interphase

Acknowledgments

This work was supported by the French National Research Agency (ANR) under contract ANR-12-JS09-0001-01.


Thibault Lejemble
email: thibault.lejemble@ensimag.grenoble-inp.fr

Camille Schreck

Amélie Fondevilla

Pierre-Luc Manteaux

Nicolas Durin

Paul G Kry
email: kry@cs.mcgill.ca

Thibault Blanc-Beyne

Marie-Paule Cani

Interactive procedural simulation of paper tearing with sound

Keywords: CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation; tearing, paper, sound, real-time

: From left to right, an example result of our method showing detail added to the torn edge of the paper, a closer view of the geometry of the tear, and an illustration of an example application of our technique using interactive tearing paper for torn paper collage creation on a multi-touch tablet.

Introduction

Although it is a common material, paper is seldom represented in virtual environments. Indeed, paper has complex behaviors that are challenging to model and simulate, such as bending, crumpling, cutting, and tearing, while having resistance to shearing and stretching. There is a collection of successful techniques in recent research that can produce these behaviors through physics-based simulation methods, but not in real-time. In this paper we present the first tech-nique in computer graphics for real-time paper tearing with a fast phenomenological model for both geometry and sound.

Our goal is an interactive method that allows easy manipulation and tearing of virtual paper for games or artistic applications, such as the virtual crafts application proposed in Figure 1. Thus, in addition to the realism of the results, we must ensure real-time computation. Physical simulations are not well suited because the high rigidity required to accurately simulate paper often causes numerical instabilities. Additionally, the number of degrees of freedom necessary to obtain a sufficient level of detail leads to prohibitive computational costs. That is why we choose to use a phenomenological approach to procedurally reproduce the behavior of paper being torn.

For interaction, we focus on the simplified case shown in Figure 1 where the movement of each hand is defined by the index and thumb and stays in a 2D plane; although the sheet of paper is initially flat, rotation of the hands produces tearing and 3D deformations due to out of plane buckling.

Paper is composed of microscopic fibers which have an important influence on the behavior of paper at a macroscopic level. To reproduce the behavior, we identify and focus on three features in our procedural model: isometric deformations, stochastic tearing, and the spectral energy of tearing sounds. Isometry of paper deformation comes from the fact that the fibers are oriented in all directions and cannot stretch or compress. Thus, upon any deformation, the surface of the paper preserves distances according to its initial flat state (i.e., the surface of the paper is developable). The stochastic nature of the shape of a tear can similarly be attributed to the random distribution of fibers. While the main direction of tearing is imposed by the movement of the hands, details along the tear exhibit variation due to tear-growth favoring the direction of fibers at the crack tip. Finally, to have realistic virtual paper, it is important to also provide sound that matches the animation. The sound of paper tearing is caused by the breaking of fibers, and we note the distribution of fibers and speed of the tear has an important effect on the sound spectrum and the creation of a plausible auditory signal.

Our contributions are a real-time solution for the challenges of modeling and animating detailed tear shapes and paper buckling, and a model for generating plausible tearing sounds. The rest of this article is organized as follows. We first summarize the related work. Then we explain the method: the isometry-preserving model of deformation, the algorithm to find the path of the tear, and the model for the sound. Finally, we evaluate the results and discuss limitations and future work.

Related work

Modeling deformation and sound is a common objective in computer animation research with interactive applications such as games. In the context of tearing paper, we can divide the related work on deformation and tearing into several categories: physically based methods, geometric methods, and hybrid methods that mix physics and geometric approaches. Equivalently, related work on interactive sound synthesis tends to be either physics-or examplebased, or more commonly, a combination of both.

Deformation and tearing

The simulation of paper is similar to that of cloth, which is a topic that has been investigated extensively in computer graphics over the last three decades. In recent research, inextensible cloth models [START_REF] Goldenthal | Efficient simulation of inextensible cloth[END_REF]] and strain limiting techniques [START_REF] Thomaszewski | Continuum-based strain limiting[END_REF][START_REF] Wang | Multi-resolution isotropic strain limiting[END_REF]] have helped increase the quality of cloth simulations. While these techniques do not address the isometry of paper deformations, it is possible to explicitly take into account the developable nature of the surface through adaptive meshing, thus improving the simulation of physically based paper crumpling [START_REF] Schreck | Non-smooth developable geometry for interactively animating paper crumpling[END_REF]].

In computer animation, tearing of cloth and thin shells has received less attention in comparison to the fracture of solids. [START_REF] Metaaphanon | Simulation of tearing cloth with frayed edges[END_REF] present a method for simulating the tearing of cloth that also produces fraying at the edges. Although it makes sense to break cloth using the structure of fibers, using lower resolution cloth simulations and discrete cuts along existing edges allows for interactive simulation with lower fidelity results [START_REF] Souza | Fast simulation of cloth tearing[END_REF]. In general, these methods do not transfer well to the randomly distributed fibers in paper, where treating the surface as a continuum and remeshing to accommodate fracture is preferable. The method of [START_REF] Busaryev | Adaptive fracture simulation of multi-layered thin plates[END_REF] adapts the mesh to the fracture, and furthermore presents high-quality off-line simulations of tearing paper that include multiple layers. Including plasticity in the simulation allows for crumpling [START_REF] Narain | Folding and crumpling adaptive sheets[END_REF], and more recently, [START_REF] Pfaff | Adaptive tearing and cracking of thin sheets[END_REF] present adaptive meshing with both plastic deformation and fracture, permitting simulations of crumpling and tearing.

Physically based simulations produce excellent results, but not without large computational cost. In contrast, our approach to producing interactive deformations is a purely geometric approach and is therefore closer in spirit to kinematic methods. For instance, [START_REF] Sumner | Mesh-based inverse kinematics[END_REF] demonstrate pose changes for a cloth sheet using an example-based inverse kinematics method that exploits as-rigid-as-possible mesh deformation. Other work on geometric methods directly addresses developability in flexible surfaces. [START_REF] Rohmer | Folded paper geometry from 2D pattern and 3D contour[END_REF] use fold lines to direct procedurally generated static geometry of folded paper that approximately respects isometry. In contrast, [START_REF] Solomon | Flexible developable surfaces. Computer Graphics Forum (CGF)[END_REF] present an interactive developable surface modeling system that enforces exact developability.

Interactive sound synthesis

Early examples of sound simulation in computer animation use physics simulation and phenomenological impact and roughness models to drive modal vibration models [START_REF] Van Den Doel | Foleyautomatic: Physically-based sound effects for interactive simulation and animation[END_REF].

Modal vibration frequencies and decays can be estimated from data [START_REF] Pai | Scanning physical interaction behavior of 3d objects[END_REF], or computed using modal analysis of finite element models [O' Brien et al. 2002;[START_REF] Picard | Advances in modal analysis using a robust and multiscale method[END_REF]]. Accounting for acoustic transfer significantly increases the quality of modal sounds [START_REF] James | Precomputed acoustic transfer: Output-sensitive, accurate sound generation for geometrically complex vibration sources[END_REF], and similarly, propagation and scattering are important in larger scale virtual environments [START_REF] Manocha | Interactive sound rendering[END_REF]; however, we do not address these effects in our work. Furthermore, we note that tearing sounds are largely due to the breaking of fibers, and the resulting sounds are not well modeled by physically based modal vibrations.

In computer animation, the fracture of brittle elastic objects has been investigated with fast physically based sound synthesis methods that make use of simplified fragment vibration models and acceleration noise [START_REF] Zheng | Rigid-body fracture sound with precomputed soundbanks[END_REF][START_REF] Chadwick | Faster acceleration noise for multibody animations using precomputed soundbanks[END_REF]]. But tearing is quite different because the object being torn divides slowly into two parts, while methods for fracture sounds typically assume that the separation of fragments occurs at a much smaller time scale.

Data-driven techniques include methods that break up recorded sounds into small chunks, and recombine them to make new sounds [START_REF] Schwarz | Concatenative sound synthesis: The early years[END_REF]]. Such methods can be described as concatenative techniques, and they have been useful in the production of clothing sounds [START_REF] An | Motiondriven concatenative synthesis of cloth sounds[END_REF]. Example based techniques could be very useful in synthesizing paper tearing sounds, but ultimately we note that the much simpler approach of matching the noise spectrum of examples provides compelling and plausible results. We note that the acoustic emissions from fractures of paper and their characteristic power-law statistics has also been studied by the physics community [START_REF] Salminen | Acoustic emission from paper fracture[END_REF].

Paper deformation model

As mentioned above, simulating the behavior of a paper sheet under tearing remains a challenge for interactive computer graphics applications. Instead of using a mechanical simulation to reproduce this behavior, we choose to create a phenomenological model based on observations. Our model is able to reproduce plausible paper tearing in real-time, thus offering the possibility of interaction.

In order to capture the key features of paper deformation, we performed several informal experiments and observed the behavior of tearing paper. The setting of our experiments (and our model) is a paper sheet placed on a planar surface and manipulated with two hands. The thumb and index finger of each hand is in contact with the paper, and together they apply a rotation that deforms and tears the paper. An example observation can be seen in Figure 2a.

We observe that paper is inextensible, and note that this implies that it is developable, i.e., it can be unfolded onto a plane without being stretched or compressed. Therefore, as soon as paper begins to tear, it must also bend in order to preserve the distances. We also observe that paper stays smooth, and we assume that we can model the surface as a C 2 developable surface consisting of a set of parts: planes, cylinders, cones, and tangent developables. As the point at crack tip represents the one and only singular point of the surface in our model, we chose to define the bending portion of the paper as a generalized cone. As shown in Figure 2b, we thus distinguish three specific regions in our model: two planar regions, and a generalized cone. The cone is delimited by the two lines formed by the crack tip and each thumb, and defined by the set of lines through the apex s at the crack tip and the points of a C 2 curve, as seen in Figure 2c.

In our model, each hand is modeled by two points that represent the index finger and the thumb. They can be placed anywhere on the paper sheet. At each time step, the hands perform a rotation, which we denote ∆θ1 for the left hand and ∆θ2 for the right. This hand motion induces tearing (as described in Section 4) and a deformation of the surface between the thumbs. The planar regions follow a rigid motion directly imposed by the hands, rotating around the crack tip at each step by the corresponding angle of hand rotation.

The process of building the generalized cone is shown in Figure 3. We approximate the geodesic curve that joins the two thumbs by two Hermite cubic splines which meet at vertex m, at an equal distance to each thumb. The constraints to compute the Hermit curves are the position of each thumb, the point at the top of the bend m, and the tangents at those vertices. We choose all tangents to be parallel to the line between the thumbs to simulate the pressure applied by the thumbs on the paper. We select in advance an appropriate magnitude for the tangent to correspond to the size and type of paper. The vertical position of the middle point m is computed using an iterative algorithm such that the length of the curve between m and each thumb matches the distance between them in the initial 2D configuration. In order to preserve the distance between crack tip s and the points of the curve, we rotate the curve such that the distance between m and s on the 3D shape corresponds to the one on the 2D pattern. The generalized cone is then defined by the apex s and the curve between the thumbs. To find each vertex b on the curved border of the paper, we compute points that preserve the distance between s and b along the corresponding line.

Tearing model

Paper tearing in real-time is a challenging problem. In physicsbased simulation, the crack direction at each time step is based on the eigenvalues of the stress tensor at the crack tip, and is perpendicular to the greatest stress direction. In the case of simulated paper, an extremely fine discretization is necessary to produce realistic small scale details of torn paper. This makes stress-based crack direction computation inappropriate for real-time paper tearing.

Instead, we build a procedural tearing algorithm that exhibits the characteristic stochastic details of torn paper that we observe in empirical experiments. Our algorithm is based on two observations. First, the main crack direction is guided by the motion of the hands. Second, the details along the crack are due to the distribution of fibers and a resistance to stress that is not uniform across the paper's surface. Our algorithm is therefore made of two steps. We first compute a main tearing direction based on hand motion, and then add tearing details based on an anisotropic fiber distribution model.

Tearing direction

We observe that tearing initiates at a random position on the edge between the two index fingers. In our model, we therefore choose to initiate tearing at a random point along the edge of the sheet, specifically, at a point between where the rays from thumb to index positions intersect the edge. Then, at each time step i, we compute the position of the next crack tip si according to the motion of the hands. For this, we define a vector t that represents the advance of the tear between si and the previous crack tip position si-1, as shown in Figure 4. In our experiments, we observe that rotating a single hand produces a tear in the direction of the thumb of that hand. Therefore we compute t as a linear combination of the vectors v1 and v2 defined by the crack tip s and thumb of each hand. We weight each vector by the corresponding rotation angle ∆θ to compute the current tearing direction

t = C ∆θ1 v1 v1 + ∆θ2 v2 v2 , ( 1
)
where C is a constant that controls the tearing speed. That is, the norm t is the speed of the tear, and it is used to parametrize the sound as described in Section 5.

Tearing details

Advancing the tear using only the tear direction computed above will not produce the rough edges we see in real torn paper. Therefore we procedurally add fine details along the path between two successive crack tip positions. Our algorithm is inspired by the work of [START_REF] Chen | Physicsinspired adaptive fracture refinement[END_REF], and makes use of two 2D functions as shown in Figure 5. A Perlin noise function P lets us represent the stress resistance of fibers at each point of the paper, while a Gaussian function models the stress field Gi produced by hand motion between the steps i -1 and i. Specifically, the Gaussian is centered on the edge between si-1 and si. We combine these two functions to obtain a 2D texture νi(x, y) = (1 -p)P (x, y) + p/Gi(x, y),

where each texel can be considered as a breaking coefficient that intuitively represents the probability that the tear will go through this texel. The parameter p lets us adjust the blend between the two functions to achieve a desired rough appearance of the tear.

From this texture, the tear detail can be chosen as the shortest path, but we prefer a more intuitive and simple approach. Going from the texel representing si-1, we iteratively choose the neighboring texel with the greatest breaking coefficient until we reach the destination si. To ensure that we actually reach si without going backward, we only consider the texels in directions that form an angle less than π 2 with the vector between the current texel and the destination si.

Sound

Realistic auditory feedback plays an important part in the creation of compelling virtual experiences. We enhance our interactive animation by synthesizing plausible sounds for the virtual paper as it is torn. We choose to create the sound with a very simple and fast procedural algorithm that has minimal computational and memory costs. Our method is based on observations made from audio recordings of real paper, and hypotheses concerning the breaking of fibers and the speed of the tear.

The sound of tearing paper is caused by fibers or links between them being broken. With the random distribution of fibers in the paper, this process is largely stochastic, and therefore sounds like (see Figure 7). Indeed, white noise makes a good starting point for tearing sound synthesis given the largely flat nature of the spectral power density of the sounds we recorded. Each fiber produces an energy burst when being torn apart. As previously noted, the fiber distribution is anisotropic across the paper surface. This produces sudden changes of energy in the sound and resembles crackling, which is visible in the slow tearing example (blue) shown in Figure 6. We reproduce the phenomenon by first generating a white noise and then multiplying it by r α , where r is a random number between 0 and 1, and the exponent α is used to tune the sound to be something between white noise and a crackling noise. When α = 1 the sound is a white noise, while for a greater value of α the sound will contain clearer peaks of energy.

In our examples, we use α = 100.

The crack propagation speed also affects the sound. First, the sound energy changes according to the speed, as can be seen the comparison between a slow tearing and a fast tearing in Figure 6. This can be explained by the fact that a faster tear involves a larger number of broken fibers in the same amount of time. Second, the amplitude changes that are caused by varying density of fibers across the sheet will also occur faster. Thus, we linearly tune the amplitude and pitch of the modified white noise according to the speed.

Results and discussion

Our system is implemented in C++, using OpenGL for graphics and the FMOD library for sound synthesis. We use scripted motions to produce different real-time simulations of tearing paper with sound. An example of our results along with a close-up image showing de-tail along the tear can be seen in Figure 1, but please see the accompanying video for additional examples that also include synthesized tearing sounds.

We performed a simple experimental evaluation of our deformation model. Figure 9 shows a comparison of the cone height produced by our model with that of a real sheet of paper. Specifically, we compare the real and simulated cone height normalized by original distance between the thumbs at different amounts of strain (i.e., compression, normalized by the original distance between the thumbs). The results suggest that our selection of the magnitude of Hermite curve tangents in our model works well for the scale of paper and typical interactions (e.g., thumb placement as shown in the figure). We note that our model could be improved by making the tangent magnitude a function of the type of paper and the original distance between the thumbs.

In Figure 8, we show spectrograms and spectral power density plots of synthesized slow and fast tearing sounds. We are able to capture the noise and crackling features of tearing, but in comparison to recorded sounds, our synthesized sounds focus only on the tear and are therefore missing several important qualities related to the cone. Specifically, the cone acts as a resonator and its size and height influence the sounds created in real tearing. Further research is nec-essary to identify how the cone parameters can be included in our procedural sound synthesis method.

Conclusions and Future Work

In this paper we present a first geometric model for paper that can be torn in real time with integrated sound. The approach has several limitations that we would like to address in future work. Notably, the paper primarily lies in a plane (e.g., a virtual table top). While parts of the paper can leave the surface of the table during interactions, the placement and motions of the hands that drive the deformation and tearing of paper are also confined to the plane. We note that this constraint on interaction is not necessarily an important limitation in the context of how we envision our method used on multi-touch tablets or tables (e.g., Microsoft PixelSense). More importantly, in our current method the hands can only rotate about the crack tip, while translation must be addressed to have a general solution for table-top tearing. Furthermore, we only track the frontier of one tear, while many paper tearing scenarios involve the creation of multiple simultaneous tears and branching [START_REF] Manteaux | Interactive detailed cutting of thin sheets[END_REF]. We also note that we could take into account the weight of the paper, and how this influences the tear and the sound. Finally, we note that multiple layers of fibers tend to give a rough edge to torn paper and we leave the addition of these details to future work.

Figure 2 :Figure 3 :

 23 Figure 2: In (a), the model for the movement of the hands. In (b), the two planar regions in yellow and the generalized cone in red. In (c), an example of the cone shape formed between the tip of the tear s and the thumbs. The position and orientation of the two flat portions of the paper on either side of the cone are driven by the motion of the hands. The cone is defined by its apex s and two Hermit curves joining, respectively, each thumb and the vertex m at the middle of the geodesic passing though both thumbs.

Figure 4 :

 4 Figure 4: Our model for determining the direction of the tear at each time step depends on the incremental rotation of each hand.

Figure 5 :Figure 6 :

 56 Figure 5: Variation in the tearing direction is produced with the combination of 2D Perlin noise (left) with a 2D Gaussian function (right) aligned with the current tear-trajectory direction.

Figure 7 :

 7 Figure7: Spectrograms of 0.05 s clips from the middle of the recorded sounds shown in Figure6, with a plot of the spectral power densities. All frequencies are present at a relatively constant speed-dependent power, with a small frequency dependent power drop for fast tearing.

Figure 8 :

 8 Figure 8: Spectrograms of 0.05 s clips of synthesized sounds, with a plot of the spectral power densities.

Figure 9 :

 9 Figure 9: Experimental evaluation shows that our cone deformation model produces results that correspond to real paper.

Acknowledgements

This work was partially funded by the advanced grant no. 291184 EXPRESSIVE from the European Research Council (ERC-2011-ADG 20110209).