Nano-filled polymers are becoming more and more important to meet the continuously growing requirements of modern engineering problems. The investigation of these composite materials at the molecular level, however, is either prohibitively expensive or just impossible. Multiscale approaches offer an elegant way to analyze such nanocomposites by significantly reducing computational costs compared to fully molecular simulations. When coupling different time and length scales, however, it is in particular important to ensure that the same material description is applied at each level of resolution. The Capriccio method [ 5], for instance, couples a particle domain modeled with molecular dynamics (MD) with a finite element based continuum description and has been used i.a. to investigate the effects of nano-sized silica additives embedded in atactic polystyrene (PS) [ 6,4]. However, a simple hyperelastic constitutive law is used so far for the continuum description which is not capable to fully match the behavior of the particle domain. To overcome this issue and to enable further optimization of the coupling scheme, the material model used for the continuum should be derived directly from pure MD simulations under thermodynamic conditions identical to those used by the Capriccio method. To this end, we analyze the material response of pure PS under uniaxial deformation using strain-controlled MD simulations [ 10]. Analogously, we perform simulations under pure shear deformation to obtain a comprehensive understanding of the material behavior. As a result, the present PS shows viscoelastic characteristics for small strains, whereas viscoplasticity is observed for larger deformations. The insights gained and data generated are used to select a suitable material model whose parameters have to be identified in a subsequent parameter optimization. Keywords molecular dynamics ·simulation of polymers ·mechanical properties of polymers ·material characterization 1 Introduction and methodology In contrast to continuum mechanical approaches, particle-based simulation methods like molecular dynamics (MD) enable the investigation of mechanisms taking place on the molecular level. However, due to the high computational costs of such methods only small systems can be examined. To overcome this, partitioned-domain multiscale techniques were developed, which allow to simulate regions of high interest on a fine resolution while treating the remaining regions on a coarse scale to reduce the overall computational effort. One example is the Capriccio method [ 5] that couples a particle domain modeled with MD with a finite element based continuum. This approach was already successfully employed to investigate the effects of nano-sized silica additives embedded in atactic polystyrene by [6] and [4]. Characterization of polystyrene using MD PREPRINT When a material is treated at different resolutions it is necessary that the individual material descriptions fully match. The material behavior on the fine scale is predefined by the interactions between the particles usually given as force fields. Consequently, the constitutive law used for the continuum has to be identified and parameterized accordingly. To this end, we introduced a methodology to characterize the mechanical behavior of a material directly from MD simulations [ 10] based on a classification scheme by [ 2] depicted in Tab. 1. We have derived four load cases from this, Table 1: Classification scheme to characterize material behavior introduced by [2] (visualization: [10]) rate dep endence y es noquasi-static hysteresisy es novisco- plasticit yplasticit y visco- elasticit yelasticit y visualized in Fig 1: time proportional loading with different strain rates to evaluate the rate dependence, time periodic tests either with varying strain amplitude or varying strain rate to quantify the inelastic effects and relaxation tests with preceding time periodic loading to identify whether a quasi-static hysteresis exists. (a) strain rate timestraintime p rop o rtional (b) timestraintime p erio dic (c) timestraintime p erio dic (d) timestrainrelaxation Figure 1: Applied deformation: (a) time proportional with varying strain rate, (b) time periodic with different strain amplitude, (c) time periodic with different maximum strain rate, (d) relaxation tests with different preceding strain amplitudes Up to now, we applied this scheme to investigate atactic polystyrene as a model system under uniaxial deformation and identified purely elastic behavior for small strains, followed by a viscoplastic regime for larger deformations [10]. In this contribution, we supplement our work with examining the same material under pure shear deformation: First, we show how this deformation can be applied within the MD framework, then we present and discuss the results obtained from the individual load cases. The acquired insights and the generated data are used by [ 11] for the equilibration of a suitable constitutive law. Fig. 2 illustrates the procedure needed to use multiscale modeling techniques to gain deeper insights into polymer nanocomposites or fracture of polymers and puts the scope of this contribution into context. 2 Preparation of MD systems As a model system atactic PS is investigated in coarse-grained resolution in which each monomer is substituted by a so-called superatom. As a result, the number of particles and thus the computational effort is reduced by a factor of 16 compared to an atomistic model. This enables us to examine samples consisting of a sufficiently large number of polymer chains at reasonable computation times. On the contrary, the reduced number of particles leads to an implicit speed-up of the dynamics of a coarse-grained model and thus all subsequently derived quantities evaluated on the coarse-grained level do not match their atomistic counterparts. However, this does not affect the significance of our results, as our aim is to reproduce the behavior of the coarse-grained particle model with a continuum description. 2 Characterization of polystyrene using MD PREPRINT sepa ration of domains enablesenablesp rovidesequili- b rates(a) (b) (c) (d)P a rticle-Based material description rECha racterization F strainstressContinuum mech. constitutive la w Ψ(C) = ... P a rtitioned-Domain multiscale mo deling Figure 2: Schematic overview: Particle-based material description (a), characterization procedure (b), continuum mechanical constitutive law (c) and partitioned-domain multiscale investigations (d), [10] The coarse-grained force fields describing the interactions between the superatoms were derived by [ 1] applying an iterative Boltzmann inversion approach [9] presented in [7]. The PS chains are placed in the simulation box via a self-avoiding random-walks algorithm implemented by [ 1]. In our case, these MD systems consist of 300 PS chains with 200 superatoms each and thus 180 000 degrees of freedom have to be considered. Subsequently, the MD samples are equilibrated with constant number of particles in three steps: Firstly, the molecules are given 5 ns to disentangle at constant volume and a temperature of 590 K (NVT ensemble) followed by an equilibration of20 ns under constant temperature and atmospheric pressure of 101.3 kPa (NPT ensemble). Secondly, the samples are cooled down to 100 K at a constant cooling rate of 5 K ns−1. This is significantly below the glass transition temperature of170 K of coarse-grained PS [ 8]. Lastly, the systems are kept at constant temperature and pressure for another 2 ns to ensure sufficient equilibration at low temperature. These steps are carried out under periodic boundary conditions with the MD solver IBIsCO developed by [ 3]. In order to obtain statistically reliable results, 10 slightly different samples are generated in this way. A more detailed description of the initialization and equilibration of the MD samples can be found in [10]. 3 Pure shear deformation within MD Embedding initial and deformed configurations of our MD sample in a standard continuum mechanical setup as depicted in Fig. 3, allows us to track the location of each particle before and after the deformation with position vectors Xand x(t), respectively. Furthermore, the overall deformation of the sample can be captured, e.g. with the Green-Lagrange strain E= 0.5· Ft·F−1 (1) derived from the deformation gradient F=∂y(X, t) ∂X(2) with deformation map y(X, t). The MD code used only supports deforming the sample normal to the coordinate axes and thus a shear deformation cannot be applied directly. However, we can exploit the transformation rules of tensors and apply stretch in ˜ex-direction, compression in ˜ez-direction while leaving the ˜ey-direction unchanged as shown in Fig. 4 (a). A rotation of 45°around the˜ey≡ey-axis results in the desired pure shear deformation cf. Equ. 3 and Fig. 4 (b). 3 Characterization of polystyrene using MD PREPRINT ez≡EZ ex≡EXey≡EYF(X, t) X x(t)u(X, t)y(X, t) Ω0 Ωt t=t0t > t 0 Figure 3: Continuum mechanical setting: Initial ( Ω0) and current configuration ( Ωt) of the MD simulation box, position vectors Xandx(t), deformation map y(X, t), displacement u(X, t)and deformation gradient F(X, t), [10] h ˜Ei =1 2  λ2−1 0 0 0 0 0 0 0 − λ2−1 rot y 45 °⇐⇒ [E] =1 2 0 0 λ2−1 0 0 0 λ2−1 0 0  (3) In the MD simulation, this deformation is applied stepwise via adjusting the box lengths in the different directions accordingly which requires an NVT ensemble. In each of these loadsteps the Cauchy stress tensor is computed as ˜σ(t) =−[p(t)−patm1] (4) from the pressure tensor pand atmospheric pressure patm. Moreover, the stress is filtered loadstep-wise to exclude numerical noise and the equilibration at the beginning of each loadstep [ 10]. Rotating the stress tensor by 45°to theei coordinate system enables an investigation of the relation between the applied shear strain Exz(t)and the resulting stress state σ(t)which will be discussed in the following. (a) ˜ey˜ez ˜exLxlx=λ·Lxlz=λ·LzLz Ly=ly (b) ˜ey=eyez ex Figure 4: Pure shear deformation: (a) in ˜eicoordinate system with stretch λi, undeformed and deformed box dimensions Liandℓi, respectively and (b) in eicoordinate system obtained by 45°- rotation around ˜ey≡eyillustrating the actual shear deformation 4 Simulation results The components of the Cauchy stress tensor σdue to time periodic loading with strain amplitude of Ea xz= 4 % and maximum strain rate of ˙Emax xz= 1 % ns−1is shown over time in Fig. 5. The stress state is clearly dominated by σxz 4 Characterization of polystyrene using MD PREPRINT which follows the applied sinusoidal deformation while the other two shear stresses equal zero. On the contrary, the normal stress components σxx,σyy,σzzcoincide and fluctuate around zero with twice the deformation frequency. Since the other stress components are considerably small compared to the stress in xz−direction, only σxzis evaluated in the following. (a) (b) 0 5 10 15 20 25−20−1001020 timet innsCauchy stress σii inMPaσxx σyy σzz 0 5 10 15 20 25−20−1001020 timet innsCauchy stress σij inMPaσxy σxz σyz Figure 5: Comparison of the components of the Cauchy stress tensor σover time tfor time periodic pure shear test in Exz, strain amplitude Ea xz= 4 % and maximum strain rate ˙Emax xz= 1 % ns−1: (a) normal stress components σxx,σyy, σzzand (b) shear stress components σxy,σxz,σyz The shear stress σxzobtained from time-proportional simulations with different strain rates is depicted in Fig. 6. Overall, the faster the load is applied the stiffer the material responds indicating a clear dependence on the load rate which is characteristic for viscous material behavior. For each individual curve, the stress increases linearly up to 1 % strain and then flattens noticeably. Once the maximum stress is reached, σseems to remain constant. 0 1 2 3 4 5 6 7 80102030˙Exz Green-Lagrange strain Exz in%Cauchy-stress σxz inMPa˙Exz=20%ns−1 ˙Exz=10%ns−1 ˙Exz=5%ns−1 ˙Exz=2%ns−1 ˙Exz=1%ns−1 ˙Exz=0.5%ns−1 ˙Exz=0.25%ns−1 ˙Exz=0.1%ns−1 Figure 6: Influence of strain rate: Cauchy stress in shear direction σxzover Green-Lagrange strain in shear direction exzfor time proportional pure shear tests with maximum strain Emax xz= 8 % and strain rates from ˙Exz= 0.1 % ns−1 to20 % ns−1 Time periodic loading leads to a stress strain hysteresis as depicted in Fig. 7 for strain amplitude Ea xz= 8 % and initial strain rate ˙Emax xz= 1 % ns−1. In the first loading cycle, a stress maximum of 20.1 MPa and a minimum of −17.6 MPa are obtained while in the subsequently cyclic softening leads to a stabilization of the hysteresis reaching σmax xz= 16.2 MPa andσmin xz=−16.0 MPa . This behavior is qualitatively observed for the other strain amplitudes investigated, too. Since a stress strain hysteresis is obtained, inelastic effects, i.e. viscosity and/or plasticity, have to be present. These can be quantified by computing the dissipated energy density per cycle as dcyc 0=Z cycS: dE=Z cyc2Sxz·dExz (5) 5 Characterization of polystyrene using MD PREPRINT −8 −6 −4 −2 0 2 4 6 8−30−20−100102030 Ea xz=8% ˙Emax xz=1%ns−1Pr Green-Lagrange strain Exz in%Cauchy stress σxz inMPa1stcycle 2ndcycle 3rdcycle 4thcycle Figure 7: Time periodic pure shear: Cauchy stress σxzover Green-Lagrange strain Exzin shear direction for strain amplitude Ea xz= 8 % and maximum strain rate ˙Emax xz= 1 % ns−1, starting points Prfor relaxation test with the Piola-Kirchhoff stress S=JF−1·σ·F−tandR cycdenoting the integral over one deformation cycle. Fig. 8 (a) indicates that for strain amplitudes up to 2 %there is almost no dissipation and thus the samples behave purely elastic. For larger deformations, however, dcyc 0and therewith the inelasticity increase quadratically. Next to this dependence on the strain amplitude, there is also a correlation with the strain rate as depicted in Fig. 8 (b). For small ˙Exzthe largest dissipation occurs approaching an asymptote at 0 % ns−1. For increasing strain rates, however, dcyc 0decreases rapidly before flattening and reaching its minimal value. A likely explanation is that for small deformation rates, the polymer chains have more time to rearrange and thus to dissipate more energy compared to high strain rates. (a) (b) 0 2 4 6 8 100246 ∼x2 ˙Emax xz=1%ns−1 Ea xz in%dcyc 0inmJ mm−3 3rdhysteresis 4thhysteresis p olynomial t 0 5 10 15 200.811.21.4 Ea xz=4% ˙Emax xz in%ns−1dcyc 0inmJ mm−33rdhysteresis 4thhysteresis rational t Figure 8: Dissipated energy density dcyc 0derived from cyclic pure shear tests for (a) different strain amplitudes with Ea xzwith ˙Emax xzand (b) different maximum strain rates ˙Emax xzwithEa xz= 4 % 5 Relaxation tests So far, we identified the presence of inelasticity for large deformations, but we could not distinguish between viscous and plastic effects. In order to do so, we apply a zero strain to the MD samples after having been subjected to one deformation cycle. Thus the starting point of the relaxation tests coincides with the end point of the first stress-strain hysteresis, labeled Pr, cf. Fig.7. Since at this point the sample is not subjected to a deformation, the elastic part of the stress has to vanish and thus elastic effects are eliminated. On the contrary, the viscous part of the stress degrades over time due to relaxation effects such as rearrangement of polymer chains, while the plastic part of the stress remains unchanged. This evolution of the shear stress σxzover the relaxation time is depicted in Fig.9 for preceding strain amplitudes Ea xz from 2 %to8 %. First, we observe an initial increase of σxzwhich is more pronounced for larger preceding deformation. AtPrthe samples are experiencing zero strain, but maximum deformation rate ˙Emax xzwhich is now instantly reduced to 6 Characterization of polystyrene using MD PREPRINT zero. Since the present material exhibits partly viscous behavior, as displayed in Fig. 6, this discontinuity has to lead to a temporary increase of the stress since the deformation is completely prescribed. Afterwards, the stress degrades over time with decreasing slope eventually approaching a constant level. This relaxation effect is more distinct for larger Ea xz, which implies a correlation between viscosity and applied load. The remaining σxzindicates that the samples experience also plastic effects even for the smallest preceding deformation. 0 10 20 30 40 50051015 PrPrPrPr relaxation time t innsCauchy stress σxz inMPaEa xz=2% Ea xz=4% Ea xz=6% Ea xz=8% Figure 9: Identification of plasticity on the basis of relaxation tests with Exz= 0 % = const with preceding cyclic loading with different strain amplitudes Ea xz 6 Summary and outlook In this contribution, we thoroughly discuss the application of a pure shear deformation within the MD framework from a continuum mechanical point of view. We investigate atactic polystyrene as a model system and successfully employ a characterization scheme for the material behavior of polymers [ 10]. From time proportional tests, we observe a rate dependent material behavior and time periodic simulations reveal an elastic regime for small deformations. However, inelastic effects are present, which could be decomposed into viscous and plastic parts by relaxation tests. On the one hand, these findings underline the applicability of the used characterization scheme. On the other hand, the obtained data and the gained insights can be utilized to further validate the constitutive law calibrated for this model material in [ 11]. Embedding this constitutive law into multiscale techniques like the Capriccio method will enable us to investigate complex polymer phenomena such as fracture or the behavior of nanocomposites. Acknowledgments We thank for the very fruitful discussions with Michael C. Böhm from the Theoretical Physical Chemistry Group at Technical University Darmstadt and with Christian Wick from the PULS Group at Friedrich-Alexander Universität Erlangen-Nürnberg. Furthermore, we would like to thank Lukas Laubert for performing the MD simulations . Funding Sebastian Pfaller is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 396414850 (Individual Research Grant ’Identifikation von Interphaseneigenschaften in Nanokompositen’). Maximilian Ries, Paul Steinmann, and Sebastian Pfaller are funded by the DFG - 37747279 (Research Training Group GRK 2423 ’Fracture across Scales - FRASCAL’). 7 Characterization of polystyrene using MD PREPRINT