Characterization of the material behavior and identification of effective elastic moduli based on molecular dynamics simulations of coarse-grained silica

The addition of fillers can significantly improve the mechanical behavior of polymers. The responsible mechanisms at the molecular level can be well assessed by particle-based simulation techniques, such as molecular dynamics. However, the high computational cost of these simulations prevents the study of macroscopic samples. Continuum-based approaches, particularly micromechanics, offer a more efficient alternative but require precise constitutive models for all constituents, which are usually unavailable at these small length scales. In this contribution, we derive a molecular-dynamics-informed constitutive law by employing a characterization strategy introduced in a previous publication. We choose silicon dioxide (silica) as an exemplary filler material used in polymer composites and perform uniaxial and shear deformation tests with molecular dynamics. The material exhibits elastoplastic behavior with a pronounced anisotropy. Based on the pseudo-experimental data, we calibrate an anisotropic elastic constitutive law and reproduce the material response for small strains accurately. The study validates the characterization strategy that facilitates the calibration of constitutive laws from molecular dynamics simulations. Furthermore, the obtained material model for coarse-grained silica forms the basis for future continuum-based investigations of polymer nanocomposites. In general, the presented transition from a fine-scale particle model to a coarse and computationally efficient continuum description adds to the body of knowledge of molecular science as well as the engineering community.


Introduction
The versatility of polymers can be further expanded by the addition of fillers. In particular, filler particles at the nanoscale have been found to significantly improve the polymers' mechanical [1,2], thermal [3,4], electrical [5], and optical [6] behavior. Experimental studies show that the presence of filler particles leads to a significant modification of the properties of the surrounding polymer [7][8][9]. This region of influence, referred to as the interphase, is considered to be responsible for the improved behavior of nanocomposites [10]. The experimental investigation of these nanocomposites is highly challenging due to the extremely small length scales and is further complicated by the large number of influencing factors [11] such as matrix-filler interaction [12] as well as filler size [13], shape [14], distribution [15], and content [16]. In order to identify promising polymer composite configurations more easily, empirical models and virtual testing strategies are increasingly resorted to.

Empirical models
Initially, attempts were made to estimate the properties of the polymer composites by empirically motivated, mostly analytical models, e.g., by Leidner and Woodhams [17] and Pukanszky [18]. These models, originally introduced for conventional composites, were extended to cover nanocomposites as well [19,20]. With these extensions, mechanical properties such as Young's modulus [21] or tensile strength [22] can be derived from the material behavior of the individual constituents. However, the complex mechanisms within the interphases responsible for the observed reinforcements can only be represented in a highly simplified way in such empirical models, e.g., by the B and Z parameters in the work by Peng et al. [23].

Molecular dynamics studies
By way of contrast, particle-based numerical approaches such as molecular dynamics (MD) inherently capture the effects within the interphase since they are solely based on the particle interactions. Therefore, they are widely used to study nanocomposites: Lu et al. [24] investigate graphene-polymer nanocomposites by means of atomistic MD simulations and identify an interphase that exhibits 1.5 times the polymer matrix stiffness. In previous publications [25,26], we have studied the interphase effects in polystyrene-silica composites using an MD-continuum coupling method [27][28][29][30]. This procedure enabled us to perform comparative pseudo-experiments on a large set of samples to identify elastoplastic effects occuring in the interphase [26]. Bocharova et al. [31] use coarse-grained (CG) MD simulations to complement their experimental findings and unravel the molecular mechanism responsible for observed reinforcements in poly(2-vinylpyridine)-silica nanocomposites. Despite these essential findings, MD simulations remain computationally too expensive to investigate engineering problems on macroscopic samples.

Micromechanics
To reduce the computational costs, micromechanics rely on efficient continuum material descriptions, while still resolving the composite at the level of the nanoparticles. To this end, representative volume elements (RVEs) that include all relevant geometrical information are considered and their mechanical behavior is homogenized to obtain the macroscopic properties of the composite [32]. However, the continuum constitutive laws of the constituents (filler and polymer matrix) necessary for the RVEs are usually not available at these extremely small length scales.
Therefore, the scope of this contribution is to present the calibration of a continuum mechanical constitutive law based on MD simulations of CG silica serving as an exemplary filler material (section 2.1). We employ a strategy introduced in our previous contribution [33] to characterize the material behavior using MD pseudo-experiments (section 2.2). To this end, we perform uniaxial tensile and simple shear tests (section 2.3) with time-proportional and time-periodic loadings. Within the limits of the small strains expected in a polymer nanocomposite, the examined silica exhibits linear elasticity with a pronounced anisotropy (sections 3.1 and 3.2). Based on these pseudo-experimental data, we calibrate linear elastic constitutive laws with different degrees of anisotropy, and thus accurately represent the material behavior (section 3.3).
This study is complementary to the characterization [33,34], and complex continuum modeling of polystyrene using a novel viscoelastic-viscoplastic constitutive law [35].

