Multimodel approach for accurate determination of industry-driven properties for Polymer Nanocomposite Materials

The need for effective and efﬁcient design and production of sophisticated materials with advanced performances on a competitive time scale is strongly driving the integration of material modeling and simulation techniques into material selection decision processes. Speciﬁcally, for complex struc- tural materials such as polymer-based nanocomposites (PNCs), there is a strong industrial demand for chemistry/physics-based models and modeling workﬂows able to predict relevant material properties in an accurate and reliable way. Under this perspective, in this work we describe the application of multiscale molecular modeling techniques for the choice of PNC materials for aerospace applications. The results are obtained in the framework of the European project Multi-scale Composite Material Selection Platform with a Seamless Integration of Materials Models and Multidisciplinary Design Framework (COMPOSELEC- TOR), funded by the European Commission within the H2020 call Advancing the integration of materials modeling in business processes to enhance effective industrial decision making and increase competitiveness . © 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The integration of modeling and simulations techniques to support material selection processes is one of the most compelling needs in advancing material industry and manufactory, due to the necessity of effective/efficient design and production of sophisticated materials, components and systems with advanced/extreme performance on a competitive time scale. In this arena and, specifically, for complex structural materials such as polymer-based nanocomposites (PNCs) there is a strong industrial demand for chemistry/physics-based models and modeling workflows able to predict relevant material properties (aka Key Performance Indicators or KPIs) in an accurate and reliable way [1,2].
Developments and improvements have been taking place in disparate communities considering different types of models and phenomena on several different length scales, with advancements in the so-called multiscale approaches, multi-disciplinary design optimization, and visualization [3][4][5][6]. Important links within a hierarchy of processing, nano/microstructure properties, and expected performance are currently available [7][8][9][10][11][12]. Nevertheless, they are far from being sufficient for materials design and selection, and suffer from a lack of integration across different types of models and related communities (especially discrete/continuum and modeling/experimental coupling and validation). Materials selection process (MSP) is fundamentally goal-oriented, aimed at identifying material structures and processing paths that deliver required KPIs. To be reliable, MSP must be built upon physical and engineering frameworks and based on methods that are systematic, effective and efficient in modeling complex, hierarchical materials such as PNCs. For material design and selection, understanding and quantifying the links between material structure at the nanoscale and their macroscopic effects is, therefore, essential. As manufacturing and testing is very time-consuming, expensive and sometimes unfeasible, the adoption of materials modeling and simulation approaches is therefore imperative. This implies the need for development and integration of models for the reliable prediction of material behavior at different scales as well as the derivation of efficient material-processing-property relationships.
For several reasons, current evaluation and selection methodologies for composites (and, thus, composites structural design) do not yet make extensive use of materials modeling tools, largely relying on empirical methods. On the one hand, since products have become more and more complex and sustainability issues coupled with a worldwide fierce economic competition are soaring, modeling and simulation -from design through material selection and production process -are more than ever needed. Their ultimate goals are shorter lead-time, better product quality, more competitive cost, and higher customer satisfaction. On the other hand, the complexity of actual technologies combined with the lack of a continuative and solid connection between industry and academia presents a condition in which a near-deadlock situation seems to emerge.
In this area, this work aims at demonstrating the use in COM-POSELECTOR of a multi-disciplinary, multi-model approach for the accurate, reliable, efficient and cost-effective industry-driven KPIs determination for PNC materials, thereby filling the gap between business processes and materials science/engineering workflows. Specifically, the results obtained for a case study pertaining to aerospace applications and the relevant KPIs determination (i.e., mechanical, viscoelastic and interfacial properties) will be presented and discussed. The high degree of modeling and data prediction integration into MSPs will then provide valuable information about current "holes in property space", enable what-if scenarios and provide mechanistic insights into the high performance of these fascinating materials.

