This contribution introduces an unconventional procedure to characterize spatial profiles of elastic and inelastic properties inside polymer interphases around nanoparticles. Interphases denote those regions in the polymer matrix whose mechanical properties are influenced by the filler surfaces and thus deviate from the bulk properties. They are of particular relevance in case of nano-sized filler particles with a comparatively large surface-to-volume ratio and hence can explain the frequent observation that the overall properties of polymer nanocomposites cannot be determined by classical mixing rules, which only consider the behavior of the individual constituents. Interphase characterization for nanocomposites poses hardly solvable challenges to the experimenter and is still an unsolved problem in many cases. Instead of real experiments, we perform pseudo exper- iments using our recently developed Capriccio method, which is an MD-FE domain-decomposition tool specifically designed for amorphous polymers. These pseudo-experimental data then serve as input for a typical inverse parameter identification. With this procedure, spatially varying mechanical properties inside the polymer are, for the first time, translated into intuitively understandable profiles of continuum mechanical parameters. As a model material, we employ silica-enforced polystyrene, for which our procedure reveals exponential saturation profiles for Young’s modulus and the yield stress inside the interphase, where the former takes about seven times the bulk value at the particle surface and the latter roughly triples. Interestingly, hardening coefficient and Poisson’s ratio of the polymer remain nearly constant inside the interphase. Besides gaining insight into the constitutive influence of filler particles, these unexpected and intriguing results also offer interesting explanatory options for the failure behavior of polymer nanocomposites. Keywords polymer nanocomposites; interphase layers; interphase properties; atomistic-continuum coupling Highlights: • Numerical study of mechanical interphases in polymer nanocomposites • Uniaxial deformation via atomistic-continuum coupling method A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT • Identification of material property gradients inside the interphase • Revealing elastoplastic behavior in the vicinity of filler particles 2 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 1 Introduction and outline The concept of interphases has been introduced based on experimental observations of the mechanical properties of adhesively bonded joints. These frequently show a strong dependence of the mechanical properties on the thickness of the adhesive layer. Thereby, interphase denotes those regions of the polymer in the substrate’s vicinity, in which the mechanical properties deviate from those of the undisturbed bulk material, i.e., far away from the substrate. These 5 property changes are caused by the interactions of the substrate surface with the polymer network and are typically highly dependent on the material pairing (epoxy on aluminum, polyurethane on steel, . . . ) [ 1,2,3,4], the manufacturing conditions (curing temperature, humidity, . . . .) [ 5,6,7], and finally the application environment (temperature cycles, aging, . . . ) [8, 9, 10]. In the case of composites, where filler particles are embedded into the polymer matrix, one can assume that similar 10 mechanisms lead to interphase-like shells surrounding each filler particle [11]. Regarding the particular case of nanocomposites, the concept of interphases may well explain the frequent observation that the overall mechanical properties cannot simply be deduced from the properties of the constituents (particles and matrix) by classical mixing rules [ 12,13]. Unlike microcomposites, for which these mixing rules have worked well, nanocomposites exhibit a considerably larger surface-to-volume ratio, i.e. the impact of the volume content of the 15 interphases surrounding the particles becomes dominant [ 14]. This size effect has been observed experimentally [ 15] and allows for the precise tailoring of the mechanical properties of the nanocomposites to the respective application, thus making nanocomposites an outstanding class of material. To exploit their full potential, physically reasonable structure simulations of nanocomposites are needed, which require detailed knowledge of the mechanical property gradients inside the interphases. 20 A sketch of the interphase situation in nanocomposites is given in Figure 1: Interphases can be described as spherical shells of finite thickness which surround the particles and show properties between those of the bulk and the particle materials. 121 2 NP IP bulk dmech. propertiesbulkNP property proﬁle? Figure 1: Schematic representation of interphases (IP) as forming between bulk polymer and nanoparticles (NP) and the corresponding mechanical property profiles. Three questions immediately arise in this regard: (i) How thick are interphases, (ii) which mechanical properties are affected at all, and (iii) according to which functional relation are the properties changing from the particle surface 25 towards the bulk polymer? Due to their particular location – one occasionally also speaks of buried layers – interphases usually pose a serious challenge to the experimenter, especially in the case of nanocomposites [16]. The outstanding properties and versatility of nanocomposites motivated many researchers to approach these questions from different directions. In the following, we present a brief overview of the methods and tools commonly used in this context. 30 3 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 1.1 Experimental studies Experimental studies using small-angle X-ray (SASX) scattering [ 17,18], and transmission electron microscopy (TEM) [ 19,20] were carried out to investigate the morphology of the interphases. Maghami et al. propose an experimental methodology to derive the interphase’s morphology, thickness, and permeability employing SAXS on polyethersulfone-silica composites [21]. 35 Broadband dielectric spectroscopy (BDS) helps to reveal the dynamic behavior of polymers and is used to identify the glass transition temperature of the neat polymer and polymer nanocomposites [22]. Another handy experimental method is Brillouin light scattering (BLS), which has been commonly used to examine interphases [ 23,24,25] and, in the case of nanocomposites, allows the measurement of mechanical quantities such as Poisson’s ratio and Young’s modulus in the submicron range [ 26,27,28]. Holt et al. [ 29] used BLS to determine the 40 bulk and shear moduli of P2VP-silica composites and compared the results for grafted and physically absorbed filler particles at comparable molecular weights. The non-grafted composite exhibits a more pronounced strengthening of the mechanical properties than the grafted one relative to the neat polymer. The authors attribute this strengthening to chain packing and stretching within the interfacial layer. Bulge tests of thin ( ≈100 nm ) polymeric films [ 30] allow for the measurement of stress-strain behavior of polymer 45 nanocomposites and thus identification of Young’s modulus, Poisson’s ratio, and yield strength [ 31,32]. However, there are no direct methods to locally measure the properties of the interphase. 1.2 Empirical models In order to predict the effective properties of polymeric composites, many empirical models, i.a. by Leidner et al. 50 [33], and Pukanszky [ 34], have been developed for micron-sized additives taking into account the influence of filler size, filler content, matrix-filler interactions, and interfacial area. These models were thoroughly extended to address the challenges of nano-sized filler particles and performed well compared to experimental results [ 35,36,37,38,39]. In contrast to classical composites, the interphase is a crucial factor in nanocomposites [ 40,41] and hence requires appropriate modeling, as in the three-phase model proposed by Ji et al. [ 42]. Using this model, Zare et al. [ 43,44] 55 subdivide the interphase into layers and conclude that a power law is best suited to predict Young’s modulus spatial profile for various polymeric nanocomposites. Due to van der Waals interactions, the nanoparticles tend to form clusters, which increases the effective filler size. As a result, this agglomeration has a considerable influence on the fillers’ effectiveness and thus on the material properties of the composite. Consequently, considerable effort is devoted to the adequate modeling of this effect [45, 46, 47]. 60 1.3 Numerical studies The outstanding properties of nanocomposites originate from complex mechanisms at the molecular scale. Since the aforementioned empirical models can only partially capture these, atomistic studies with molecular dynamics (MD) have been conducted to unravel the interfacial region. Shandiz et al. [ 48] evaluate the density distribution within the 65 interphase to generate homogeneous samples at various density levels. Subjecting these samples to uniaxial tension and identifying the respective Young’s modulus enables them to approximate the interphase’s stiffness profile. Odegard et al. [49, 50] introduced the equivalent-continuum modeling: A continuum mechanical constitutive law is calibrated by matching its energy under deformation to the energy obtained in preceding MD simulations under identical boundary conditions. 70 Due to MD’s spatial and temporal limitations, concurrent multiscale approaches have been increasingly developed, aiming to combine the accuracy of fine-scale methods with the efficiency of macroscopic descriptions. These approaches can be categorized into hierarchical and partitioned-domain methods [ 51]. In hierarchical methods, the coarse-scale provides the fine-scale boundary conditions and, conversely, the fine-scale provides the constitutive law for the coarse- scale. Here, both resolutions are applied simultaneously to the entire domain, resulting in the macroscopic material 75 response as an average of the fine-scale, thus also referred to as homogenization. A well-known example is the widely used FE2method, where a fine-scale representative volume element (RVE) is evaluated at every quadrature point of the coarse-scale finite element (FE) model [ 52,53]. Recent advances in computational homogenization for nanocomposites have been made, e.g., by Kumar et al., by comparing an interface-energetics-enhanced approach with a graded-interphase-enhanced formulation to capture the size-effect [54]. 80 Partitioned-domain multiscale approaches have been introduced for crystalline materials most prominently with the Quasicontinuum method [ 55,56], the Arlequin method [ 57,58] or the Bridging Domain method [ 59,60]. A thorough overview of these atomistic-to-continuum coupling strategies can be found in [ 51]. Based on the preceding methods, 4 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT Pfaller et al. introduced the C APRICCIO method [61, 62, 63, 64, 65, 66, 67] for the treatment of amorphous polymers, which we discuss in more detail in Section 2.1. 85 Another noteworthy atomistic-to-continuum coupling approach was introduced by Luan et. al [ 68]. They transfer local displacements between the two descriptions in an overlap region and are able to reproduce atomistic results without refining the continuum mesh to atomisitic resolution. A comprehensive overview of multiscale studies of polymer nanocomposites is given in [69, 70, 71]. 90 1.4 Outline In this contribution, we propose a rather unconventional procedure to determine the spatially varying mechanical properties inside the interphase in a continuum setting, which is a precondition to simulate systems with more than a few nanoparticles. To this end, we translate alterations of structural features at the atomistic or molecular level as caused by the polymer-particle interactions into continuum mechanical parameter profiles intuitively interpretable 95 and applicable to engineering problems. Figure 2 sketches the main steps of our study: Contrary to the frequently used empirical models, we rely on an MD-based model of the polymer composite, which integrates particle-based descriptions of the bulk polymer and the filler particles at the same scale (Figure 2 a and b). To perform physically sound deformation simulations of a polystyrene matrix containing silica nanoparticles, we apply the aforementioned Capriccio method (cf. Section 2 and Figure 2 c). In particular, we consider systems with two particles located at 100 several different initial distances, which turns out to provide a local strain measure well-suited to characterize polymer interphases. In other words, we use multiscale simulations of specifically designed nanocomposite systems to substitute real experiments. Due to the underlying MD description at the polymer-particle contact, these pseudo-experimental data can be assumed to naturally capture both the formation and the effects of the interphases. Based on these data, a finite element (FE) based inverse parameter identification procedure is applied afterward to identify and quantitatively 105 determine the gradients of the inelastic mechanical properties around the nanoparticles in the continuum setting (cf. Section 3 and Figure 2 d). The novelties of this contribution are •introduction of the interparticle strain as a suitable measure to detect interphases around nanoparticles and to characterize their mechanical properties, 110 • translation of MD structure into mechanical properties via a reverse engineering approach, • determination of the inelastic property profiles of the interphases, and •providing a methodology to calibrate continuum mechanical material models based on molecular pseudo- experiments. The idea to utilize concurrent multiscale simulations of two-particle systems as described above has already been 115 introduced in an earlier contribution of some of the authors [ 72]. However, by that time, the setup and equilibration of the nanocomposites’ initial MD description suffered from some minor unphysical constraints, which have meanwhile been resolved. Our previous observations are mainly reproduced by the new data presented in the following, except for systems with very close nanoparticles. 5 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT material descriptioninitial configurationsrepresentative samples a) Particle-Based material descriptionb) Sample preparation c) FEMD-pseudo- experiments d) Identification of IP properties FE coupled FE-MD MD Figure 2: Graphical outline: The behavior of the investigated polymer is prescribed via a particle-based material description, which enables the preparation of our samples. The resulting initial configurations are used in subsequent, coupled MD-FE pseudo-experiments. Representative samples are chosen to identify interphase properties in an FE framework. 2 Coupled MD-FE pseudo-experiments 120 2.1 The Capriccio method The Capriccio method is a recent multiscale technique, which is particularly designed for amorphous polymers. It has been developed in close collaboration between chemists and engineers [ 62,63] and has recently been thoroughly optimized to reduce artifacts arising from the coupling of very different material modeling approaches [ 65]. The Capriccio method is, according to [ 51], a partitioned-domain technique that concurrently couples a discrete domain 125 treated by MD and a continuum domain resolved with the Finite-Element-Method. At the outer boundary ( 2 nm ) of the particle domain, the MD equations of motion are supplemented by the friction and random terms of dissipative particle dynamics (DPD) [ 73,74]. On the one hand, the DPD terms act as a thermostat via the fluctuation-dissipation theorem. On the other hand, the random terms compensate for the missing interactions of the MD beads with particles outside the box and thus allow a reproduction of the dynamical behavior as observed 130 under periodic boundary conditions. In this way, coupling artifacts such as the reflection of deformation waves at the atomistic-continuum interface are suppressed [61]. A crucial aspect of the Capriccio method is introducing so-called anchor points, which are auxiliary particles responsible for the information transfer between the domains, cf. Figure 3 (left). FE nodes to APs MD beads to APsenergy blending bridging domain equilibrate FE with APs freeequilibrate ME with APs ﬁxed undeformed systemdeformed system A P reactionforces APpositions Figure 3: Capriccio method overview: Coupling of MD and FE domains through anchor points and elastic springs (left), information transfer between domains and staggered solution scheme (right). 6 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT From the MD point of view, these anchor points are spatially fixed and confine the MD particles to their designated 135 domain. To this end, harmonic potentials are employed to define an interaction between anchor points and MD particles. By doing so, classical periodic boundary conditions (PBC) can be substituted by so-called stochastic boundary conditions (SBC). This concept has been described in detail in [ 61] and enables deforming the particle domain via the anchor points. To this end, the anchor points are coupled to the continuum in a weak sense and may change their spatial positions according to the continuum deformation. This anchor point to continuum coupling bases on the Arlequin 140 method, introduced in the late 1990s by Ben Dhia [ 75], and realizes the coupling by blending the energies of the different domains in a bridging domain. A kinematic coupling constraint relates the continuum displacement field and the particle displacements in a weak sense. The Arlequin method, which was initially intended to couple FE domains of different resolutions, has later been extended to the coupling of continua with discrete, crystalline domains [ 57,76]. The anchor point density increases exponentially towards the atomistic-continuum interface to ensure a smooth transition 145 between the two domains [61]. The staggered algorithm used to solve the coupled problem is sketched in Figure 3 (right): Upon applying the first load step in the FE domain, the anchor points change their position according to the deformation of the continuum. Thus, they disturb the equilibrium in the particle domain, which has to be re-equilibrated subsequently. This, however, renders updated forces the MD particles exert on the auxiliary anchor points, such that the following FE step requires 150 adjusting the anchor point positions. After a certain number of MD-FE iteration steps, the system achieves equilibrium, and another load step can be applied. An appropriate coupling of MD and continuum domains requires that the material behavior of both descriptions fits together. To this end, a thorough characterization of the particle domain’s mechanical properties is required to enable choice and proper calibration of the constitutive law for the continuum. Based on MD pseudo-experiments as detailed 155 in [77,78,79], appropriate constitutive descriptions have been chosen, and an inverse parameter identification has been carried out [ 80]. A crucial aspect of this procedure is to carefully choose the load range in which the continuum model is valid: In general, larger deformations require more sophisticated models but enable to reduce the size of regions to be treated entirely on the particle scale. Accordingly, the numerical effort can be reduced. 2.2 Sample preparation and simulation setup 160 In this contribution, an atactic polystyrene matrix with two embedded nano-sized silica particles is subjected to uniaxial tension. To this end, well-established and computationally efficient coarse-grained force fields for both the polymer and the nanoparticles are used, which were derived via an Iterative Boltzmann-Inversion [ 81,82]. Initially, the two nanoparticles, each consisting of 873 beads resulting in a diameter of DNP= 4 nm , are placed at different initial distances LNPof 6 to 12 nm relative to each other. The surrounding 300 polystyrene molecules are generated with 200 165 beads each using a self-avoiding random-walk algorithm [ 81]. Subsequently, a preparation procedure described in more detail in [ 77] is applied, where the samples are first equilibrated at high temperature and then cooled down below glass transition temperature using the MD software, IBIsCO [83]. The MD box size is chosen large enough to provide sufficient space for the filler particles to arrange during the cooling down and thus not constrain them in an unphysical way. Consequently, the nanoparticles occupy rather arbitrary 170 locations and distances inside the matrix polymer. Note that the resulting samples are all approximately the same size. In contrast to the above steps, which are simulated under periodic boundary conditions, the coupled MD-FE scheme requires so-called stochastic boundary conditions [ 63,61]. Therefore, polymer chains protruding beyond the simulation box are cut and relocated inside. Additionally, anchor points are inserted in the bridging domain to confine the MD particles to the simulation box and enable the data transfer between FE and MD computations [64, 65]. 175 The samples are then embedded in a continuum cuboid of edge length LFE= 50 nm , referred to as the FE domain. To compensate for the nanoparticles’ arbitrary orientation, we rotate the MD system such that the connection vector of the two nanoparticles aligns with the global x-axis. Furthermore, a translation ensures that the center between the two particles coincides with the global coordinate origin resulting in a setup as shown in Figure 4. To counteract the imbalance caused by the transition to stochastic boundary conditions, the samples are equilibrated 180 once more, but now in the coupled MD-FE configuration and without any external loads. 7 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT z xy 50 nm50 nm 50 nmcontinuum bridging domain (width: 1 .6 nm) particle region (edge length: 17 .8 nm) nano particle (radius: 2 nm) Figure 4: MD-FE simulation setup: A particle region comprising two nanoparticles and a polystyrene matrix, as well as its surrounding bridging domain are integrated into a continuum cuboid in such a way that the nanoparticles are aligned on the x-axis and their center coincides with the coordinate origin. In total, 116 samples are prepared in this way and then uniaxially deformed in x-direction up to 6 %strain. The output quantities of interest are the overall tensile force FFEand deformed particle distance ℓNPas shown exemplarily in Figure 5. /lscriptFE /lscriptNP DNP FFEMDNPFE CGPS bead anchor point Figure 5: Deformed pseudo-experimental setup: Two spherical silica nanoparticles of diameter DNP= 4 nm inside a polystyrene matrix, enclosed by anchor points and surrounding FE mesh, with deformed edge length ℓFE, resulting external force FFEand interparticle distance ℓNP, depicted in ’ideal’ conformation (no rotation of MD-box required). 2.3 Results 185 Figure 6 visualizes the effects of the initial nanoparticle distance LNPon the resulting force FFE(right) as well as on theinterparticle strain εNP xx:=ℓNP− L NP LNP(1) depicted in Figure 6 (left), for the 116 samples investigated. The overall force-response appears to be linear elastic and is obviously independent of the initial particle distance, which is a reasonable result to be expected since particle volume ratio and layout of the specimens are the same for all specimens. Thus, FFEis not sensitive to LNPand hence 190 not suited for the characterization of interphases. 8 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT However, interparticle strain εNP xxreveals a clear dependence on initial particle distance LNP: Samples with initially close particles, i.e., with LNP−values between 5 to 6 nm, provide (blue) curves that are clustered around small values and show only a slight linear increase with the applied overall deformation. This means that the particles stick together and do not follow the surrounding polymer’s deformation, like if they are connected by a material much stiffer than 195 the bulk – or by stiff interphases. Contrary to this, systems with more distant particles show larger initial values for εNP xx, and the (yellow to red) curves increase significantly and nonlinearly with the applied overall deformation. In these cases, the soft bulk-like polymer is also located between the particles and increasingly deforms like the outer material. Consequently, interparticle strain εNP xxis well suited to detect and quantify differences in the deformation behavior that originates from the initial particle distance. It further indicates the ratio of bulk and interphase material located between 200 the particles. This sensitivity predestines this measure for the characterization of interphases. 6 8 10 12initial distance of NPs LNPinnm 0 1 2 3 4 5 60510 strain of FE domain εFE xxin %interp. strain εNP xxin % 0 1 2 3 4 5 6050100 strain of FE domain εFE xxin %ext. forces FFEin nN 0 1 2 3 4 5 6050100 strain of FE domain εFE xxin %ext. forces FFEin nN Figure 6: Influence of initial distance of NPs on uniaxial pseudo-experiments: Resulting external forces in deformation direction FFE(right) and interparticle strain εNP xx(left), versus strain of FE domain εFE xx, respectively. 3 Continuum-based identification of interphase property profiles To exploit the above-presented results of the Capriccio pseudo-experiments for the characterization of interphases, purely continuum mechanical, FE discretized surrogate models of the two-particle nanocomposite specimens have been developed. By applying a typical parameter optimization scheme, these FE models then enable a quantitative 205 determination of mechanical property profiles inside the interphases around the particles, based on the εNP xx(εFE xx)-data given in Figure 6 (left). The FE deformation simulations are thereby repeated with iterative adjustment of the material parameters inside the interphases until the interparticle strain profiles for two selected, representative configurations with initially close and distant nanoparticles are reproduced sufficiently well. 3.1 FE models and parameter optimization procedure 210 All of our FE models are fully three-dimensional, and since the filler particles are spheres and the simulation box is cubic, we may reduce the computational effort by exploiting symmetry, i.e., we consider only one octant of the hybrid MD-FE simulation box. Figure 7 depicts a close-up of the nanoparticle(s) (yellow, 4 nm diameter) embedded in the polymer matrix (grey). We choose four separate, equally thick interphase layers as a compromise between computational effort for and robustness of the inverse parameter identification and an appropriate discretisation of the 215 parameter profiles in the interphase. These interphase layers are denoted from IP1 (blue) next to the bulk polymer to IP4 (red) at the particle surface. The thickness of each of these layers is chosen as 0.5 nm , adding up to an overall interphase thickness of 2 nm , which is in line with results from Ndoro et al. [ 84], Riccardi et al. [ 85], and Liu et al. [ 86], who reported disturbed structural properties of the polymer for only the first two nanometers away from the nanoparticle surface. Finer discretizations of the interphase with more layers, i.e. more sampling points to approximate the actual 220 property profiles, are of course possible but imply correspondingly higher optimization efforts and longer computation times. For our purposes, four values are initially sufficient to capture the transition between the properties of bulk material and polymer at the particle surface. Bulk polymer, interphase layers, and particles are discretized with quadratic tetrahedral elements (Abaqus: C3D10) and such that conforming meshes result, i.e. there are no hanging nodes between the different domains, but the neighboring 225 9 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT Figure 7: FE surrogate model of the Capriccio simulations: Close-up to the nanoparticle (yellow) and the surrounding four interphase layers IP1 at the bulk polymer (grey) to IP4 at the particle surface. Symmetry is exploited with respect to the three front planes via appropriate boundary conditions. domains share all nodes at the interfaces. Consequently, all interfaces between bulk polymer, interphase layers, and nanoparticles are tied, i.e. no contact description is required, and there are no relative lateral displacements, friction, or separations, which correctly reproduces the situation of the Capriccio pseudo-experiments. The selected finite element size or mesh density naturally is a compromise between computation time and accuracy. All FE models used here consist of approximately 10 000 elements with about 17 000 nodes, see Figure 7. It has been 230 verified within the scope of a small convergency study that the relevant quantities FFEand interparticle strain εNP xx improve only in the low per mille range when even finer meshes are used. The comparably stiff silica nanoparticles are treated by small strain linear elasticity with Young’s modulus ENP= 162.1 GPa and Poisson’s ratio νNP= 0.272, which is in line with corresponding values available in the relevant literature. For the bulk polymer, the simplest hyperelastic large strain formulation, i.e. the NeoHookean material model 235 [87], is initially used. This constitutive law is natively implemented in Abaqus and considered sufficient for the purposes relevant here, although there are, of course, much more sophisticated models available for polystyrene. The material parameters of the bulk polymer have been obtained from MD simulation goal data using a least-squares minimization: Bulk polymer boxes have been deformed using pure MD, and the material parameters of FE surrogate models have been adjusted such that the stress-strain behavior of the FE simulation reproduces that of the MD simulations sufficiently 240 well (in a least-squares sense). The resulting elasticity constants of the bulk polymer E0= 842 .756 MPa andν0= 0.3 translate to c10= 162 .068 MPa andd1= 0.002848 , which are the corresponding material parameters required by Abaqus. The four interphase layers are treated by the same material model as the bulk polymer, but each have individual material parameters E1...4andν1...4, which are subject to optimization with respect to the Capriccio εNP xx(εFE xx)-data. In the following, two selected but representative cases for a system with initially very close and a system with initially 245 more distant nanoparticles are chosen out of the available Capriccio data presented in Subsection 2.3, Figure 6 (left). In the close case , the particle centers have an undeformed distance of LNP= 5.18 nm , i.e. the two interphases are strongly overlapping, compare Figure 8 (left), which shows the corresponding surrogate FE model. In the distant case , initial particle distance LNP= 7.95 nm and the two interphases touch each other only slightly, see Figure 8 (right). Both surrogate models cover an overall simulation domain of [50 nm]3and precisely reproduce the deformation histories 250 of the corresponding Capriccio simulations and provide the same output quantities, i.e. the overall reaction force FFE in the horizontal direction and the deformed interparticle distance ℓNP, which allows the computation of interparticle strain histories εNP xx(εFE xx)according to Equation (1). For the sake of simplicity, the latter are now normalized with respect to the applied overall strain εFE xx= [ℓFE− L FE]L−1 FEwithLFE= 50 nm , i.e. normalized interparticle strain histories 255 ˜εNP=εNP xx εFExx∈[0,1] (2) 10 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT are used to build the objective (or goal) function required by the optimizer to determine interphase property profiles. In the limiting case ˜εNP= 1, the two particles exactly follow the deformation of the surrounding polymer, like it would be the case when there is enough undisturbed bulk polymer in between them to compensate for possibly much stiffer interphases. Vice versa, in the limiting case ˜εNP= 0, the particles do not move at all with the surrounding polymer, as if they were glued together due to overlapping interphases with a much larger stiffness. 260 Figure 9 (left) depicts the resulting ˜εNP-profiles for the considered close and distant cases. While the former shows a roughly constant value of around 0.25, indicating particles that are mainly sticking together, the latter is more or less continuously growing from 0.6 to about 0.9, i.e. the, particles deform increasingly better with the surrounding polymer. To ensure the roughly linear elastic force response as measured by the Capriccio pseudo-experiments, the goal function additionally contains one of the two corresponding, almost identical force-strain curves as depicted in Figure 9 (right). 265 Both FE surrogate models (Figure 8) and the objective functions (Figure 9) are now used in an inverse parameter identification scheme to optimize the mechanical properties of the four interphase layers IP1...4. In an iteration loop, the two FE simulations are simultaneously carried out with Abaqus, the reaction forces and interparticle strains are subsequently extracted using a Python script and then transferred to an optimization routine available in Matlab. (lsqcurvefit , a gradient descent method optimizer with finite differences for the computation of the sensitivities) 270 iteratively updates the interphase layers’ material parameters until the deviation from the objective function falls below a prescribed threshold. The arithmetic mean ∆ =1 nnX i=1rh ˜εNP i,goal−˜εNP i,fiti2 (3) quantifies the deviation between the goal and the fit curves of the interparticle strain with number of data points n comprising the close and distant case. The results of four of these procedures with increasingly complex constitutive descriptions of the interphases are discussed in the following subsections. 275 3.2 Homogeneous elastic material without interphases The first result considers a homogeneous, elastic polymer without any interphases, i.e. we do not perform a parameter optimization but only simulate the two cases using the bulk parameters constant in all four interphase layers: E1...4=E0= 842 .756 MPa , ν 1...4=ν0= 0.3 (4) Figure 10 compares the resulting ˜εNP-profiles and force responses to the corresponding goal data. While the force curves are perfectly matching, see Figure 10 (right), the interparticle strain profiles are both missed by far, quantitatively 280 Figure 8: FE surrogate models for two representative cases of the Capriccio pseudo-experiments: 2D-cut of the close (left) and distant (right) case, cut through the center plane. Deformation and force are applied and measured in horizontal direction, respectively. 11 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 1 2 3 4 5 600.20.40.60.81 εFE xxin%˜εNPclose goal distant goal 0 1 2 3 4 5 6020406080100120 εFE xxin%FFEinnNclose force goal distant force goal Figure 9: Components of the objective function used by the optimizer to determine mechanical property profiles inside the interphases: Normalized interparticle strain histories ˜εNP(left) for the close and distant cases and corresponding overall reaction force curves (right). 2 4 600.20.40.60.81 ∆ = 0 .353 εFE xxin %˜εNPclose goal close sim distant goal distant sim 0 2 4 6020406080100120 εFE xxin %FFEin nN close goal close sim distant goal distant sim Figure 10: FE simulation results vs. Capriccio data for the close and distant cases, assuming homogeneous elastic material without interphases: Profiles of normalized interparticle strain (left) and overall force (right). as well as qualitatively ( ∆ = 0 .353). Both cases provide roughly straight lines, which are independent of the overall strain. One can obviously identify a purely geometrical effect of the different initial particle distances already for a homogeneous polymer matrix, which is reflected by the offset of about 15% in level between the two curves. Even without any interphases present, the initially closer particles tend to stick together more than their counterparts being farther away. The latter deform similarly to the outer matrix, i.e. ˜εNPis close to one, whereas the former reach only 285 about 80% of this value. Despite this purely geometric effect, both ˜εNP-curves are still far away from their goals, not only in terms of magnitude but in the distant case also in terms of shape. This implies the necessity to increase the structural and constitutive complexity of the models. 3.3 Elastic material with linear property profiles inside the interphases The polymer is still treated as hyperelastic, but we now allow different elasticity parameters inside the four interphase 290 layers. It turns out that a simultaneous, free optimization of all eight interphase parameters works poorly and is highly dependent on the selected starting values. For the sake of simplicity, we thus assume that Young’s modulus and Poisson’s ratio can change linearly from the bulk values towards the particle surface, whereby both an increase or a 12 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT decrease are allowed. The optimizer hence has to determine only the two slopes for Eandν, resulting in the following optima: 295 E4= 5.30197·E0= 4468 .27 MPa , ν 4= 1.6·ν0= 0.48. (5) The increase in νhas been artificially limited to a factor of 1.6 to avoid locking problems. An increase of Young’s modulus by a factor of 5.3 turns out to be the best option to capture both cases simultaneously as well as possible. A separate optimization for the close and distant cases would, of course, provide different values, which then fail even more in reproducing the complementary case. Figure 11 depicts the corresponding linear profiles of Eandνinside the four interphase layers, in each of which both values are piecewise constant. The increase in stiffness is an expected 300 result in view of the pseudo-experimental Capriccio data, which indicated that the particles are sticking together stronger the smaller their initial distance. A growth of Poisson’s ratio towards the surface can also be considered reasonable since lateral contractions of the polymer are increasingly hindered in the vicinity of the stiff particles. ...bulk 1 2 3 4 particle...012345162.1 interphase layerYoung’s modulus Ein GPa 00.10.20.30.40.5 Poisson’s ratio ν Figure 11: Optimal linear profiles for the elasticity parameters inside the interphases. Figure 12 compares the resulting ˜εNP-profiles and force curves to the corresponding goal data. The force responses are again perfectly reproduced, see Figure 12 (right). Linear profiles for the elasticity parameters obviously have no effect, 305 which can be considered reasonable due to the negligible volume ratio of ≈0.37 % occupied by the interphases. On the other hand, the normalized interparticle strain profiles are now much closer to their respective goals ( ∆ = 0 .093) and more separated from each other than before. Still, both ˜εNP-curves remain straight and independent of the overall strain. Furthermore, they are not able to simultaneously reproduce both goal data, at least on average. Marginal improvements can be achieved in this regard when more complex, nonlinear profiles for the elasticity parameters are considered for 310 the interphases, but the required dependence on the overall strain cannot be realized. Altogether, this implies that the current assumption of a purely elastic material is still too simple. 13 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 2 4 600.20.40.60.81 ∆ = 0 .093 εFE xxin %˜εNPclose goal close ﬁt distant goal distant ﬁt 0 2 4 6020406080100120 εFE xxin %FFEin nN close goal close ﬁt distant goal distant ﬁt Figure 12: FE simulation results vs. Capriccio data for the close and distant cases with optimized linear profiles of the elasticity parameters inside the interphases: Profiles of normalized interparticle strain (left) and overall force (right). 3.4 Elastoplastic material with linear property profiles inside the interphases The constitutive treatment of bulk and interphase polymer is now extended by a standard von Mises plasticity formulation with linear isotropic hardening as natively available in Abaqus. To this end, a yield stress σyneeds to be specified. For 315 the bulk polymer, a value of σy 0= 50.6 MPa has been identified from MD simulations of our model polystyrene under uniaxial tension [ 88]. A linear hardening can then be realized in Abaqus by simply prescribing the relation between plastic strains and stresses in the following tabular format, each for the bulk polymer and the four interphase layers, σpl εpl σy i 0 σy i+hi1 εpl0σ σy i 1hi Figure 13: Realization of plastic hardening in Abaqus: Relation between plastic strains εpland stresses σplfor each interface layer i= 1, ...,4. i.e. for i= 0,1, ...,4and whereby the hardening coefficients hidetermine the slopes of the plastic yield curves (cf. Figure 13). In fact, a hyperelastoplastic approach of this form is not fully adequate for polystyrene. A hyperelastic- 320 viscoplastic material law would be more realistic instead. However, due to the current lack of rate dependence in the Capriccio data, such a formulation cannot be calibrated appropriately yet and will therefore be treated in a subsequent contribution. The preliminary limitation to simple elastoplasticity nevertheless provides new and crucial insights to explain the mechanical interphase properties and thus understand nanocomposites’ behavior, as the following results demonstrate. 325 Like in the previous subsection, linear profiles for the elasticity and plasticity parameters throughout the interphases are now to be optimized. The yield stress of the bulk polymer is thereby kept fixed at the above mentioned value. Different to this, the hardening coefficient h0cannot yet be derived directly from bulk MD simulations, due to the limited strain range available. Therefore, the bulk hardening coefficient h0is left to optimization which provides h0= 0.3269 MPa , (6) i.e. a value only slightly off ideal plasticity. Furthermore, the optimizer has to determine the four slopes for E,ν,σy330 andh, which leads to the following increases relative to the bulk values, whereby also a reduction of the material 14 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT parameters would have been possible: E4= 8 .32406·E0= 7015 .15 MPa (7) ν4= 1 .56253·ν0= 0.469 (8) σy 4= 2 .28587·σy 0= 115 .67 MPa (9) h4= 21 .4272·h0= 7.01 MPa (10) Figure 14 compares the resulting ˜εNP-profiles and force curves to the goal data. The linear force responses are again perfectly reproduced (cf. Figure 14 (right)), although plastic deformations are now possible everywhere in the polymer, but obviously, arise only in a volume fraction small enough not to contribute significantly. The interparticle 335 strain profiles are reproduced very well now, including the required dependence on the overall strain resulting in ∆ = 0 .050. Still, some improvement is desirable, but considering plasticity in the material model is essential to capture the interparticle polymer’s deformation. Figures 15 and 16 depict the corresponding linear profiles of Eandνrespectively σyandhinside the four interphase layers. The increases in stiffness and Poisson’s ratio are comparable to the hyperelastic case, whereby the former is 340 considerably more pronounced now. Also, the optimized yield stresses and the hardening coefficients reveal significant increases towards the particle, i.e. the polymer is not only stiffer at the surface but also more resilient to inelastic deformations. This is an intriguing result concerning the frequently observed increase in material strength through the addition of nanofillers. A drastically increased yield point at the particle surface may very well explain the deceleration or even stop of crack propagation there. 345 Figure 17 contains contour plots of the equivalent plastic strain distributions around the particles and at maximum overall deformation for the close and distant cases. The volume fraction of plastically strained material is very small, which explains the apparently elastic overall responses in Figure 14. The region most strained with about a 23 % peak turns out to be the bulk polymer between the particles in the distant case, whereas only about one-half of this strain level is reached by the bulk outside the particles in the close case. The plastic zones naturally form in the tensile direction but 350 are shielded by the particle and the surrounding, significantly stiffer and stronger interphase layers, effectively impeding potential crack propagation. Finally, to provide better visualization of how the material properties change between the different interphase layers, one-dimensional stress-strain responses according to the different material parameter sets are depicted in Figure 18. Starting from the rather soft and nearly ideal elastoplastic curve of the bulk polymer at the bottom, stiffness and yield 355 stress clearly increase for the layers closer to the particle surface, whereas the hardening behavior remains more or less unaffected, at least in the eye-ball norm. 2 4 600.20.40.60.81 ∆ = 0 .050 εFE xxin %˜εNPclose goal close ﬁt distant goal distant ﬁt 0 2 4 6020406080100120 εFE xxin %FFEin nN close goal close ﬁt distant goal distant ﬁt Figure 14: FE simulation results vs. Capriccio data for the close and distant cases with optimized linear profiles of the elastoplastic parameters inside the interphases: Profiles of normalized interparticle strain (left) and overall force (right). 15 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT ...bulk 1 2 3 4 particle...01234567162.1 interphase layerYoung’s modulus Ein GPa 00.10.20.30.40.5 Poisson’s ratio ν Figure 15: Optimal linear profiles for the elasticity parameters inside the interphases. 3.5 Elastoplastic material with exponential saturation profiles For this final optimization, all previously linear property profiles, except that for Poisson’s ratio ν, are now replaced by exponential saturation functions, which is inspired by natural growth processes. Elasticity parameters and yield stress 360 of the bulk are identical to the previous values, while its hardening coefficient h0is again left to optimization, which provides h0= 3.53385 MPa , (11) i.e. a value about ten times higher than before, but still rather small compared to the yield stress. The material parameters Ei,σy iandhifor the different interphase layers i= 1, ...,4are simultaneously optimized according to an exponential saturation function of the format 365 xi=x0+αxx0 1 +e5−βxi, i= 1, ...,4, x =E, σy, h. (12) ...bulk 1 2 3 4 particle...0255075100125 ∞ interphase layeryield stress σyin MPa 0246810 hardening coefficient hin MPa Figure 16: Optimal linear profiles for the plasticity parameters inside the interphases. 16 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT Hence, the interphase parameters grow exponentially towards the particle, from the bulk value as the lower bound to a multiple αxof this as the upper saturation value, whereby the ’velocity’ of this transition is governed by another parameter βx. Together with h0and the slope of the linear increase of ν, eight variables in total are left to the optimizer, which then provides the following results E4= 6582 .56 MPa with αE= 6.81224 , βE= 3.35703 (13) ν4= 0 .998806 ·ν0= 0.299642 ≈ν0 (14) σy 4= 170 .060 MPa with ασy= 3.36544 , βσy= 1.46362 (15) h4= 3 .53849≈h0with αh= 0.06, βh= 0.3 (16) Figures 19 and 20 depict the corresponding profiles of Eandνrespectively σyandhin the four interphase layers. 370 Poisson’s ratio νand hardening coefficient hshow nearly constant profiles, as if the optimizer could not derive any meaningful benefit from modifying these values. Different to this, Young’s modulus Equickly takes its significantly larger saturation value, whereas, the yield stress σyincreases much slower and does not reach a true plateau value. All property increases are in ranges comparable to those of the previous linear case, and the much more sophisticated exponential saturation is obviously suited to separate the material parameters into important ( E, σy) and less relevant 375 (ν, h) for the different deformation characteristics of the polymer between close or distant particles. The one-dimensional elastoplastic stress-strain responses of the four interphase layers in Figure 21 provide another visualization of the property changes and emphasize the nonlinearity of the optimal parameter profiles. Figure 22 finally compares the resulting ˜εNP-profiles and force curves to the goal data. The linear force responses are still perfectly reproduced, and also the normalized interparticle strain profiles fit extremely well now with ∆ = 0 .035 380 which is only a tenth of the initial ∆ = 0 .353obtained without interphases (cf. Figure 10). Figure 17: Equivalent plastic strain distributions around the particles at maximum overall strain εFE= 6% for the close (left) and distant (right) situation. 17 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 0 0.04 0.08 0.12 0.16 0.2020406080100120 εσin MPa bulk IP-layer 1 IP-layer 2 IP-layer 3 IP-layer 4 Figure 18: One-dimensional stress-strain curves for bulk polymer and the four interphase layers as resulting from the optimization of elastoplastic material behavior with linear property gradients. ...bulk 1 2 3 4 particle...01234567162.1 interphase layerYoung’s modulus Ein GPa 00.10.20.30.40.5 Poisson’s ratio ν Figure 19: Optimal profiles for the elasticity parameters inside the interphases. 18 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT ...bulk 1 2 3 4 particle...0255075100125150175∞ interphase layeryield stress σyin MPa 01234567 hardening coefficient hin MPa Figure 20: Optimal exponential saturation profiles for the interphase plasticity parameters. 0 0.04 0.08 0.12 0.16 0.2050100150 εσin MPa bulk IP-layer 1 IP-layer 2 IP-layer 3 IP-layer 4 Figure 21: One-dimensional stress-strain curves for elastoplastic material with exponentially saturating profiles. 19 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 2 4 600.20.40.60.81 ∆ = 0 .035 εFE xxin %˜εNPclose goal close ﬁt distant goal distant ﬁt 0 2 4 6020406080100120 εFE xxin %FFEin nN close goal close ﬁt distant goal distant ﬁt Figure 22: FE simulation results vs. Capriccio data for the close and distant cases with optimized exponential saturation profiles of the elastoplastic parameters inside the interphases: Profiles of normalized interparticle strain (left) and overall force (right). 20 A methodology to characterize mechanical interphases ACCEPTED MANUSCRIPT 4 Conclusions This contribution introduces a method to numerically determine mechanical property profiles inside polymeric inter- phases in nanocomposites. To this end, we employ the CAPRICCIO method – a concurrent atomistic-to-continuum coupling technique – to perform MD-FE pseudo-experiments. We apply uniaxial tension to particularly designed 385 specimens consisting of a model polystyrene matrix surrounding two silica nanoparticles at various initial distances. The normalized interparticle strain ˜εNP, which quantifies how strong the particles are bound to each other by the polymer in between them, proves to be a well-suited measure to determine elastic and inelastic mechanical property profiles in the vicinity of the particles. To this end, a classical inverse parameter identification scheme based on pure FE surrogate models of the pseudo-experiments is carried out. Here, we select a close and distant particle setting as 390 suitable representatives for slightly and strongly overlapping interphases. The interphases themselves are modeled by 4 layers of equal thickness with individual material parameters. It turns out that a purely elastic material description is not sufficient to reproduce the pseudo-experimental data. Instead, an elastoplastic constitutive law is required where exponential saturation profiles for the material parameters deliver the best results. Interestingly, stiffness and yield stress are playing the dominant role, whereas Poisson’s ratio and the 395 hardening behavior are less important and remain constant within the interphase. The most intriguing result is that even an apparently linear elastic behavior requires elastoplasticity and gradients of the interphases’ plasticity parameters. This can provide entirely new options to explain the differences in the failure behavior of nano- and microcomposites. In this regard, our unconventional approach can be used as a toolbox to get insights into regions that are not accessible by experiments. 400 Extensions of the here presented method could investigate cross-linking polymers at flat surfaces to generate comparable results for adhesion applications. Furthermore, a homogenization of the interphase profiles determined with our method is indispensable to enable proper structure simulations of real nanocomposite volumes. Acknowledgments In memory of our dear friend Michael C. Böhm, whose continuous enthusiasm and support in numerous fruitful 405 discussions is sorely missed. 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 410 GRK 2423 ’Fracture across Scales - FRASCAL’).