Crystal nucleating agents offer an effective strategy for controlling the morphology, dimensional stability and rate of solidification of polymers during processing. Molecular dynamics (MD) simulation can shed light on nucleation behavior at the nanoscopic length and time scales over which nucleation occurs. In this work, crystal nucleation of a polyethylene oligomer, n-pentacontane, on three graphene-like substrates, hexagonal boron nitride (hBN), molybdenum disulfide (MoS2), and tungsten disulfide (WS2), was simulated, and the thermodynamic efficiencies of these substrates as nucleating agents were determined. Experimental measurements of heterogeneous nucleation of a high-density polyethylene on nanoparticles of these three graphene-like materials were performed using the method of dispersed microdroplets in an immiscible polystyrene matrix. Qualitative agreement between simulations and experiments was observed for trends in nucleation rate, J, and interfacial free energy difference, Δσ, with \(J_{\text{hBN}} > J_{\text{MoS}_{2}} > J_{\text{WS}_{2}}\). The simulations are then used to gain additional insight into the mechanisms of nucleation. Epitaxy is confirmed in all systems, with small mismatches in lattice spacing being accommodated by strain in the oligomer crystal. However, epitaxy alone is insufficient to explain the observed trends. The strength of interaction between the nucleating agent and the polyethylene oligomer is found to be the strongest predictor of nucleating agent efficiency. The strength of interaction is in turn related to the density of interaction sites at the interface: hBN has the highest density, and thus the fastest nucleation rate.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Introduction
In polymers, the kinetics of crystallization is largely responsible for the characteristic semi-crystalline morphology, which in turn imparts many favorable properties to these materials [1]. These kinetics are typically dependent upon the cooling rate and/or temperature at which the polymer is processed, which in turn affect the relative rates of crystal nucleation and crystal growth as well as the number and size of crystallites in a representative volume. Importantly, the nucleation rate and number density of crystal nuclei also depend on the presence of certain foreign impurities or proprietary additives that act as nucleating agents (NAs). Such nucleating agents facilitate the formation of stable polymer crystal nuclei, a process known as heterogeneous nucleation. By stabilizing the formation of small crystal nuclei, NAs reduce the kinetic barrier to nucleation and can increase the rate of nucleation by orders of magnitude compared to homogeneous nucleation [2]. What constitutes an effective NA, however, is still largely a matter of empiricism [3]. Continuum theories adapted from the nucleation of bubbles and droplets invoke contact angle phenomena based on differences in interfacial energies [4]. On the other hand, molecular theories typically emphasize the roles of epitaxy and grapho-epitaxy. Epitaxy involves the matching of crystal lattice dimensions between the NA and the polymer to maximize favorable interactions [5‐8], while in grapho-epitaxy, topological features like surface roughness or edges stabilize crystalline fragments of polymer chains through increased numbers of interactions [9‐11]. Recent molecular simulations have shown that not only the commensuration of interactions between polymer crystal and NA are important [12, 13], but also the surface compliance of the polymer crystal or NA can play a role, allowing the accommodation of mismatches in epitaxy through small crystallographic strain [12‐14].
Molecular simulations are an excellent tool for investigating polymer nucleation. Simulations allow for the direct observation (in silico) of events that occur on nanoscopic spatiotemporal scales. Epitaxy and crystallographic strain to improve epitaxy can be quantified and correlated to the values of force field parameters that are used to model different materials. For example, Bourque et al. examined heterogeneous nucleation of polyethylene oligomers on monatomic tetrahedrally coordinated (“diamond-like”) and hexagonal lattice (“graphene-like”) materials in precisely this way [12, 13]. Using different values for the 3-body bond angle, both types of materials were modelled using the Stillinger–Weber (SW) force field, originally developed for silicon [15]. The three key parameters in the force field are ε, λ, and σ, which control the lattice rigidity and atomic spacing of the simulated NA. Realistic behavior for a variety of materials can be reproduced by tuning these parameters [15‐18]. In the work of Bourque et al., these parameters were scanned over ranges of up to ± 40% relative to certain reference materials (diamond and graphene, respectively), to simulate 150 different materials, some of which have analogues in the real world but many of which are purely theoretical [12, 13]. From these parametric scans, useful heuristics were extracted and recommendations for improvements in the design of nucleating agents were made.
Anzeige
In particular, among the graphene-like, platelet materials, Bourque et al. predicted that hexagonal boron nitride (hBN) would be an even better NA for high-density polyethylene (HDPE) than graphene (GN) [13], which is itself well-known to be an excellent NA for HDPE [19‐22]. Often seen as an alternative to GN, hBN has excellent thermal and oxidative stability [23] and has been shown to improve the mechanical [24, 25], thermal [24‐27], gas barrier [28], and electrical properties [29] of HDPE+hBN nanocomposites over those of the base polymer. Similar improvements have been noted for other polymers as well [30‐34]. Only recently, however, has hBN received attention as a crystallization modifier [35, 36]. Ayoob et al. [29] observed a significant change in the morphology, crystallinity, and crystallization temperature of HDPE when crystallized in the presence of hBN, but the mechanism(s) behind such changes were not identified. Evidence for epitaxy has been experimentally observed with HDPE and reduced graphene oxide [19]. Because there is a close structural similarity between GN and hBN, epitaxy is also likely to be present in HDPE+hBN nanocomposites.
Hexagonal boron nitride has been studied with molecular simulations using a variety of force fields [37‐40], and also as a nanocomposite with HDPE [35, 41]. Reactive force fields are generally the most accurate, but also are computationally expensive. Relative to those force fields, the SW potential presents a computationally efficient model that makes long nucleation simulations feasible.
In addition to GN and hBN, two materials that are closely related to the hexagonal lattice but have non-planar bond angles are the transition metal dichalcogenides (TMDs), molybdenum disulfide (MoS2) and tungsten disulfide (WS2). In each monolayer of these TMDs, the transition metal atoms are sandwiched between two layers of sulfur atoms. These structures are contrasted with hBN in Fig. 1. Like hBN, TMDs are of interest as additives. Improvements in mechanical [42‐47], thermal [42, 44, 46, 48], and electrical properties [46] have been noted for nanocomposites containing MoS2 and WS2, with the potential for novel applications arising from polymer nanocomposites with self-healing [49], electrical storage capabilities [50], antibacterial properties [51], and biodegradability and biocompatibility[52]. Finally, SW parameterizations have already been developed for TMDs [53‐55] and have been used to model nanoplatelets undergoing mechanical processes [56‐58].
Figure 1
Structure of hexagonal boron nitride (a, b) and transition metal dichalcogenides (c, d). (a) and (c) are side views of a single monolayer, and (b) and (d) are top-down views. ahBN and aTMD are the unit cell parameters for the crystal lattice, and dZZ and dAC are the distances between parallel “zigzag” and “armchair” lattice structures, which are present in both types of materials.
×
Experimentally, the method of immiscible blending (IB) [59] has been found to be particularly powerful for the study of heterogeneous nucleation in polymers [22, 60‐63]. This method is a variation of the original droplet nucleation experiments of Vonnegut [64]. In the IB method, a crystallizable polymer is melt-blended as the minority component with a second, immiscible polymer, resulting in millions of micron-sized domains that serve as independent vessels for the study of nucleation phenomena. The domains are small enough that their complete crystallization is nearly instantaneous once nucleation occurs, but nucleation/crystallization within each domain is completely independent of that in any other domain. The result is an ensemble of nucleation events that release an exothermic heat flow that can be detected by differential scanning calorimetry (DSC). When an additive is pre-mixed with the crystallizable polymer prior to melt-blending, heterogeneous nucleation can occur. Heterogeneous nucleation of isotactic polypropylene (iPP) in the presence of various additives and melt-blended with polystyrene has been reported [60]. Heterogeneous nucleation has been reported for HDPE melt-blended with both iPP [61] and polystyrene [62], as well as when pre-mixed with certain pigments [63] and GN [22] as additives. This method is robust and can be extended easily to other additives, as well as other crystallizable polymers.
Anzeige
In this work, we study the heterogeneous nucleation of semicrystalline HDPE on hBN, MoS2, and WS2, both computationally by molecular simulation and experimentally by the method of immiscible blending and DSC. Computationally, the nucleation rate for each NA is measured, and the effects of epitaxy, lattice rigidity of the NA, and NA—polymer interaction strength are examined. Experimentally, nucleation rates are calculated over a range of crystallization temperatures and NA loading, and the thermodynamic nucleation efficiency for each material is quantified. Finally, a comparison between these two methods is made.
Methods
Materials
All materials were obtained from Millipore-Sigma. High-density polyethylene with a melt flow index of 12 g per 10 min (190 °C, 2.16 kg load) and an as-received density of 0.952 g mL−1 was used. The majority component for immiscible blending was atactic polystyrene (PS) with Mw = 192 kDa. Hexagonal boron nitride (hBN) nanoplatelets were obtained as a powder with a nominal particle diameter of 137 nm. MoS2 and WS2 nanoplatelets were also received in powder form, both with nominal particle diameters of 90 nm.
Nucleation rate
Molecular simulation
n-Pentacontane (C50) is an oligomer of polyethylene and was used as a model for HDPE in molecular simulations of heterogeneous nucleation. C50 was modeled as a chain of united atoms (UA), with each UA representing a CH2 or CH3 group. The force field developed by Paul, Yoon, and Smith (PYS) [65] and subsequently modified by Waheed et al. [66, 67] was used for C50. This force field comprises valence interactions and two-body Lennard–Jones interactions for nonbonded UAs. It has been shown to reproduce the melting temperature and enthalpy of fusion of n-alkanes [68], including the melting temperature of C50 [69]. This modified PYS force field was also used previously in studies of the heterogeneous nucleation of C50 on polyethylene [69], and diamond-like [12] and GN-like materials [13].
For NAs, the Stillinger–Weber (SW) force field was used to calculate interactions within each monolayer sheet [15]. The overall potential for the SW force field is given by
where i, j, and k index distinct atoms, rij is the distance between atoms i and j, θijk is the bond angle formed by atoms, i, j, and k, and ϕ2 and ϕ3 are the two- and three-body interactions, given by
The parameters that are most relevant to this study are the strength of the two-body interaction (εSW), the strength of the three-body interaction (λSW), and the atomic spacing (σSW). In order to compare to previous work with GN-like materials [13], a new parameterization for hBN was optimized empirically, keeping many of the parameters in the SW potential constant, and only varying εSW, λSW, and σSW. Bourque and Rutledge parameterized the SW force field for GN by optimizing parameter values in order to reproduce experimental elastic constants with molecular simulations [18]. A similar method was followed here for hBN. Details are provided in the Appendix. Due to the complex structure of TMDs relative to planar hBN, a new parameterization for TMDs was not made. Instead, SW parameterizations for single-layer TMDs were obtained from the literature [54]. Parameters for MoS2 were adjusted to reproduce experimental phonon dispersion curves [70]; subsequent calculations of the Young’s modulus were in reasonable agreement with experimental values [71, 72]. The SW force field was parametrized for WS2 to reproduce phonon dispersion curves obtained from ab initio calculations [73], and resulted in simulated values of the Young’s modulus that were in agreement with ab initio calculations [74, 75]. For both of these force fields, the pair_style sw/mod in LAMMPS was used to compute three-body contributions to the force field over proper angle types, as described by Jiang [53]. For this modification, cosθijk − cosθ0,ijk in Eq. (3) is defined as δ, and is scaled by a switching factor fc, given by
where δ1 = 0.25 and δ2 = 0.35 for TMDs. This modification reduces the number of atom types needed to uniquely define all bond angles in TMDs from twelve atom types per sheet to only two atom types per sheet [54].
Multi-layered 2D materials were constructed so that the NA thickness was larger than the cutoff distance of C50–C50 intermolecular interactions. The interaction between sheets was modeled using a Lennard–Jones (LJ) potential, as was done previously to model multi-layered sheets of GN [18], given by
Note that σLJ and εLJ are distinct from σSW and εSW. This potential was applied to both boron and nitrogen atoms in hBN and was optimized empirically as described in the Appendix. However, to avoid extensive modification of the force field for TMDs, the inter-sheet LJ potential was applied only between the metal atoms in TMDs, which is sufficient to maintain the experimentally observed inter-sheet spacing between the midpoints of each sheet. This was done for the sake of simplicity; additional potentials between sulfur atoms would have increased the number of parameters and were not necessary to maintain the correct distance. σLJ was set to reproduce the experimentally observed inter-sheet spacing, d = 21/6σLJ [76, 77]. The strength of the inter-sheet interactions, εLJ, is small compared to that controlling the stiffness of the 2D platelet, and small variations in the inter-sheet distance are expected to have little influence on the interface between C50 and the NA substrate. Therefore, approximate values were taken from GN [18] and applied to MoS2 and WS2, while the value for hBN was obtained by the force field fitting procedure described in the Appendix.
Intermolecular interactions between C50 and all substrate atoms were also modeled with LJ interactions, following Eq. (5). LJ parameter values for each substrate atom type were obtained from the Universal force field (UFF) [78], and Lorentz–Berthelot mixing rules were applied. The same cutoff distance was used for all LJ interactions between UAs in C50 and between UAs and substrate atoms. Tables 1, 2 and 3 contain the values of all relevant force field parameters.
Table 1
Force field parameters for non-bonded Lennard–Jones interaction potential
Two-atom force field parameters for Stillinger–Weber interaction potential
BN
MoS2a
WS2a
σSW [Å]
1.28
1.252
0.889
εSW [kcal mol−1]
68.63
23
23
A
5.8341
6.918
5.664
B
0.6022
7.223
24.525
p, q
4, 0
4, 0
4, 0
γ, aSW
1.2, 1.8
1, 2.523
1, 3.558
aReference [54], where bond angle interactions have been re-defined to reproduce TMD bond angles
Table 3
Three-atom force field parameters for Stillinger–Weber interaction potential
B–N–B and N–B–N
Mo–S–S(a)
S–Mo–Mo(a)
W–S–S(a)
S–W–W(a)
λSW
27.36
67.88
62.45
37.69
33.55
θ0 [°]
120
81.78
81.78
81.78
81.78
aReference [54], where bond angle interactions have been re-defined to reproduce TMD bond angles
Sheets of each NA were built using crystal files from the Springer Materials database [79‐81]. Lateral substrate dimensions were approximately 6.4 × 6.4 nm2, chosen to be longer than fully extended C50 chains so that chains would not impinge upon their periodic image when crystallized, but not so large that multiple nuclei could form and frustrate crystal growth, as previously noted in a few cases for GN [13]. The number of sheets for each NA material was chosen to be thicker than the cutoff distance for C50 interactions: 4 sheets for hBN, 2 sheets for both MoS2 and WS2. Three hundred crystalline C50 chains were placed on top of the NA substrates. Periodic boundary conditions were used in all three directions, so that the C50 domain was bounded on either side in the z-direction by the NA substrate. A representative example of the simulation box after melt equilibration is shown in Fig. 2a, and the nucleation and growth of a crystalline phase advancing in time is shown in Fig. 2b–d.
Figure 2
Representative simulation box containing C50 (black atoms) bounded by substrate on two sides (blue and white atoms). The simulation box is periodic in all directions. In (a) all atoms are displayed, while in (b–d) only substrate and crystalline C50 atoms are displayed, with snapshots taken at increasing amounts of time from the start of the quench to Tc: (b) 1 ns, (c) 10 ns, (d) 70 ns.
×
Simulations were performed in LAMMPS [82]. All simulations were performed in the NPT ensemble at 1 atm, maintained by a Nosé–Hoover thermostat and barostat (damping frequencies: ωT = 1/100 Δt, ωP = 1/1000 Δt, Δt = 3 fs). The simulation box was periodic in all dimensions, and each dimension was varied independently. Forces were integrated with a rRESPA multi-timescale integrator, with an inner time step of 1.5 fs for bond length, bond angle, and bond dihedral interactions, and an outer time step of 3 fs for LJ and SW interactions.
The equilibrium melting temperature Tm for C50 using the PYS force field is 370 K (97 °C) [69]. Each C50+NA system was first equilibrated at 500 K (227 °C) for 20 ns. A simulation of 10 ns was sufficient to completely melt the starting crystalline configuration for C50 [69], and configurations were saved every 1 ns subsequently, resulting in 10 independent melt configurations for each C50+NA system, based on an estimated Rouse time τR < 1 ns [13]. Crystallization simulations were then run by quenching each melt configuration to Tc = 360 K (87 °C), for a supercooling (1 − Tc/Tm) of 2.7%.
To quantify crystallization, the crystalline growth front was monitored as a function of time, following the methods of Bourque et al. [12, 13]. Specifically, UAs were determined to be crystalline or amorphous using the local p2 order parameter given by
where θij is the angle between the vectors spanning the midpoints of the bonds to either side of the ith and jth atoms, and the average was taken over all UAs within a cutoff distance of 2.5σLJ from the ith UA; UAs with p2 > 0.4 were marked as crystalline, following the analysis of Yi and Rutledge [83] and in accordance with previous work on heterogeneous nucleation in C50 [12, 13]. The C50 domain was sliced into 0.4 nm thick layers normal to the direction of growth, starting at the interface with the NA substrate and moving outwards into the C50 domain. Layers of this thickness correspond to molecular layers within the C50 crystal [69]. The crystallinity within each layer, X(z), was calculated as the number of crystalline atoms divided by the total number of atoms within the layer. The location of the crystal growth front, D, was calculated using a Gibbs dividing surface:
Following the methods of Santana and Müller [84], a mixture of HDPE and NA was melt-blended with PS using a twin blade blender (Brabender ATR Plasticorder with 3-piece mixer) under the same conditions reported previously [22]. Specifically, each NA was mixed with HDPE at three loadings (0.5, 2, 5 wt%) at 180 °C and 100 rpm for 10 min. This mixture (HDPE+NA) was then melt-blended at 15 wt% HDPE with PS at 180 °C and 100 rpm for 30 min. This final blend (PS/HDPE+NA) was hot pressed (Carver Model C) at 115 °C (which is above the glass transition temperature, Tg, of PS but below the melting temperature, Tm, of HDPE) into a sheet roughly 1 mm thick. A 3 mm diameter disc was punched out for thermal analysis with DSC. Each sample was prepared in duplicate.
Samples were analyzed by DSC (TA Instruments DSC 2500) using two procedures: nonisothermal crystallization at constant cooling rate (10 °C min−1) from the melt (equilibrated for 2 min at 180 °C) and isothermal crystallization at discrete crystallization temperatures (Tc = 121.5–124 °C at 0.5 °C increments) following a rapid quench from the melt (30 °C min−1). A PS sample of similar thermal mass to the test sample was used as the reference in all experiments to reduce the effect of an endothermic hook upon rapid quenching to Tc, as noted elsewhere [22]. The exothermic heat flow of the sample was used to quantify crystallization.
Results
Simulated nucleation
Figure 3a–c shows trajectories of the crystal growth front, D, for each NA, averaged over all replicates. After a short period of slow rearrangement leading to organization of C50 chain segments near the interface (below the dotted line), a stable crystallite forms and subsequently grows at a steady rate. The linear growth regime was extrapolated backwards to z = 0 nm to identify an induction time for nucleation, τ [12, 13], as shown by the dashed lines in Fig. 3a–c. The induction time for each NA is shown in Fig. 3d, with the corresponding nucleation rate (1/τ) being shown on the right ordinate axis. The interfacial area between C50 and substrate was the same, ~ 41 nm2, for all of the nucleation simulations. From this data, the nucleation rate is highest for hBN (τ = 6.3 ± 1.6 ns), followed by MoS2 (τ = 14.1 ± 3.3 ns), and lowest for WS2 (τ = 19.2 ± 2.6 ns). The average growth rate, G, for all NAs was 0.033 ± 0.002 nm ns−1, which compares well with previous measurements of the growth rate of C50 in the presence of tetrahedral NAs (G = 0.031 nm ns−1) [12] and graphene-like (G = 0.043 nm ns−1) [13] NAs.
Figure 3
Trajectories of crystal growth front (D) from molecular simulations, averaged over 10 independent starting configurations (2 growth front trajectories per configuration): (a) hBN, (b) MoS2, (c) WS2; (d) induction times (left ordinate axis) and nucleation rates (right ordinate axis) calculated from these trajectories (error bars denote the 95% confidence interval). (a–c) Averages (solid lines) are taken over 20 growth front trajectories, and the grey region in each plot denotes the 95% confidence interval around the average. The dashed lines denote the fit to the linear growth regime, and the dotted lines denote the lower and upper bounds for fitting of the linear growth regime. The lower bound is limited by subcritical fluctuations, while the upper bound is limited by impingement with the growth front propagating from the opposing interface in the simulation. The induction time, τ, is found at the intersection of the fitted line and the time axis.
×
As an aside, we note that roughly 20–25% of the trajectories for MoS2 and WS2 exhibited the formation of a second large nucleus at an orientation that could not readily merge with the first, or largest, nucleus, resulting in an abnormal slowdown of the growth front. This behavior is manifested in the observation of larger confidence intervals about the average progression of the crystal/melt interface in the growth regime for these two NAs in Fig. 3b and c. Excluding these trajectories from the analysis results in a modest 5–10% reduction in induction time, which does not affect any of our conclusions, and a reduction in the confidence intervals about the growth front to similar magnitudes as observed with hBN. This behavior is not seen with hBN (Fig. 3a) due to the more rapid spreading of the crystalline cluster on the surface and completion of a crystal layer before a second heterogeneous nucleation event occurs.
Thermodynamically, the interfacial free energy difference, Δσ, is responsible for the decreased free energy barrier to heterogeneous nucleation relative to homogeneous nucleation. This parameter is of particular significance for comparing the effectiveness of different NAs, and has been used to define the thermodynamic efficiency of a polymer plus NA pair [22],
$$E = 1 - \Delta \sigma /\sigma$$
(8)
where σ is the analogue to Δσ for homogeneous nucleation. In this definition, the best NAs yield efficiency values near 1.
From classical nucleation theory [85], a nucleus in the shape of a rectangular prism has a free energy, ΔG, given by
where a is the width of the nucleus, l is the length of deposited chain stems, b is the height of the nucleus, and σ and σe are the interfacial free energies for the lateral and stem end surfaces of the C50 crystal in contact with the melt. Δgv is the bulk free energy given by Δgv = ΔhfΔT/Tm0, where the experimental value of Δhf for C50 reported by Hammami et al. (185 kJ/mol) [86, 87], the experimental density of C50 at 93 °C reported by Sinnatt (0.794 g cm−3) [88, 89], and the melting temperature of C50 using the PYS force field reported by Bourque et al. (370 K) [69] were used to calculate Δgv. Δσ is more specifically defined as σNA-crystal + σcrystal-melt − σNA-melt and is unique to each polymer plus NA pair.
Setting derivatives of ΔG with respect to each dimensional variable to zero gives the dimensions of the critical nucleus,
where A* is the area of the critical nucleus. A* can be estimated by the area of the first layer of crystalline C50 UAs at the induction time and should be independent of substrate, according to Eq. (10) where A* is not a function of Δσ. b* was calculated as the height of a prism of equivalent volume to the actual nucleus observed in simulation (i.e. b* = V*/A*, where V* is the volume of the critical nucleus measured in simulation). By estimating b* in this way, Δσ can be calculated from Eq. 10; the results are shown in Fig. 4. By this analysis, MoS2 is the most efficient (lowest Δσ) NA, followed by hBN and WS2, which are approximately equivalent.
Figure 4
Interfacial free energy difference, Δσ, calculated from critical nucleus sizes measured by molecular dynamics simulations. Data were averaged over 20 trajectories, and error bars denote the 95% confidence interval.
×
Experimental nucleation
The temperature range of interest for crystallization of HDPE in the presence of NAs was determined using nonisothermal crystallization experiments. The exotherms from DSC are shown in Fig. 5. Crystallization was observed to occur around 120 °C for all samples containing NAs. We have previously shown that samples produced using these methods with no added NAs exhibit crystallization peaks near 80 °C, consistent with prior work in the literature [22, 62, 90]. The upward shift of the peak by 40 °C can be attributed to the presence of NAs. Crystallization in this temperature range is also consistent with other crystallization experiments reported for HDPE+hBN and HDPE+MoS2 nanocomposites [29, 45].
Figure 5
Crystallization of samples at a constant cooling rate from the melt. Data are organized by NA: hBN (blue), MoS2 (orange), WS2 (green). Percentages after NA type indicate the loading of NA in the sample. Curves have been vertically shifted to allow for comparison.
×
The presence of two high temperature peaks between 110 and 120 °C is noticeable in the DSC traces for almost all of the samples; this behavior was also observed in similar experiments with HDPE micro-domains nucleated in the presence of GN nanoplatelets [22]. The lower temperature peak is more prominent in samples with lower loadings of NAs, while the higher temperature peak becomes more prominent as loading is increased. The more effective nucleation process beginning above 120 °C was selected for in-depth analysis using isothermal crystallization experiments.
Isothermal crystallization experiments were run at six crystallization temperatures (Tc). The lower bound of Tc is limited by the onset of rapid crystallization before the sample can be fully cooled to Tc. The upper bound of Tc is limited by crystallization that is too slow to register with the DSC instrument used, without resorting to more exotic methods (e.g. isothermal step crystallization) [60, 61]. The evolution of crystallinity as a function of time is shown in Fig. 6 for hBN-0.5%, while the data for the remainder of samples are shown in Fig. S1 of the Supplemental Information. Crystallinity was calculated as the cumulative heat flow, Δh(t), scaled by the total enthalpy of crystallization, given by
$$X = \frac{\Delta h\left( t \right)}{{\int\limits_{0}^{\infty } {\Delta h\left( t \right){\text{d}}t} }}$$
(11)
Figure 6
Evolution of crystallinity at different Tc for sample hBN-0.5%. Each curve is the average of two replicates, and the error bars denote the high and low values of the replicates at evenly spaced intervals.
×
Avrami model
The evolution of crystallinity was first analyzed using the Avrami equation:
where t0 is the induction time prior to the onset of crystallization, K is the overall crystallization rate, and n is the Avrami index, consisting of two components, n = nG + nN. nG is the dimensionality of growth, and nN is the dimensionality of nucleation. nG typically takes values between 0 (no growth) and 3 (growth in 3 independent directions). nN is often found to be a number between 0 (instantaneous nucleation) and 1 (sporadic nucleation), with 0.5 corresponding to a special case of diffusion-limited nucleation. Non-integer values of nN are often explained as having a combination of both sporadic and instantaneous nucleation [91].
Equation (12) was fit to the data for crystallinity versus time in Fig. 6 and Fig. S1 for conversions between 0.03 and 0.20. Secondary crystallization effects are insignificant in this region, allowing the most accurate characterization of primary crystallization processes [92]. The fitted Avrami indices n for all samples are shown in Fig. 7. Avrami indices between 1 and 2 are observed, with no apparent dependence on Tc. While a nucleation-limited process would be expected to exhibit an Avrami index near 1 [59], similar indices of 1 to 2 were reported for experiments with HDPE microdomains nucleated in the presence of GN nanoplatelets [22]. The argument for a nucleation-limited process is examined further in the Discussion.
Figure 7
Avrami index n, obtained from fitting of the DSC data for crystallization vs time to the Avrami equation, Eq. (12), in the range 0.03 < X < 0.20. Data was averaged over both replicates and all Tc for each NA loading. Error bars denote the 95% confidence interval.
×
First-order nucleation model
For a nucleation-limited process, the Avrami equation can be approximated as
where J0 is a constant related to the number of potential nucleation sites and is relatively temperature insensitive. In the second term, U* is the activation energy for polymer chain reptation in the melt, R is the universal gas constant, and T∞ is the temperature at which all reptation stops, approximately 30 K below Tg [94]. In the third term, Tm0 is the equilibrium melting temperature for PE, k is the Boltzmann constant, ΔT = Tm0 − Tc is the degree of supercooling below Tm0, and Δhf is the enthalpy of fusion. f is a correction factor to account for the temperature dependence of Δhf at large ΔT [95], given by
$$f = \frac{{2T_{c} }}{{T_{c} + T_{m}^{0} }}$$
(15)
Values reported in the literature for parameters in these equations are shown in Table 4.
The nucleation rate as a function of the thermal driving force, as described by Eq. (14), is plotted in Fig. 8 for all samples. For all NA types, trends observed are in line with expectation: increasing NA loading increases the nucleation rate, while increasing the temperature decreases the nucleation rate. Comparing between NA types, hBN exhibits the highest nucleation rates, while the nucleation rates for MoS2 and WS2 are roughly the same.
Figure 8
Nucleation rates according to classical nucleation theory. Symbols denote the average of two replicates at each NA loading, and error bars denote the high and low values; in most cases, the difference is smaller than the size of the symbol itself.
×
According to nucleation theory, the slope of the line through each of the sets of data in Fig. 8 is directly related to Δσ. Because the size of the critical nucleus cannot be measured with the experimental methods here, this analysis provides an alternative to Eq. (10) to calculate Δσ. Values for Δσ and J0 from linear regression on this data are shown in Fig. 9. For each NA type, the calculated values for Δσ agree very well between 2 and 5% loading, with a small deviation at the lowest loading; averaging over 2–5% loadings, the Δσ values (and corresponding thermodynamic efficiencies) are as follows: hBN: 0.943 ± 0.050 erg cm−2 (E = 0.902), MoS2: 0.811 ± 0.085 erg cm−2 (E = 0.916), WS2: 0.944 ± 0.078 erg cm−2 (E = 0.902); for 0.5% loading, the efficiencies are slightly lower: EhBN = 0.895, EMoS2 = 0.904, EWS2 = 0.901. Comparing between NA types, MoS2 has the lowest value for Δσ and thus the highest nucleation efficiency (i.e. lowest barrier to nucleation); however, the difference between NA types is not large enough to be statistically significant. Comparing J0 across NA types, hBN > WS2 > MoS2, without any clear trend as a function of NA loading. These results are examined further in the Discussion section.
Figure 9
(a) Fitted values of interfacial free energy difference, Δσ (left axis), and nucleation efficiency, E (right axis, reversed axis); (b) nucleation prefactor, J0. Error bars denote the 95% confidence interval.
×
Discussion
Simulated nucleation
Effects of force field parameters
In their simulations of heterogeneous nucleation of C50 on GN-like materials, Bourque and Rutledge found an induction time of 14 ± 2 ns for C50+GN for a 50 nm2 substrate [13]. Compared to the results above, after scaling the induction time for nucleation on a slightly larger GN substrate, we find that τhBN < τMoS2 < τGN < τWS2. Bourque and Rutledge also identified correlations between induction time and SW force field parameters: εSW, λSW, and σSW. εSW and λSW describe two- and three-body interaction strengths, respectively, and σSW describes atomic spacing. Varying these three SW parameters was sufficient to vary uniquely the rigidity of the NA substrate and the degree of epitaxial match to the C50 crystal structure, as all the other parameters in the SW force field were unchanged. However, in the parameterization for TMDs used in this work, other SW parameters were changed in addition to εSW, λSW, and σSW. Therefore, to better quantify two- and three-body interaction strengths and atomic spacing between different NAs, parameter combinations that appear in Eqs. (2) and (3) are used here. Specifically, (AεSW) is the prefactor in Eq. (2) for two-body interactions, (λSWεSW) is the prefactor in Eq. (3) for three-body interactions, and (aSWσSW) is the atomic spacing that appears in the final term of Eq. (2).
These alternative definitions are used to compare NAs relative to GN in Fig. 10, using the parameter values reported in Tables 2 and 3. Relative values are given as x* = xNA/xGN, where x is one of the foregoing force field parameter combinations. In both two- and three-body interactions, the same trend is observed, with GN being the most rigid substrate, followed by hBN, MoS2, and WS2 being the softest. hBN has a similar atomic spacing relative to GN, while both TMDs have larger spacings between atoms in the substrate. According to general trends noted by Bourque and Rutledge, weaker two- and three-body interactions result in increased lattice flexibility that better accommodates mismatches in epitaxy between the NA and C50; however, if these interactions are too weak, the two- and three-body interactions hinder crystallization because thermal fluctuations in the NA substrate become too strong [13]. This compromise is reflected by the simulation results in this work, with a reduced induction time for hBN relative to GN (attributed to flexibility accommodating lattice mismatch) as predicted by Bourque and Rutledge, but a similar or slightly larger induction time for TMDs (attributed to poorer lattice matching and excessively strong thermal fluctuations).
Figure 10
Comparison of force field parameters, (a) two-body interaction strength for NAs (AεSW), (b) three-body interaction strength for NAs (λSWεSW) where three-body interactions for TMDs are defined as in ref. 53, (c) atomic spacing of NAs (aSWσSW). Values for each are scaled relative to GN, e.g. (AεSW)* = (ANAεSW,NA)/(AGNεSW,GN); AGN = 5.8341, εSW,GN = 125.42 kcal mol−1, λSW,GN = 40, σSW,GN = 1.28 Å, aSW,GN = 1.8 [18].
×
A fourth parameter that was examined in the parametric sweep of Bourque and Rutledge was εAD. εAD quantifies the strength of the intermolecular interactions between C50 and the NA using LJ interactions, calculated by εAD,X = (εLJ,C50εLJ,X)0.5, where X is an atom in the NA. The combined interaction strength for the NA, εAD, was averaged over the atoms in the NA according to atomic ratios. A comparison is made in Fig. 11a, where we find that εAD,GN < εAD,hBN < \(\varepsilon_{\text{AD,MoS}_{2}}\) ~ \(\varepsilon_{\text{AD,WS}_{2}}\).
Figure 11
Comparison of force field parameters, (a) interaction strength between NA and C50 (εAD), and (b) areal density of the energy of interaction (ΦAD/AI) between NA and C50 at the distance of the first crystalline layer of C50. Values for each are scaled relative to GN, e.g. εAD* = (εAD,NA) / (εAD,GN). In (a), εAD are averaged over the atoms in the NA according their atomic ratios, e.g. εAD,MoS2 = (εAD,Mo + 2εAD,S)/3. The values for εAD and (ΦAD/AI)* for GN are 0.1034 kcal mol−1 [18], and − 2.566 kcal mol−1 nm−2, respectively.
×
However, this comparison oversimplifies the energetic landscape between C50 and the NA substrates. The density of atoms is far greater in each hBN sheet compared to the two TMD sheets due to the smaller atomic spacing in hBN (Fig. 10c). Additionally, two hBN sheets are within the LJ cutoff distance of the first layer of C50, while only a single TMD sheet is within this cutoff distance, due to the thicker structure of TMD sheets and the larger distance between sheets (c.f. σLJ in Table 1 for inter-sheet interactions set by experimental inter-sheet spacing). Therefore, the total interaction potential between all NA atoms and C50 UAs, ΦAD, at the distance of the first crystalline layer of C50 was calculated using
where ELJ,k is the LJ potential given by Eq. (5) between the kth NA atom and C50, evaluated at the distance between that atom and the first crystalline layer of C50 UAs, Rk. ELJ,k is only evaluated for Rk < rcutoff, LJ. This total potential strength was scaled by the interfacial area, AI, to obtain an areal density of the energy of interaction between the NA and C50 in the first crystalline layer. The results are shown in Fig. 11b, where a stark difference is observed between hBN and the TMDs: hBN > GN > MoS2 ~ WS2. Bourque and Rutledge found that increasing εAD reduced induction time, with no observed extremum [13]. This negative correlation is consistent with the induction times observed for the NAs here, where hBN exhibits the strongest NA-C50 interactions and the shortest induction time, while WS2 has the weakest NA-C50 interactions and the longest induction time.
While all three of these NAs have SW parameters that are outside the bounds of the parametric sweeps of ref. [13], a few inferences can be made based on extrapolated trends from those earlier results. For comparison here, the alternative definitions using combined parameters for atomic spacing and two-body interaction strengths are used as surrogates for σSW* and εSW* from ref. [13], respectively. For all three NAs, εAD* ≈ 1.2, so the induction times of these NAs are compared in Fig. 12 to the results reported in Fig. 2a of ref. [13]. Both TMDs fall in a region of (AεSW)* < 0.25, (aSWσSW)* > 1.35 where induction times appear to be increasing strongly with distance from the minimum; assuming this trend in induction times continues for lower (AεSW)*, it is expected that τTMD ≥ 20 ns [13]. In contrast, hBN falls in line with the extended “valley” between 1.00 < (aSWσSW)* < 1.10 for all values of (AεSW)* > 0.30, with an extrapolated τ ~ 9 ns [13]. Despite the crudeness of the extrapolations, these estimates for induction times based on the data from ref. [13] are in rough quantitative agreement with the induction times calculated for the three NAs tested here.
Figure 12
Comparison of SW parameter values for hBN (blue X), MoS2 (orange X), and WS2 (green X) to parametric studies of GN-like NAs with εAD* = 1.2 (ref. [13], Fig. 2a). Induction times for hBN, MoS2, and WS2 are from MD simulation results reported in Fig. 3d. Colored portion of plot reprinted from ref. [13], Copyright 2018, with permission from Elsevier.
×
Epitaxy
Unlike the case for the experimental measurements, epitaxy is easily quantified by simulation. Hexagonal unit cells have been reported for these NAs, with the following lattice parameters: ahBN = 2.50 Å [79], aMoS2 = 3.16 Å [80], and aWS2 = 3.15 Å [81] (c.f. Figure 1). For C50, crystals always form in layers that are spaced about 4 Å apart, which corresponds well to the distance (aC50/2) between layers of chains in the pseudo-orthorhombic unit cell of C50 when simulated using the PYS model; thus, the inter-chain spacing within the first layer that is relevant for epitaxy is the unit cell length bC50,PYS = 4.9 Å [12, 13].
Unit cell parameters at the interface were measured in simulations using the radial distribution function for NAs (between interfacial B-B atoms in hBN, and between interfacial S–S atoms in TMDs, Fig. 13a), and by averaging the distance between the longest C50 stem and its neighbors in a crystal cluster, bC50,obs (Fig. 13b). The lattice strain in C50 at the interface was then calculated as εstrain = (bC50,obs − bC50,PYS)/bC50,PYS × 100%; lattice strain in the NA was calculated similarly using aNA and the location of the first peak in the radial distribution function. These strains are shown in Fig. 13c for C50 and the NAs. Differing amounts of strain in the C50 unit cell are observed for each C50 plus NA system, but in all cases the C50 strain is larger than that in the NA and compressive in nature. C50+hBN exhibits the largest strain (εstrain = − 8%), followed by C50+MoS2 and C50+WS2 (εstrain = − 7%). The magnitude of the fluctuations in the lattice strain for NAs is captured by the error bars in Fig. 13c, where TMDs exhibit ~ 40% larger fluctuations than hBN. This agrees with the earlier hypothesis that TMD lattices suffer from excessively strong thermal fluctuations that hinder crystallization, due to weak two- and three-body interactions (Fig. 10a, b).
Figure 13
Calculation of unit cell parameters from simulation. (a) The radial distribution function for each NA, averaged over the final 10 timesteps, where aNA is calculated as the first peak in g(r). (b) Example calculation of bC50 from inter-chain distance between the longest crystalline stem and its neighbors in the crystal for a single MD trajectory (blue). The induction time, τ (dotted red line) and the average for bC50 for t > τ (dashed black line) are also shown. (c) Strain, εstrain, calculated as deviation of unit cell parameters measured in simulation from the reported values for C50 (PYS, dark grey) and NA (experimental, light grey) crystals; error bars denote one standard deviation.
×
For these hexagonal NAs, there are two significant crystallographic spacings: The first is that between parallel “zigzag” (ZZ) directions, and the second is that between parallel “armchair” (AC) directions, as shown in Fig. 1. The spacing between parallel ZZ segments, dZZ, is given by 0.866aNA, while the spacing between parallel AC segments, dAC, is given by aNA. Epitaxial mismatches are quantified as a percentage using
where dNA can be either dZZ or dAC for a NA [6, 100, 101]. Because the unit cells of the NAs are significantly smaller than the inter-chain spacing of C50, supercells of the NAs exhibit the best epitaxial matching with C50. mC50 and mNA in Eq. (17) are integers corresponding to the number of unit cells required in the supercells to get the best epitaxial matches (smallest Δ) between the crystal lattices of C50 and the NA. For hBN, a 2:1 supercell combination produces the best epitaxial match, while a 3:2 supercell combination is best for both TMDs. Figure 14a shows the nominal values for Δ for two supercell combinations (2dNA: dC50, and 3dNA: 2dC50), calculated using bC50,PYS and the literature values for aNA, while Fig. 14b shows the values for Δstrained using the same supercells for the NAs but the interfacial strained value for C50, bC50, obs. The range of acceptable epitaxial mismatch is generally considered to be |Δ|≤ 15% [6], demarcated by the dashed black lines in Fig. 14. Based on the nominal value for Δ, the AC direction should be favored for all NAs; however, in simulation, neither orientation is strongly favored over the other (50–60% of nucleation events were aligned along ZZ direction for each NA). This observation is consistent with the values for Δstrained, where epitaxy is satisfied for both ZZ and AC directions.
Figure 14
Lattice mismatches between NA and C50 crystals: (a) Δ, the nominal epitaxial mismatch using literature values for aNA (experimental) and bC50,PYS, and (b) Δstrained, the lattice mismatch using strained unit cells observed at the interface in simulation. Lattice matching along the zigzag (ZZ) direction are shown in red, while blue corresponds to the armchair (AC) direction. A 2:1 supercell (2dNA: dC50) is used for hBN, while a 3:2 supercell (3dNA: 2dC50) is used for TMDs.
×
Even though the TMDs exhibit the lowest Δstrained, hBN exhibits the lower induction times. This discrepancy suggests that factors other than epitaxy are at work. In a study of heterogeneous nucleation of acetominophen, Chadwick et al. found that registration of specific surface functionalities between acetaminophen and crystalline substrates was more important in reducing induction time than epitaxial matching [102]. The importance of specific chemical interactions was also noted by Olmsted and Ward with other small organic crystals [103]. Although C50 lacks any specific chemical interactions with the NAs studies here, nevertheless the density of interactions mentioned earlier gives rise to a stronger C50–NA attraction in hBN that is at least as important as epitaxy in determining the efficiency of the NA.
Experimental nucleation
For a nucleation-limited crystallization process, the Avrami index should be near 1, and the Turnbull number, tgrow/tnucl, should be much less than 1 [59]. Here, tnucl is the characteristic time to nucleate a crystallite within a micro-domain, tnucl = 1/J, and tgrow is the characteristic time for the resulting crystallite to grow to fill the microdomain, tgrow = d/(2G), where d is the nominal size of the NA (c.f. Materials section) and G is the crystal growth rate of the HDPE resin used in this work [22]. As was the case for similar experiments with GN as the NA [22], Avrami indices 1 < n < 2 were observed for all NAs in this work. However, the Turnbull numbers for these samples range from 10–5 to 10–4. This calculation implies that the crystallization process is indeed nucleation-limited. The nucleation-limited nature of the process was confirmed qualitatively using polarized optical microscopy, where induction times on the order of seconds to minutes were observed, but growth appeared instantaneous to the eye (< < 1 s) [22].
Another interesting feature in the DSC data (Fig. 5) is the presence of two crystallization peaks observed during non-isothermal crystallization experiments, characterized by a transition in prevalence from a high-temperature peak at the two higher NA loadings (2 and 5%) to a second peak 2–3 °C lower at the lower NA loading (0.5%). The role of impurities in the HDPE can be eliminated based on nonisothermal crystallization experiments in the absence of NA, where nucleation does not take place until the system has cooled to nearly 80 °C [22]. The role of impurities in the NAs themselves is unlikely, given the certificates of 99% purity provided by the manufacturer; additionally, the concentration of such impurities should be in proportion to the NA, so no difference would be observed with decreasing NA concentration. We hypothesize that this bimodal crystallization behavior is a consequence of NA aggregation: as the loading increases, the potential for contact between and aggregation of NA platelets in the samples increases. Lack of perfect registration between the particles in such aggregates could create edges that promote a grapho-epitaxial mechanism and nucleation at a higher temperature. In samples with low NA loadings, such aggregation would be less common, resulting in a smaller high temperature peak. Further evidence for aggregation of NA particles is given by the lack of trends in J0 as a function of NA loading (Fig. 9b). J0 should be positively correlated with NA loading, but this is only weakly observed with WS2, and not observed for either hBN or MoS2.
This difference in mechanism also may explain the difference in thermodynamic nucleation efficiency, E, that is observed between the higher NA loadings (2, 5%) and the lower NA loading (0.5%). For 2% and 5% loading, the thermodynamic efficiency could be higher due to the grapho-epitaxial mechanism for heterogeneous nucleation, while the thermodynamic efficiency for 0.5% loading could be lower due to a significant percentage of the sample nucleating by conventional epitaxy. Interestingly, WS2 did not exhibit lower E at 0.5% loading than at 2% and 5% loadings, but still shows the two-peak behavior during non-isothermal crystallization.
Compared to other NAs for PE, the NAs studied here perform well (E = 0.902–0.916). Using the same HDPE resin with GN nanoplatelets serving as the NA, a comparable value of E = 0.913 was reported previously [22]. These values for thermodynamic nucleation efficiency are relatively high, indicating these materials are all effective as NAs for HDPE when compared to many other polymer plus NA pairs that have been published in the literature [22].
Comparison between simulations and experiments
The nucleation rates observed experimentally are orders of magnitude slower than those calculated in silico. This discrepancy is readily explained by the significantly faster dynamics of C50 compared to experimental samples of HDPE, where entanglements play a dominant role. Nucleation is an activated process, described classically by both a pre-factor and an activation energy. The dependence of the pre-factor on molecular weight in experiments is well documented, through its influence on viscous or diffusive processes, whereas the activation energy typically is not dependent on molecular weight [104, 105]. This effect on the comparison between experiments and simulations has been discussed previously [106]. For this reason, a comparison of absolute values for nucleation rate between HDPE and C50 is not advised. Additionally, we also do not have accurate values for the interfacial area of the NAs in the experimental work, as would be required to compare nucleation rates in the appropriate units of s−1 m−2. Nevertheless, some qualitatively similar trends are observed between simulations and experiments. Figure 15 shows relative nucleation rates from both simulation and experiment; agreement is not quantitative, but the trend of \(J_{\text{hBN}} > J_{\text{MoS}_{2}} > J_{\text{WS}_{2}}\) is observed in both cases.
Figure 15
Nucleation rates from experiments and simulations, relative to hBN. Experimental nucleation rates were calculated as JNA/JhBN averaged over all Tc.
×
By contrast, Δσ contributes to the activation energy for nucleation and is therefore not strongly affected by molecular weight. A comparison can be made between Δσ values that were calculated from experiments and those that were estimated from simulations using critical nucleus sizes. The values for Δσ from simulation are approximately double the values calculated from experiments (compare Figs. 4 and 9a). There are several plausible explanations for this discrepancy, including the difference in methods employed to evaluate Δσ. A second possible explanation is the onset of a grapho-epitaxial mechanism hypothesized to occur in the experiments, whereas idealized flat substrates were used for simulation. While the estimates from simulation are fairly noisy, qualitatively similar trends are observed in both experiments and simulations: ΔσMoS2 < ΔσhBN ~ ΔσWS2. This result provides a measure of confidence in the simulations and their ability to predict accurately the relative nucleation rates and interfacial energies in the real materials.
Despite MoS2 having the smallest critical free energy barrier to nucleation (ΔG* ~ Δσ), hBN demonstrates the highest nucleation rates, observed both in experiments and simulations. The kinetic prefactor (J0 in Eq. 14, Fig. 9b) is related to the concentration of nucleation sites. hBN has a smaller crystal unit cell than both TMDs, providing more sites for nucleation to occur per unit area. The higher density of interactions is reflected by the larger experimental value of J0. Since the free energy barriers are very low and roughly the same across the three NAs, the concentration of nucleation sites and resulting stronger interactions are the significant factors that result in JhBN > JMoS2 ~ JWS2.
Conclusions
A combination of molecular dynamics simulations and experimental measurements are used here to characterize heterogeneous nucleation rates for HDPE on hexagonal 2D nucleating agents. Despite several limitations of the simulations, good agreement is found between simulations and experiments regarding qualitative trends in relative nucleation rates and relative values of the interfacial free energy difference, Δσ. Of the NAs examined, hBN exhibits the highest nucleation rate in both simulations and experiments, in agreement with earlier predictions from simulations of GN-like nucleating agents. Simulations furthermore permit a more detailed examination of the mechanisms underlying heterogeneous nucleation. Strain in the organic crystal permits improvement in epitaxy and enables alignment of the chains along either the zigzag or armchair directions of the substrates, both of which facilitate nucleation. The strength of interaction between the nucleating agent and the oligomer was also found to be a key factor in promoting nucleating agent efficiency, with stronger interactions positively correlating with nucleation rate. With this validation of the predictive capabilities of molecular simulations for polymer nucleation, simulations can be used with greater confidence to examine other nucleating agents or to design novel nucleating agents.
Acknowledgements
Financial support for this work was provided by the Designing Materials to Revolutionize and Engineer our Future (DMREF) and Grant Opportunities for Academic Liaison with Industry (GOALI) programs of the National Science Foundation, Division of Civil, Mechanical and Manufacturing Innovation, Award #1729304. This work made use of the Institute for Soldier Nanotechnologies Shared Experimental Facilities at MIT.
Declarations
Conflict of interest
The authors declare no competing interests or relationships that could influence or bias this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
In order to compare to previous work with graphene-like materials [13], a force field for bulk hBN was desired that uses the Stillinger–Weber (SW) potential for intra-sheet interaction and the Lennard–Jones (LJ) potential for inter-sheet interactions. In their work on graphene (GN), Bourque and Rutledge optimized parameter values in the force field to reproduce experimental elastic constants for GN by molecular simulations [18]. In this work, many of the parameters of the SW potential were kept the same as those used previously for both GN [18] and silicene [107]. To model the elastic constants for hBN, the inter-sheet interaction strength, εLJ, the intra-sheet two-body interaction strength, εSW, and the intra-sheet three-body interaction strength, λSW were optimized. The inter-sheet distance, σLJ, and the intra-sheet bond length, σSW, were set by the experimental values (dinter-sheet = 3.33 Å, dbond = 1.45 Å, σ = d/21/6) [108], and the bond angle, θ0, was fixed at 120° to replicate the symmetry of the hBN lattice.
hBN was built as 6 × 6 × 2 repeats of the unit cell structure (4 layers) [79]. The a- and c-axes of the unit cell were aligned with the Cartesian x- and z-axes, respectively, of the simulation box, and periodic boundary conditions were employed in all Cartesian directions. Simulations were performed in LAMMPS [82]. hBN was equilibrated from its initial configuration in the NPT ensemble (300 K, 1 atm), maintained by a Nosé-Hoover thermostat and barostat (damping frequencies: ωT = 1/100Δt, ωP = 1/1000Δt, Δt = 2 fs). Each dimension was allowed to vary independently. A Verlet time step integrator with a time step of 2 fs was used to integrate the forces. After equilibration, the simulation box was deformed along the x-, y-, or z- dimension at a constant engineering strain rate of 107 s−1, while the orthogonal dimensions were held fixed. Temperature was held constant at 300 K by the thermostat, with velocity contributions due to deformation removed using the “remap x” command. The stress as a function of strain was used to calculate elastic stiffnesses over the linear regime at low deformation (0–5%). A linear model was fit to the simulated elastic stiffnesses, and a new set of parameter values was proposed to reduce the error between simulated elastic constants and the experimental values. This process was iterated several times until simulations reproduced the experimental values. Three independent simulations using the final parameter set were used to characterize the variances in elastic constants.
Results
The stress response as a function of strain for the final parameter set is shown in Fig.
Figure 16
Stress, σ, due to deformation of hBN nanoplatelets as a function of strain, e. (a) Stress along the x-axis (σ1, left axis, black curves) and the y-axis (σ2, right axis, grey curves) for deformations along the x-axis (e1) and along the y-axis (e2). (b) Stress for deformation along the z-axis.
16. As was noted previously for GN [18], armchair and zigzag directions are also present in hBN, and lead to some anisotropy between C11 and C22, and between C12 and C21. Therefore, the simulated values were averaged for comparison to the experimental values, and ambiguity between stiffness and modulus for the latter was disregarded. A comparison of experimental vs. simulated elastic constants is made in Table
Table 5
Comparison of elastic constants determined experimentally with those reproduced through simulation using the combined SW-LJ force field
Die im Laufe eines Jahres in der „adhäsion“ veröffentlichten Marktübersichten helfen Anwendern verschiedenster Branchen, sich einen gezielten Überblick über Lieferantenangebote zu verschaffen.