The industrial case
While the use of the so-called "reinforced" plastics (RPs) has become a mainstay in aircraft production, the choice of polymeric matrix and, especially, (nano)fillers/fibers remains a largely unexplored scenario. Currently, many industrial aircraft manufactures are considering substituting thermoset RPs -most commonly employed in the air transport sector -with thermoplastic RPs. This last family of materials offers several key advantages over the former, including excellent fatigue and damage tolerance properties, shorter manufacturing cycles, lower moisture absorption, and almost complete recyclability. This last property makes thermoplastic composites in aerospace applications very attractive since any major aircraft industry and their suppliers produce hundreds of tons of waste material each year [13]. In the quest of new nano-RPs for aircraft specific applications, in COMPOSELECTOR two main material goals are pursued: a) the use of poly(ether ethyl ketone) (PEEK) as thermoplastic polymeric matrix in place of the widely employed thermosetting (e.g., epoxy) resins; b) the replacement of actual carbon fibers with multiwall carbon nanotubes (MWC-NTs). Accordingly, the following industrial KPIs determination were identified as tests cases for MSP implementation: i) glass transition temperature (T g ), ii) mechanical moduli, and iii) viscoelastic properties of the PEEK/MWCNT nanocomposite systems as function of the filler content. Furthermore, since the interfacial stress transfer between the CNT and the polymer matrix in the nanocomposite mainly depends on the joining strength among the nanofiller and the polymer, and the pullout of CNTs from the matrix is one of the main interfacial mechanisms experimentally observed, an additional KPI was selected, namely the estimation of the pulling forces of a CNT from the PEEK matrix. Accordingly, computational models on different scales were devised, run, and the predict values for the above KPIs where validated against experimental data, as described below.

General details
All simulations were carried out using the parallel version of the open source code LAMMPS (http://lammps.sandia.gov/), running on our own CPU/GPU cluster, using the COMPASS force field [14]. Starting from the optimized PEEK monomer unit, 3 independent configurations of a PEEK chain were built using the Rotational Isomeric State (RIS) method [15] at each desired temperature. For each of these chains, a three-dimensional (3D) simulation box with periodic boundary conditions (PBCs) was built, either using only PEEK chains or 1 multi-wall carbon nanotube surrounded by a number of PEEK chains required to achieve a specific weight fraction (vide infra), and the resulting cell was relaxed by energy minimization. Accordingly, 3 different 3D simulation boxes were obtained for each RIS generated PEEK chain, and the computational recipe for the determination of a given property was applied to each of these boxes (see Fig. 1).
Thus, each resulting property is expressed as the average value over 3 distinct simulations. The Berendsen thermostat/barostat [16] with coupling time constants T = 0.1 s and P = 0.5 ps, was used for temperature and pressure control during MD runs. Short-range intermolecular interactions of full atomistic polymer models were accounted for by the Lennard-Jones 9-6 potential function with geometric mean combination rules. Long-range electrostatic interactions of polymer models were handled with the particle mesh Ewald (PME) [17] method (accuracy 0.001 kcal/mol). The cut-off distance for truncation of all non-bonding interactions was set equal to half of the simulation cell length. In all MD simulations, integration of Newton's equations of motion was achieved using the velocity Verlet algorithm with an integration time step of 1 fs, if not otherwise states. The length of each MD simulation and the number of frames collected along each MD trajectory varied depending on the specific system and the relevant property to be estimated. Therefore, details will be given in the specific sections below.

Computational recipe for the prediction of glass transition temperature in PEEK/MWCNT nanocomposites
The estimation of the T g for the pure PEEK matrix and its MWCNT nanocomposites requires the knowledge of the corresponding system density values as function of the temperature. To the purpose, the following molecular dynamics (MD)-based procedure was adopted. First, for each PEEK-based system 3 3D simulation cells were built at 298 K as described in Section 2.2.1. MWNCT different weight fractions were realized by varying the number of polymer chains with respect to the fixed nanotube content. The resulting cell configurations were optimized via the following refinement protocol: i) initial system energy relaxation ii) brief MD simulation run (0.5 ns) in the constant volume-constant temperature (i.e., the canonical or NVT) ensemble at 298 K, and iii) final energy minimization via combined conjugate gradient/Newton-Raphson methods using 10 −3 kcal/(mol Å) as the convergence criterion for the energy gradient. The relaxed simulation boxes were then subjected to MD runs in the isochoric-isothermal (NPT) ensemble for 2.5 ns to reach realistic density values at the desired temperature and pressure (1 atm). The simulation cell was further relaxed by carrying out an NPT-based MD simulation annealing protocol as follows: each system was first heated from the desired temperature T (e.g., 298 K) to the temperature T + 200 K (e.g., 498 K) in increments of 20 K, and then cooled back to the initial temperature with the same temperature variation scheme. The total duration of the annealing process was 0.2 ns. At the end of the annealing process, each cell was again energy minimized, according to the following protocol: first, a combination of steepest descent/conjugate gradient methods was applied until the energy derivative on any atom decreased to 1000 kcal/mol Å. Next, Newton-Raphson minimization was performed until the energy derivative became less than 0.1 kcal/mol Å. For the final equilibrium density determination, starting from the well-relaxed simulation cells obtained at the end of the simulated annealing process further 1 ns NPT MD simulation run was performed for data acquisition at the desired temperature. The first 0.5 ns were discarded and only the last 0.5 ns of the equilibrated NPT MD trajectory were retained for data analysis. MD frames were saved every 5 ps, for a total number of frames of 100. The entire sequence described above, from PEEK chain generation via RIS to data production NPT runs, was repeated to estimate the density values for all PEK systems in the temperature range of interest (i.e., 298-500 K). From the equilibrium density values, the specific volume (reciprocal of density, V sp ) at each temperature was next determined, and the corresponding T g was estimated from the intersection of the linear regressions of the V sp vs. T curves in the two (glassy and rubbery) regimes.