Coarse-grained silica
In this study, we rely on CG force fields introduced by Ghanbari et al. [36], who substitute each SiO 2 unit by one CG bead with an aggregated molar mass of 64.0 g mol21. The CG interaction potentials were derived via iterative Boltzmann inversion [37] and comprise bond, angle, and nonbonded contributions. The obtained force fields were used and validated in the context of silica-enforced polystyrene in numerous studies [25,26,[38][39][40]. The CG beads form a crystalline structure due to the applied physically motivated force fields [36]. We obtain a triclinic unit cell with the lattice parameters a = 7:8Å, b = 19:1Å, and c = 5:3Å and the corresponding angles a = 116:08, b = 97:48, and g = 71:98. Since we are interested in the overall behavior of the CG silica, we are not considering finite spherical filler particles as in previous studies [26], but an infinite sample under periodic boundary conditions (PBC). Due to the triclinic structure, we set up a triclinic periodic simulation cell containing 1350 CG beads to ensure a box length of at least twice the cut-off radius in a single crystal arrangement, as shown in Figure 1.
Due to the regular structure of the single crystal, the resulting sample is not subjected to stochastic variations. Therefore, a statistical validation using different initial configurations, as carried out in comparable studies of polymers [33,34,41], is not necessary. After the initial placement of the beads, the sample undergoes an equilibration procedure (cf. section ''Sample preparation details'' in Appendix 1), resulting in a mass density of 2.32 g cm 23 which fits well with the literature [42].

Material classification scheme and characterization strategy
The mechanical response of a material to an external load can be assigned to one of the four paradigmatic material models: elasticity, plasticity, viscoelasticity, and viscoplasticity [43]. Merely, the presence of rate dependence and quasi-static hysteresis are decisive factors [44] as shown in Figure 2.
A practical approach to apply this classification scheme based on MD simulations and to quantify the individual contributions was presented as a characterization strategy in the work by Ries et al. [33]: First, time-proportional loading tests with different strain rates are performed to identify whether the material exhibits a rate dependence, i.e., viscous behavior. Second, time-periodic tests reveal whether the material is reversible by means of a quasi-static hysteresis. The application of quasi-static loads is not possible with MD due to the required computational capacities. However, we were able to show in the previous works [33,34] that the occurrence of a stable equilibrium hysteresis as a result of several loading-unloading cycles permits to evaluate the dissipation and thus the inelasticity of the material. Based on the obtained results, relaxation and creep tests may be carried out to decompose the inelasticity into its viscous and plastic components.

Applied deformation
In order to fully capture the constitutive behavior of the investigated material, both uniaxial and shear deformation, as depicted in Figure 3, need to be considered. In the case of an anisotropic, i.e., directiondependent, material response, these principal load cases have to be applied for all three spatial directions. The silica considered here is typically embedded into a softer polymer matrix, e.g. in the form of spherical particles [11]. Due to the significant difference in the stiffness of the two phases, it can be assumed that the silica is only undergoing small deformations. In this case, linearized kinematics can be used, described by the Cauchy strain tensor relying on the assumptions of the infinitesimal strain theory [44].

Uniaxial tension.
For uniaxial tension the tensile strain e ii is given by with displacement u i , initial and deformed box length in tensile direction L i and l i , respectively. In the lateral directions, the simulation cell is subjected to the atmospheric pressure p atm , as shown in Figure  3(a). This setup allows for unconstrained lateral contraction, comparable to real experiments. In MD, this is realized by combining the isothermal-isobaric ensemble (NPT) in lateral directions with a predefined deformation in tensile direction, which can be interpreted as a one-dimensional canonical ensemble (NVT) [33].

Simple shear.
Simple shear, on the other hand, is induced by the parallel shifting of planes proportional to their position and is, therefore, by definition isochoric [46]. The associated shear strain e ij is related to the displacement u i and the inital orthogonal box length L j , or the shear angle g ij via The employed parallelepiped-shaped simulation cell necessitates the use of triclinic boxes (cf. Figure  3(b)). These feature additional angular degrees of freedom g xy , g xz , and g yz , compared to conventional orthogonal cells. This simplifies the application of simple shear in MD: The requested shear angle is controlled directly according to equation (3), while all other degrees of freedom are fixed via an NVT ensemble.

Simulation details.
We choose a Nose´-Hoover thermostat [47] and barostat [48] with a temperature and pressure coupling time of 500 fs and 500 fs, respectively. Contrary to our previous studies [33,34,49], the time step is kept constant at Dt = 5 fs. We apply the deformation continuously over the evaluated time t either in time-proportional or time-periodic manner given by respectively. The former represents a linear ramp-up with a constant strain rate _ e, while the latter follows a sinusoidal function with strain amplitude e a and maximum strain rate _ e max . Based on the findings of previous studies [33,34], we evaluate four load cycles for the time-periodic tests.

Stress measures
We evaluate the internal stress due to the external load using the virial stress tensor [50] with volume V , number of atoms n A , atom position vectors r I and r J , atom velocity y I , and atom mass m I . The pair-wise interaction force f IJ is obtained, as usual in MD, by derivation of the interaction potentials. Note that for the present systems, the stress is dominated by the virial components (minuend), while the kinetic contributions (subtrahend) only represent the samples' thermal fluctuations, as demonstrated in section ''Stress evaluation'' in Appendix 1. According to Subramaniyan and Sun [51], the time-averaged virial stress coincides with the Cauchy stress commonly used in continuum mechanics. This hypothesis is, however, controversially discussed [52,53]. Since the stress obtained in MD simulations is usually subject to strong fluctuations, a Savitzky-Golay filter [54] is used for clearer visualization in all figures of section 3.

Time-proportional loading
The approximately linear stress-strain curves due to uniaxial tension in xx, yy, and zz directions are shown in Figure 4(a). The significantly softer material behavior in the yy direction indicates anisotropy, as expected due to the crystalline microstructure. Furthermore, the investigated silica exhibits strain rate independence and thus does not feature viscous components. Very similar results in terms of viscosity and anisotropy are observed for the simple shear deformation in Figure 4(b). Initially, we observe a slightly stiffer material response in the yz direction. However, at strains larger than 1.8 %, substantial instantaneous stress drops occur that are caused by the sudden transition (''snapping'') of the crystal lattice to an energetically more favorable state. Since our sample comprises only a single crystal, these stress drops are particularly pronounced. Since we can already exclude viscosity at this point, the material is either purely elastic or elastoplastic (cf. Figure 2).
By evaluating the slope in the stress-strain diagram Figure 4(a), the effective moduli E x = 198 GPa, E y = 142 GPa, and E z = 200 GPa are identified. These quantities do not yet represent a comprehensive constitutive description as usually determined from experiments. In comparison with the reported range of 46-92 GPa [55,56], it is evident that our silica behaves distinctly stiffer. This is probably due to the shift in dynamics induced by coarse graining [57] and emphasizes the importance of a thorough investigation of the mechanical properties of MD models. A similarly high Young's modulus was reported by Yu et al. [58] for silica nanobeams studied with MD simulations.

Time-periodic loading
The stress response induced by sinusoidal uniaxial tension-compression loading with a strain amplitude of 4% is plotted in Figure 5. For deformation in the xx and yy direction, we observe stress-strain hystereses, induced under compression and disappearing with an increasing tensile load. Interestingly, there are equivalent paths for this inelastic behavior evident in Figure 5(a): While in the first and fourth load cycle, the stress already returns to the reversible path in the compression range; in the second and third cycle, this occurs only in the tension range. This is observed analogously for deformation in yy direction (cf. Figure 5(b)). The area enclosed by the stress-strain curves in cyclic tests is equivalent to the dissipated energy. Since dissipation only occurs for inelastic behavior and viscosity has already been excluded, the discussed hysteresis indicates plastic effects. In contrast, the zz direction displays no hysteresis and thus no plasticity, evident in cf. Figure 5(c). As shown by the insets in Figure 5(a)-(c), plasticity does not yet occur at a strain amplitude of 1.5 %. Analogously, the stress-strain curves for sinusoidal simple shear are combined in Figure 6. In contrast to uniaxial deformation, hystereses occur only in the tensile region for deformation in xz and yz direction. The underlying plasticity is attributed to the abovementioned snapping of the crystal lattice (cf. section 3.1). No inelastic effects occur for deformation in xy direction (cf. Figure 6(a)) or for strain amplitudes smaller than e a = 1:5 %, as shown in the insets in Figure 6(a)-(c).
In summary, the investigated silica can be characterized as elastoplastic with pronounced anisotropy. However, there is a purely elastic range up to 1.5% strain, which is limited by the lattice snapping under shear deformation.

Calibration of the constitutive model
In the nanocomposite, the nanoparticles are randomly oriented, so the overall directional dependence is expected to disappear. For a precise continuum mechanical model of the neat filler, the observed anisotropy nevertheless needs to be accounted for. When used as a filler for nanocomposites, silica is embedded into a softer polymer matrix, e.g., polystyrene, and therefore experiences comparatively small deformations. Therefore, we aim to accurately model the observed material behavior in the elastic range up to 1.5 % strain by means of a constitutive law based on the generalized Hooke's law, cf. equation (6)  in section ''Elastic material models'' in Appendix 1. To this end, optimal elastic constants have to be identified for the elasticity matrix C (7). For this optimization problem, the complete deformation state obtained in terms of e(t) from the MD simulations with time-proportional loading serves as input while the corresponding stress state s(t), i.e., all stress entries are used as goal function. The deformation constraint thus does not coincide with the MD boundary conditions, which seems counterintuitive, but is admissible due to the stress-strain duality of pure elasticity. Taking into account only the tensile range effectively implies tension-compression symmetry, which is, however, supported by our pseudoexperimental data in the range of 61:5 %, cf. insets of Figures 5 and 6. A gradient-based, steepest-descent optimization [59] (lsqcurvefit) is used to identify a set of parameters capable of representing all six load cases simultaneously. To facilitate the optimization, four support points at e = ½0, 0:5, 1:0, 1:5 t are considered per load case, which is sufficient due to the approximately linear curves.
We introduce D as the root-mean-square deviation (RMSD) of the stress averaged over all support points. As expected, large deviations of D iso = 328 MPa occur for an isotropic constitutive law. A significant reduction is achieved with transverse isotropy in xy, xz, and yz direction where we obtain D tran, xy = 168 MPa, D tran, xz = 166 MPa, and D tran, yz = 179 MPa, respectively. However, the fact that all transverse orientations yield similar results contradicts transversal isotropy.
Considering an orthothropic or a fully anisotropic constitutive law further decreases the deviation to D orth = 156 MPa and D aniso = 99 MPa, respectively. We compare their performance in detail by means of the stress-strain diagrams for uniaxial tension and simple shear in Figures 7 and 8. Both approaches are capable of predicting the important stress in the loading direction accurately. For uniaxial tension, deviations are observed mainly for shear stresses, which partly deviate distinctly from the goal of 0 MPa. For simple shear, the fit for both constitutive laws performs very well. The only exception is s xz for deformation in xy direction (cf. Figure 8(a)), which cannot be captured at all in the orthotropic case and only partially in the fully anisotropic model. Evidently, the fully anisotropic elasticity matrix is better suited for modeling the observed stress response, but requires 21 elastic constants, and thus significantly more than the orthotropic model with only nine parameters. The resulting optimal elastic coefficients are given in equations (8) and (9), respectively. In crystalline materials, the unit cell defines the structure of the crystal lattice and thus its symmetries. These symmetries are also present in the overall elastic properties for a monocrystalline sample, and thus the degree of anisotropy depends directly on the unit cell structure. In our case, the unit cell is triclinic, i.e., fully asymmetric, which is consistent with the sample's fully anisotropic elasticity.

Conclusion and outlook
In this paper, we have approximated the mechanical behavior of silica at the nanoscale by a continuum mechanical constitutive law. For this purpose, a CG sample was examined using MD deformation tests under PBC. Based on uniaxial and shear tests, time-proportional as well as time-periodic, we were able to characterize the present material as elastoplastic with pronounced anisotropy. For the range of small strains up to 1.5%, which is relevant in the context of polymer nanocomposites, we have identified the elastic constants of the generalized Hooke model. These can accurately predict the material behavior in this range. This study demonstrates that a constitutive law on the nanoscale, as it is necessary for the modeling of an RVE in micromechanics, can be derived on the basis of MD simulations. This contribution is complementary to an earlier characterization and modeling of polystyrene as an exemplary matrix material [33][34][35]. Together with our extensive efforts to characterize the interphase between polymer and filler particles [25,26], we can now provide continuum mechanical descriptions of all constituents of a nanocomposite. This enables us to set up RVEs to apply the molecular-level findings on the microscale. familiar temperatures obtained in atomistic simulations. A final relaxation run at T = 100 K for 40 ns ensures convergence to an energetic minimum (cf. Figure 9(c)). We choose a Nose´-Hoover thermostat [47] and barostat [48] with a temperature and pressure coupling time of 500 fs and 5000 fs, respectively, and a time step Dt = 5 fs.

Stress evaluation
To verify our stress analysis, we investigate the influence of the kinetic component of the virial stress tensor in equation (5). To this end, Figure 10 displays the stress-strain curves of simple shear simulations in xy, xz, and yz directions for the largest strain rate of 10% ns -1 . It is evident that the only virial curves coincide with the virial and kinetic curves which are filtered as described in section 2.4. Thus, the observed fluctuations are attributed to the kinetic contribution and thus are merely thermal noise. This shows that even for the largest shear rates investigated in this paper, no significant background streaming velocity occurs and thus verifies the stress evaluation employed here. (c) (a) (b) Figure 9. Equilibration: time evolution of (a) temperature T, (b) density r, and (c) potential energy of the silica sample.