Computational recipe for the prediction of mechanical moduli of PEEK/MWCNT nanocomposites
To determine stress-strain relationship for the different PEEK-MWCNT PNC systems, starting from the equilibrated systems obtained as described in Section 2.2.2, the first step requires achieving the minimum initial stress state for the corresponding periodic box. To the purpose, each simulation box was consecutively subjected to MD simulations in the canonical (NVT) and microcanonical (constant volume-constant energy, NVE) ensembles for 5 ns and 2 ns, respectively. Next, each unit cell size was adjusted to minimize the initial stresses using NPT MD simulations (simulation time up to several thousand steps, depending on the molecular model and cell size), followed by further equilibration (2 ns) in the NVE ensemble.
For each equilibrated system, a uniform strain field (0.5%) along the required directions was applied by proportionately scaling the corresponding unit cell and the atomic positions of the PEEK chains, while keeping the positions of the MWCNT unchanged, in order to preserve an un-deformed shape of the nanotubes. By applying this computational technique, each simulation cell was strained to calculated unidirectional tension/compression and hydrostatic tension/compression, respectively. In the first case, strain was applied along one direction a time, and the procedure was repeated for the two other directions starting from the unstrained cell. In the case of hydrostatic tension/compression simulations, the same strain field was applied simultaneously to all normal directions.
At the atomic level, stress is expressed in a virial form: (1) in which V is the volume of the simulation cell, given by the sum of each atom volume V a . v a i and v a j are the i-th and j-th components of the velocity of atom a, F ab i is the i-th component of the force between atoms a and b, and r ab j is the distance between the same two atoms. The negative sign is to introduced to express tensile stress as a positive quantity.
In the framework of continuum mechanics, the assumptions that i) both the polymeric matrix and the relevant nanocomposites are endowed with material symmetry and ii) the relation between stress and strain is linear both hold. Accordingly, the generalized constitutive relation of the equivalent continuum reads: in which ij and ij are the stress and strain components, respectively, and C ij are the elastic constants.
As mentioned above and in line with other works on similar systems [18][19][20], to predict stress-strain behavior under unidirectional tension or compression, a 0.5% strain was applied to the un-deformed cell along one direction a time. Under this condition, and taking direction 1 as the first choice, the only non-zero strain component is 11 . Expanding the first row of Eq. (2), we get: Utterly similar expressions can be written by deforming the cell along direction 2 or 3. Hence, for each simulated tension or compression three different values of the elastic constant C 11 are obtained. Since the strain is applied uniformly in each direction ( 11 = 22 = 33 ), C 11 can be further averaged.
Under conditions of hydrostatic tension and compression, the elastic bulk modulus K is given by Eq. (4): Accordingly, by recalling that, in this case, ii = 0.05 in all directions, K can be easily evaluated by means of Eq. (4) by inserting the normal components of stress as derived from the corresponding simulations.
Finally, the fundamental relationships that link the bulk modulus K and the elastic constant C 11 to the Poisson's ratio and the elastic (Young) modulus E,i.e., allow the determination of the value of and E once K and C 11 are known as a result of the simulation procedures reported above.

2.2.4.
Computational recipe for the prediction of the viscoelastic behavior of PEEK/MWCNT nanocomposites 2.2.4.1. Dissipative particle dynamics (DPD) theory in brief. To provide insights into the viscoelastic properties of PEEK/MWCNT nanocomposite systems, we resorted to Dissipative Particle Dynamics (DPD) simulations [21] and computed the relevant frequency-dependent storage (G') and loss moduli (G"). Briefly, in a DPD simulation the actual material is modelled as a collection of spherical particles (called beads), each representing a group of small molecules or extensive molecular fragments, which interact by pairwise additive forces expressed by a conservative, dissipative, and random potential. The overall force acting on a bead i can be expressed as and is calculated by summation over all other particles within a certain cut-off radius, r c , which represents the intrinsic length scale of the DPD model. The conservative force represents the excluded volume interactions between particles i and j in the dimensionless form F C ij = a ij 1 − r ij r ij , where r ij = r i − r j , r ij = |r ij |,r ij ij = r ij / r ij , and a ij is the maximum repulsion between particles i and j. The dissipative and random (F D ij = ω r ij r ij /( t) −1/2 ) forces act as heat sink and source, respectively, and the combined effect of the two forces performs as a thermostat, where ␥ is a friction coefficient related to the thermal noise amplitude via the fluctuation-dissipation theorem, 2 = 2␥k B T . (r) is a weight function, is a normally distributed random variable with zero mean and unit variance that is uncorrelated for different particle pairs, t is the time step of an integration scheme, and v i = v i − v j is the relative velocity of the i th and the j th particles. The equations of particle motion, dr i /dt = v i and dv i /dt = F i , are solved using as integration scheme the velocity-Verlet algorithm. When modeling polymer chains or complex molecules, typically two additional forces are acting between bonded beads: a harmonic spring connecting two adjacent particles i and j: F B ij = k b r ij − r 0 r ij , where k b is a spring constant and r 0 the equilibrium distance between the particles, and F A ijz = 1/2k Â sin Â 0 − Â 0 , where k is a spring constant and 0 the equilibrium angle between adjacent beads triples ijz in a row.

2.2.4.2.
Coarse-grained molecular models of PEEK and MWCNTs. The first necessary step of a coarse-grain (CG) calculation is the construction of the required CG models. This is usually performed taking into account the relevant structural/thermodynamic properties of the molecules under investigation and matching these with the corresponding atomistic features in order to assure the ability of the CG simulation to reproduce them. Accordingly, in this work PEEK and MWCNT CG models were derived by a direct comparison of the appropriate atomistic and DPD pair-pair correlation functions, according to a procedure extensively reported in our previous work [10,[22][23][24][25][26][27][28][29][30][31]. The next, important issue in DPD simulations is to capture the essential intra-and inter-molecular interactions taking place among all molecular actors of the mesoscopic simulations as expressed by the values of the conservative parameter a ij . This quantity accounts for the underlying chemistry of the system considered. In this work, we employed a well-validated strategy that correlates the interaction energies estimated from a lower scale (atomistic MD) simulations to the mesoscale a ij parameter values. Following this computational recipe, the interaction energies among the PEEK/MWCNT system components estimated from MD simulations were rescaled onto the corresponding mesoscale segments adapting the procedure described in our work [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36].

DPD simulation input parameters.
To study the viscoelastic properties of PEEK/MWCNT nanocomposite systems by DPD, the values of the simulation parameters were assigned as follows. The bond length between each DPD bead was set to 0.5 DPD units, while a bond harmonic constant of 100 k b T (with k b T = 1 in DPD) was attributed both to the CNT and the polymer. To incorporate an appropriate degree of bending rigidity in the CNT, both bond and angle constants were set equal to 100 k b T. This angle constant should allow to correctly account for bending stiffness and Young Modulus of real CNT/Polymer-based system. The value of the friction coefficient ␥ was set to 6.75 [37]. To simulate realistic systems with different content of MWCNTs by weight, each MWCNT and PEEK chain were modelled as 20 and 60 DPD beads, and the number of chains and MWCNTs in each system was varied, as shown in Table 1.
The repulsion parameters for equal beads (a ii ) was set to 25; the corresponding the interaction parameter between unequal beads, a ij , obtained from the procedure described above, was found to be equal to 34.
The cut-off radius r c is related to the volume of a DPD bead by the relationship r c = ( V b ) 1/3 , in which DPD number density = 3 and V b = 614 Å 3 is the volume of a PEEK monomer as used in atomistic MD simulations. This yields r c = 1.23 nm. As the a DPD simulation cell length of 44 DPD units, the corresponding real space dimension of the box side was ∼54 nm. Finally, the relation between a DPD time unit and a real time unit t is given by the relationship ∼25.7N 5/3 m [38,39], in which N m is the number of water molecules (V H2O = 29.7 Å 3 ) filling up a DPD bead. Accordingly, in our system N m ∼ 21, so that ∼ 4108 ps. The simulation time step adopted was 0.02, that is 82.2 ps, and each DPD simulation was run for 20 × 10 6 steps. Accordingly, each total simulation run corresponded to 1.640 ms.

Calculation of the viscoelastic properties.
The calculation of the frequency-dependent storage modulus G'() and loss modulus G" () was obtained converting the time-dependent shear modulus G(t) into its ()-dependent form [40]. In turn, G(t) was computed from the stress autocorrelation function given by Eq. (7): in which V is the volume of the simulation box, the angle brackets represent ensemble average, and ␣ˇ( with ␣ / = ␤) is the ␣-␤ component of the stress tensor, which is given by the virial expression: where m i , m i␣ , and v j␤ are the mass and the ␣and ␤-velocity components of bead i. The r ij␣ and F ij␤ in Eq. (8) represent the ␣component separation distance and the ␤-component force acting between beads i and j, respectively. To obtain more accurate results, i) smoother estimates of the stress autocorrelation functions were obtained by averaging the individual curves from the three stresses, given the equivalency of the off-diagonal elements of the stress tensor and ii) stress values were acquired every simulation time step.

Fig. 2.
Interface models for determining F coh (left) and Fpo (right) in separating a MWCNT from the PEEK matrix in normal and shear directions, respectively. Color legend as in Fig. 1.
Accordingly, G'(w) and G"(w) are given by: Once G' and G" are known from Eqs. (10) and (11), the corresponding dissipative factor tan␦ can be simply obtained as the ratio G"/G'. In order to obtain reliable estimates of the rheological properties of the PEEK/MWCNT nanocomposite systems, three independent simulation runs were conducted, each starting from independent configurations, and the values of G' and G" reported as average values.

Computational recipe for the prediction of the interfacial behavior in PEEK/MWCNT nanocomposites
Two different interface models were prepared to investigate the load transfer between a MWCNT and the PEEK chains, as shown in Fig. 2. Accordingly, both normal and shear separation forces can be calculated upon slowly displacing the MWCNT in normal and shear directions, respectively. Indeed, once the nanofiller is moved away from the polymer matrix, it experiences either normal (or cohesive) forces F coh , or shear (aka pullout) forces F po , depending on the motion direction of the MWCNT. In the first case, the polymeric chains are placed on top of the MWCNT ( Fig. 2A) whereas in shear separation the nanofiller is embedded within the polymeric matrix Fig. 2B). Accordingly, the nanotube is pulled in the z direction for normal separation whereas movement along the x direction is realize shear separation. After each molecular system building, 1 ns equilibration of MD simulations in the NVT ensemble was performed. After each system equilibration was achieved, in order to realize the two different MWCNT movement, an average displacement rate of 10 −4 -10 −5 Å fs −1 for normal separation and 10 −3 Å fs −1 for shear separation, respectively, were adopted. Fig. 3A shows the behavior of the specific volume as a function of temperature T for the pure PEEK and the corresponding PNC with 6% wt of MWCNTs as an example. The procedure to determine T g is initially based on visual data inspection, to identify the kink point in the V sp vs. T curves, followed by the curve division into two parts, corresponding to the regions T < T g (glassy regime) and T > T g (rubbery state). Next, each data set is fitted to a straight line by linear regression, and Tg is determined from the intersections of the two lines (Fig. 3A). The value of T g for all PEEK/MWCNTs systems predicted with this computational procedure, along with the relevant experimental data are compared in Fig. 3B. As can be seen in Fig. 3B, the predicted glass transition temperature of the thermoplastic PEEK (T g,sim (PEEK) = 414 ± 18 K, T g,exp (PEEK) = 416 K) is unaffected by the MWCNT dispersion up to nanofiller weight fraction of 15%, in agreement with experimental observations [41][42][43]. These data highlight a first, important difference between the currently employed thermosetting polymeric matrices and the thermoplastic PEEK. In fact, in thermoset (e.g., epoxy) polymers T g is seen to increase with increasing MWCNT concentration [44,45]. This opposite behavior is mainly ascribable to the formation of partial covalent bonding phenomena between the crosslinking polymer and the functional groups populating the surface of the carbon nanotubes. Ultimately, this evidence supports the notion that little (if any) chemical interactions develop between the thermoplastic PEEK matrix and the dispersed MWCNTs, their mutual interface being dominated almost exclusively by van der Waals forces.

Prediction of stress-strain behavior for PEEK/MWCNTs nanocomposite systems
The linear elastic stress-strain relations for the pure PEEK matrix are reported in Table 2.
The results in Table 2 allow for some considerations. First of all, the simulated values of the bulk modulus K nicely fall in the range of available experimental data for pure PEEK material [46,47]. Next, both K and C 11 values derived from compression methods are slightly higher than those obtained with the alternative, tension tests. To further check the validity of the results presented in Table 2, the corresponding values of the Young modulus E and the Poisson's ratio have been calculated via Eqs. (5) and (6). Accordingly, using the C 11 and K values obtained from compression tests (7.509 GPa and 5.656 GPa, Table 2), the magnitudes of E and are 3.869 GPa and 0.386, respectively; concomitantly, those obtained using C 11 and K values derived from tension simulations (7.227 GPa and 5.428 GPa, Table 2) are equal to 3.745 GPa and 0.385 for E and , respectively. Both sets of E and values are again in excellent agreement with the corresponding literature range values of 3.76 − 3.95 GPa for E and 0.3779-0.3931 for [46,47]. Overall, these results confirm that a nice correlation exists between values obtained from unidirectional and hydrostatic simulations, and that both computational procedures yield mechanical moduli in line with experimental data. Additionally, despite an applied strain of 0.5%, which is relatively high if compared with strain values adopted in real experiments, the behavior of PEEK still falls within the linear elastic regime.
Considering the same PEEK-MWCNT composite systems for which the glass transition temperature has been determined in Section 3.1, compared with pure PEEK a progressive addition of MWCNTs results in a corresponding increase of the elastic properties in both type (i.e., unidirectional and hydrostatic) of simulations. Specifically, the average value of the elastic constant C 11 (C 11 = 7.368 GPa, as obtained by averaging the two C 11 values in Table 2) increases from that of the pristine PEEK matrix by approximately 0.4, 14, 36, 62, and 146% with the addition of 3, 6, 9, 12, and 15% by weight of MWCNTs (see Table 2), while the average value of the bulk modulusK (=5.542 GPa, Table 2) grows approximately by 10, 24, 46, 64, and 132% with the same amount of added nanofiller, as seen from Table 3. Table 3 Values of the average elastic constant C11, average bulk modulusK, average Poisson modulus¯ and average Young modulusĒ for the pristine PEEK matrix and PEEK/MWCNTs nanocomposites as a function of the nanofiller concentration as derived from unidirectional compression/tension and hydrostatic compression/tension simulations, respectively. All simulations were run in triplicate. The last two columns of Table 3 list the values of the average Young modulusĒ and of the average Poisson ratio¯ , calculated inserting the corresponding values of C 11 andK in Eqs. (5) and (6). As we see,Ē also increases by increasing MWCNT concentration in the nanocomposite, reaching an improvement of nearly 136% at 15% wt with respect to the pure polymeric matrix. A comparison with available experimental data reveals an excellent agreement between simulated and measured values ofĒ [42], as shown in Fig. 4A.
Finally, the change of Poisson's ratio as a function of the MWCNT weight fraction is shown in Fig. 4B, from which we can notice that there is no specific trend for the relation between the Poisson's ratio and the concentration of the nanofillers.

Interface behavior in PEEK/MWCNTs nanocomposite systems
In polymer-based nanocomposites, interfaces play a determinant role in load transfer mechanisms between the polymeric matrix chains and the dispersed nanoparticles. In order to investigate these phenomena in PEEK/MWCNT based nanocomposites, both the normal and shear displacement at the interface between a MWCNT and PEEK chains were simulated upon strain application. As mentioned in Section 2.2.5, displacing the MWCNT in a normal direction yields information on the cohesive forces between the nanofiller and the polymer, while shear displacement mimics a real pullout experiment. By conducting molecular dynamics experiments under controlled displacement, the systems develops a reaction force at the polymer/MWCNT interface by virtue of underlying non-bonded interactions, which varies as the nanofiller drifts away from the polymeric matrix. Fig. 5A shows the behavior of the cohesive force F coh at the polymer/nanofiller interface as the MWCNT is displaced away from the PEEK chains along the normal direction.
As seen from the graph, the initial behavior of F coh is linear, then it reaches a peak and finally slowly declines. The value of the maximum in the F coh /displacement curve is estimated to be equal to 1131 pN at the corresponding displacement of 0.0019 nm. The corresponding profile of the pullout force F po vs. displacement is reported in Fig. 5B: the maximum value of F po locates at a displacement of 0.2, and is equal to 120 pN. Overall, these experiments support the idea that cohesive forces play a major role at determining the interface properties in PEEK/MWCNTs nanocomposite systems.   Table 2 (patterned bars), and the corresponding experimental measured data (filled bars) for pure PEEK and PEEK/MWCNT nanocomposites as a function of MWCNT concentration (% wt). Experimental data from [38]. (B) Average Poisson's ratio¯ values, calculated using data in Table 2 for pure PEEK and PEEK/MWCNT nanocomposites as a function of MWCNT concentration (% wt).

Viscoelastic properties of PEEK/MWCNTs nanocomposite systems
In silico dynamic frequency sweep tests were performed to investigate the viscoelastic properties of neat PEEK and its MWCNT nanocomposite systems. The results of these simulations are reported in Fig. 6.
As inferred from Fig. 6, the addition of the nanotubes to the polymeric matrix exert a stronger influence on the frequency dependence of the viscoelastic moduli at lower values. As the nanofiller concentration increases, the magnitude of both G'() ad G"() increases, while the gradient decreases. According to the theory of linear viscoelasticity, for homogenous polymer systems G'() and G"() should obey a scaling law behavior in the low w regime, characterized by slope values of 2 and 1, respectively [40]. Nonetheless, polymer matrices filled with nanoparticles often deviates from this ideal behavior. From Fig. 6 we see that the simulated data for the pristine PEEK matrix indeed follow the linear viscoelastic scaling laws. Yet, for MWCNT content of 3% wt onward this terminal behavior disappears, and both moduli tend toward a plateau-like regime, indicative of a transition from liquid-like to solid-like behavior. According to these data, a MWCNT concentration of 3% wt can constitute the so-called rheological percolation threshold, i.e., the concentration at which a continuous, three-dimensional network of nanotubes originates. Above this threshold, filler-filler interactions prevail over polymer-filler interactions, and a 3D filler network is fully established. These data are in agreement with the experimental evidence provided for PEEK/MWCNT nanocomposites in a concentration range encompassing the values considered in this work [48]. For such real systems, the behavior of G'() and G"() is in excellent qualitative agreement with the simulated data reported in Fig. 6; the percolation threshold in that case was found to be equal to 1% wt, again in good agreement with the one predicted by the present DPD simulations.

Outlook
In parallel, high performance materials design not only requires comprehensive material properties modeling but also understanding of risks, costs, and business opportunities for a range of decisions, from material selection to designing functional structural components and systems and process optimization. Last but not least, design and selection of nanocomposites must also accommodate societal requirements for health and sustainability (see Scheme 1).
The ultimate objective of the project COMPOSELECTOR, started January 2017 and funded by the European Commission within the Horizon 2020 (H2020) Program (under the Call NMBP-23-2016 − Advancing the integration of materials modeling in business processes to enhance effective industrial decision making and increase competitiveness), is to address all these aspects by developing and testing an integrated software platform − based on a multi-disciplinary, multi-model and multi-field approach − for decision-making in the selection, design and fabrication of PNC materials. One of the most innovative and ambitious features of COMPOSELECTOR is the establishment of a new paradigm in PNCs selection according to which the cornerstone of successful new, high-performance nanocomposite materials for advanced applications is the identification of all concurrent material-selecting factors along the entire value chain. In particular, the adoption of an approach that advocates the heavy use of computational materials science and physics rather than, e.g., property cross-plots which, albeit commonly used, lead to less competitive products, is the key factor to reach COMPOSELECTOR final goals.
The available set of composite materials and related manufacturing processes is growing in number and type. Therefore, for nanocomposite materials designers and manufacturers understanding and evaluating multiple materials and processing scenarios in the early stages of product selection/design are no longer an option. Information gained during the conceptual phase creates a product knowledge foundation for improved decisionmaking, thereby increasing the potential of innovation within new material and structure concepts, as sketched in Scheme 2.
In this scenario, the European Commission, within the H2020 Framework Programme, has dedicated substantial economical efforts to foster research projects aiming at advancing the integration of material molecular simulations in industrial business decision processes. In particular, the EU fully recognizes the need for a Business Decision Support System (BDSS) that supports the selection of the optimal material and process taking into account the implementation costs but also the associated risks, uncertainties and costs related to the modeling and simulation; a priority, especially for small medium enterprises (SMEs). Such a BDSS should enable the integration of materials modeling and business tools and databases into a single work-flow, allowing for flexibility of application to different industrial sectors. Specifically, it should be based on a framework that allows the flexible integration of existing or future discrete and continuum materials models with structured and unstructured data from multiple data bases containing materials data, commercial data and information on market trends, pricing, customer needs and demands. The BDSS should then i) enable a multi-criteria optimization over all stages of product development by allowing the end-user to mirror the operational framework of their company, ii) be structured in such a way to allow back-engineering from the end-goal, and iii) be designed such as to optimize the integration of humans in new more complex industrial environments. The results yielded by the BDSS should ultimately be validated against measurements, existing data and real financial arguments.

Conclusions
In a very broad sense, molecular modeling and simulation can be defined as the use of computational methods to describe material behavior at the atomistic or molecular level, contrarily to continuum-based modeling, in which atomic-level phenomena are neglected. In the engineering field, molecular simulations can be Scheme 2. Implementation view of the integrated COMPOSELECTOR business-decision support system (BDSS). The BDSS gathers information along the entire process to generate user informed decisions. used in two modes: i) the discovery mode, in which new phenomena are predicted that have not yet been observed experimentally, or explanations are sought for known phenomena that are not understood, and ii) the data mode, which substantially consist in the accurate prediction of materials properties relying on little or no experimental information. The former way of performing molecular simulations has benefitted several industrial sectors, including the aerospace and nuclear power industries. The alternative mode has an already established return on investment in industry. The perhaps most widely quoted example is that measuring the heat of formation for a single molecule cost $100,000 in 2000, about 50 times the cost of an accurate quantum chemical calculation [2]. Nonetheless, as mentioned in the introduction, the integration of modeling and simulations techniques to support material selection processes remains one of the most compelling needs in advancing material industry and manufactory, due to the necessity of effective/efficient design and production of sophisticated materials, components and systems with advanced/extreme performance on a competitive time scale. And this is particularly true in the case of the design and optimization of high performance polymer-based nanocomposites for advanced applications.
Accordingly, in this paper we positioned ourselves in the Material layer of the COMPOSELECTOR BDSS (Scheme 2), and developed a series of molecular simulation recipes for the fast and reliable prediction of some polymer nanocomposite key performance properties (KPIs) in the aerospace industry, namely the glass transition temperature, T g , the mechanical moduli (K, E, C 11 ), the interface forces F coh and F po , and the viscoelastic properties G'() and G"() as a function of the added nanofiller concentration. The results of such predictions, validated against available experimental data, will be fed both to a Material Data Base and, most importantly, to the business layer where they concur, together with all other relevant required information (e.g., process parameters, costs, sustainability, etc.), to generate user informed decisions.
Specifically, in this work we estimated these properties for a high-value polymer, PEEK, and the relevant MWCNT composites. Atomistic-based simulations were employed to predict the first three sets of KPIs, whereas a multiscale approach, involving both atomistic and mesoscale calculations, was adopted to investigate the viscoelastic behavior of these materials. The direct comparison of computer-based predictions with available experimental data, where available, has revealed a good to excellent agreement, thereby validating these in silico experiments and enabling them to be included as a tool to support business decision process in the development of further, new polymer-based nanocomposite materials.