Nonequilibrium Molecular Dynamics Research Articles

Orientation effects on shock-induced plastic deformation in FeNiCoCu high entropy alloy

FeNiCoCu high-entropy alloys (HEAs) demonstrate promising potential for widespread use in structural and functional applications. However, a thorough understanding of dynamic deformation processes in FeNiCoCu HEA is limited due to technological constraints in detecting real-time microstructural developments at the atomic level. This study examines the shock-induced plastic deformations in the equiatomic FeNiCoCu HEA, focusing on crystallographic orientation and particle velocity, using nonequilibrium molecular dynamics simulations. We obtained the P−V/V0, P−T, P−Up, and Us−Up Hugoniot relations and evaluated their anisotropy. The shock velocity, stress, and shear stress exhibit orientation dependence due to the differences in the plastic deformation mechanism. For shock loading along [100] orientations, dislocation dominates at lower shock intensities. However, a phase transition from face-centered-cubic (FCC) to body-centered-cubic becomes the primary plastic deformation at high shock intensity. For shock loading along [112¯] and [111] orientations, the generation of disordered structures and dislocation activities is revealed to play an important role in the development of localized plastic deformation. Moreover, the competition of disordered and hexagonal-close-packed (HCP) atoms is observed. The transition from FCC to disordered atoms provides nucleation sites for dislocations, and the slip of dislocations around disordered atoms leads to the formation of HCP structures. These findings are very helpful for learning the dynamic deformation behavior of FeNCoCu HEA.

Investigation of thermal transport mechanism of silicone-modified phenolic matrix nanocomposites with different pyrolysis degrees

Polymeric nanocomposites with low thermal conductivity show promising applications for next-generation thermal protection materials used in re-entry vehicles due to their lightweight, high char yield, and excellent ablation-oxidation resistance. However, the thermal conductivity of polymeric nanocomposites varies with the pyrolysis degree of the polymer matrix in aerodynamic environments, which significantly affects thermal protection and structural applications but is challenging to identify experimentally. Herein, non-equilibrium molecular dynamics simulations combined with experiments were implemented to determine the dependence of thermal conductivities on pyrolysis degree and microstructures for polymeric nanocomposites. We further explore the thermal transport mechanism through various contributions to the morphology. The results show that the thermal conductivity of the polymer matrix can be increased by a factor of 4.44 (from 0.27 W/m/K to 1.47 W/m/K) as the pyrolysis degree increases from 0 to 100%, and the thermal conductivity depends nonlinearly on the pyrolysis degree and temperature. Molecular dynamics simulations found that the side chains of the polymer matrix are rapidly scissored with the increasing pyrolysis degrees, and the structural ordering of the residual solids containing sp2 hybridization is enhanced, exhibiting graphene-like microtopological features, which reduces phonon scattering and makes thermal transport more efficient. This work provides insight into the linkage between the thermal transport properties and the pyrolysis degree of polymeric nanocomposites, which is valuable for improving the thermal transport performance and modeling ablation response for polymeric nanocomposites.

Effect of the advancing contact angle on coal particle–bubble detachment at the mesoscopic scale

Improving particle–bubble adhesion stability and reducing particle–bubble detachment are crucial in froth flotation processes. In this study, a nonequilibrium molecular dynamics approach is applied to investigate the dynamic detachment process of low-rank coal particles adhering with/without a collector from bubbles under different pull forces at the mesoscopic scale. Results indicate that a collector oil film on the surface of a low-rank coal particle can considerably increase the advancing contact angle and counteract the contraction of the three-phase contact line, enabling the particle to withstand a greater detachment force. Conversely, the residual water clusters on the particle surface tend to cause contraction of the contact line. Force analysis of the liquid surrounding the particle and the contact line indicates that the maximum static capillary force between the particle and the bubble is equivalent to the capillary force when the bending angle of the liquid surface reaches the advancing contact angle. When the bending angle of the liquid surface exceeds the forward contact angle, the force on the liquid around the contact line becomes unbalanced, causing the contact line to move to decrease the bending angle of the liquid surface and restore balance to the surrounding liquid. Hence, on mesoscopic scales of ten to tens of nanometers, the advancing contact angle still determines particle–bubble detachment. The study of particle–bubble detachment mechanisms at the mesoscopic scale is expected to establish a relation between microscopic mechanisms at the molecular scale and macroscopic experimental studies at the micrometer to millimeter scale

Effect of charge inversion on the electrokinetic transport of nanoconfined multivalent ionic solutions

Understanding the effects of phenomena occurring at electrically charged interfaces, such as charge inversion (CI), is crucial for enabling electroosmosis as an efficient transport mechanism in nanodevices. Here, we employ molecular dynamics (MD) simulations to systematically analyze the effect of CI on the electrokinetic transport of multivalent ionic solutions confined in amorphous silica nanochannels. We employ mixtures of monovalent and multivalent counterions while fixing the total ionic concentration to establish correlations between observed phenomena and the amount of multivalent ionic species in the electrolyte solution. The results show that the development of CI is related to a decrease in the mobility of the fluid layers adjacent to the charged surface. In addition, we observe that interfacial overcharging disrupts the water molecular orientation in the fluid layers adjacent to the channel walls. From the non-equilibrium MD simulations of electro-osmotic flow, we disclose the influence of phenomena related to the presence of CI. In particular, flow reversal occurs in scenarios involving CI due to increased local viscosity and a higher concentration of coions within the hydrodynamically mobile and electrokinetically active region of the charged interface. We also find that the magnitude of the wall zeta (ζ) potential displays a monotonic increase with the development of CI in the system. Moreover, we explain why positioning the wall ζ potential at an imaginary (slip) plane, which separates the hydrodynamically mobile and immobile fluid, is misleading.

Application of a thermal circuit model for the prediction of interfacial thermal resistance between water and a nanostructure surface using molecular dynamics simulations

The present study proposes various thermal circuit models (TCMs) to predict interfacial thermal resistances (ITRs) between water and a nanostructure surface. The effects of water models, nanopillar widths on nanostructure surfaces, and composite surfaces on the relationship between ITRs calculated by non-equilibrium molecular dynamics simulations and ITRs predicted by TCMs are investigated in water-copper and water-graphene-copper systems. The results reveal that the ITRs predicted by some TCMs are slightly affected by water pressure, while MD-calculated ITRs of the nanostructure surfaces in the Wenzel states agree with most of the predicted ITRs. Additionally, it is demonstrated that water molecules in the groove of a nanostructure surface play a crucial role in TCMs in the Cassie-Baxter states. Considering the effective contact region of water molecules in the groove in the Cassie-Baxter states, the use of TCMs can effectively reduce the prediction error of ITRs.

Effect of Electric Fields on the Decomposition of Phosphate Esters.

Phosphate esters decompose on metal surfaces and form protective polyphosphate films. For many applications, such as in lubricants for electric vehicles and wind turbines, an understanding of the effect of electric fields on molecular decomposition is urgently required. Experimental investigations have yielded contradictory results, with some suggesting that electric fields improve tribological performance, while others have reported the opposite effect. Here, we use nonequilibrium molecular dynamics (NEMD) simulations to study the decomposition of tri-n-butyl phosphate (TNBP) molecules nanoconfined between ferrous surfaces (iron and iron oxide) under electrostatic fields. The reactive force field (ReaxFF) method is used to model the effects of chemical bonding and molecular dissociation. We show that the charge transfer with the polarization current equalization (QTPIE) method gives more realistic behavior compared to the standard charge equilibration (QEq) method under applied electrostatic fields. The rate of TNBP decomposition via carbon-oxygen bond dissociation is faster in the nanoconfined systems than that in the bulk due to the catalytic action of the surfaces. In all cases, the application of an electric field accelerates TNBP decomposition. When electric fields are applied to the confined systems, the phosphate anions are pulled toward the surface with high electric potential, while the alkyl cations are pulled to the surface with lower potential, leading to asymmetric film growth. Analysis of the temperature- and electric field strength-dependent dissociation rate constants using the Arrhenius equation suggests that, on reactive iron surfaces, the increased reactivity under an applied electric field is driven mostly by an increase in the pre-exponential factor, which is linked to the number of molecule-surface collisions. Conversely, the accelerated decomposition of TNBP on iron oxide surfaces can be attributed to a reduction in the activation energy with increasing electric field strength. Single-molecule nudged-elastic band (NEB) calculations also show a linear reduction in the energy barrier for carbon-oxygen bond breaking with electric field strength, due to stabilization of the charged transition state. The simulation results are consistent with experimental observations of enhanced and asymmetric tribofilm growth under electrostatic fields.

Genetic algorithm optimization of langevin thermostat and thermal properties of graphene-aluminum nanocomposites: a molecular dynamics

Abstract The thermal properties of a laminated structure of graphene-coated aluminum composite nanomaterial were investigated through non-equilibrium molecular dynamics (NEMD) simulations to address the problem of temperature deviation in the thermostat volume applied. This paper presents a new insight into the best values of timestep and Langevin thermostat damping parameters for each atom in the nanomaterial with different size configurations using the genetic algorithm (GA) method by considering the timestep and thermostat damping parameters for each atom type, as well as the thickness of the nanomaterial, the thermostat, buffer, and heat flow lengths. The initial population results indicate that the thermostat temperature deviation increases with higher thermostat damping coefficients and timestep. However, the deviation decreases significantly with increased heat flow and thermostat lengths. Variations in buffer length and aluminum thickness do not have a significant effect on temperature. The application of a GA for optimization leads to a decrease in thermostat temperature deviation. The optimized parameters resulted in better thermostat temperature deviations when analyzing the temperature, aluminum thickness, and both buffer and thermostat lengths. Additionally, the thermal conductivity of aluminum-graphene nanomaterial decreases with increasing temperature, buffer length, and aluminum thickness, but increases by up to 9.85% with increasing thermostat length.

Correlation effect between the capillary flow and evaporation rate of water in a 2-D photothermal membrane

To shed light on the correlation effect between the capillary flow and phase transition of water in photothermal evaporation, the non-equilibrium molecular dynamics (NEMD) simulations were employed to investigate the evaporation of water in Multilayered graphene oxide membranes. It was initially found that water molecules need to overcome strong energy barriers to migrate upward by layers before evaporation, with the energy barrier increasing with the reduction in nanogap, hydroxyl content, or increase in layer spacing of graphene sheets. The strong energy barrier accordingly results in dramatic sensible heat conversion and high evaporation temperature. It was shown that the enhanced temperature and hydrophilicity both weaken the hydrogen bond strength of water molecules confined in 2-dimensional channels, contributing to the reduced evaporation enthalpy. The preferential pathways in graphene channels or the weakened hydrogen bonds in hydroxyl-functionalized graphene of water molecules were found to be important influences on capillary flow energy consumption, leading to differences in evaporation energy supply. In combination with these effects, the evaporation rate could be increased by up to 2 folds and an interfacial enthalpy of evaporation below 1500 kJ kg−1 could be achieved by properly tuning the configuration and chemical properties of the graphene membrane. Our findings lay a theoretical foundation for providing a strategy to improve the interfacial evaporation performance in desalination.

Mechanisms determining water permeability and ion selectivity of multilayered glycine-functionalized graphene oxide membranes in saline water reverse osmosis processes: A computational investigation

Non-Equilibrium Molecular Dynamics simulations were performed to examine in detail the mechanisms involved in a reverse osmosis process for the removal of monovalent (Na+, Cl−) and divalent (Mg2+) ions from saline water, using multilayered glycine-functionalized graphene oxide (GO) membranes. We varied parameters such as the degree of functionalization of the GO flakes and the structural flexibility of the membranes and examined their performance in a wide range of applied pressure gradients. In all models examined, water permeability was found to be almost 2 orders of magnitude larger compared to that in conventional reverse osmosis (RO) membranes. High structural flexibility of the membranes was found to increase water permeability due to the formation of additional nanochannels that favor faster water transport. Functionalization of the GO flakes by the glycinate groups resulted in a decrease in water permeability due to the formation of a denser hydrogen-bonding network within the membranes. Three main mechanisms were identified regarding ion selectivity. Due to the overall negative charge of the membranes, the strong adsorption of the Mg2+ ions onto the GO flakes resulted in their full rejection. Weaker adsorption of the Na+ ions led only to their partial rejection. At the same time, the electrostatic repulsion of the Cl− ions by the membranes was found to result in almost full removal levels at high pressure gradients. The rejection of the monovalent ions was found to decrease with the structural flexibility of the membranes and increase with the functionalization of the GO flakes, in response to the changes observed in water permeability on these membrane modifications. The pressure-dependent deformability of the structurally flexible membranes combined with charge accumulation close to the membrane's entrance, resulted in a non-monotonic dependence of the monovalent ion rejection on the applied pressure difference.

Investigation of the effect of salts on dye wastewater separation performance in TpPa-1 membrane via molecular dynamics simulation

The rapid progress in science and technology increases the demand for dyes and results in the pollution of dye wastewater. Recently, covalent organic frameworks (COFs) have been suggested as a promising candidate for wastewater treatment. The influence of inorganic salt ions on the transport and separation processes of dye wastewater is poorly understood. In this work, the transport and separation mechanisms of dyes are investigated via nonequilibrium molecular dynamics simulations (NEMD). It is found that the inorganic salts enhance the rejection rate of different dyes in TpPa-1 membranes. The fundamental reason for this phenomenon is the hindrance caused by inorganic salt ions, which change the orientation of dye molecules as they pass through the pore channels of the TpPa-1 membrane. As a result, the dye molecules show a greater tendency to form π–π interaction with the TpPa-1 membrane surface, which accelerates their adsorption.

Thermal Transport and Thermal Polarization of Water in the Supercooled Regime.

The behavior of water in the deep supercooled regime has attracted significant interest, motivated by the hypothesis of the second critical point of water. Previous studies indicated the existence of a water anomaly, characterized by a minimum in the thermal conductivity of water. Here, we employ nonequilibrium molecular dynamics computer simulation and the TIP4P/2005 water force field to investigate the thermal conductivity of supercooled water targeting four different isobars, 1, 200, 700, and 1200 bar. We demonstrate using NEMD simulations the existence of minima in thermal conductivity associated with the maximum isothermal compressibility and the minimum speed of sound in water, hence establishing a firm connection with the second critical point of liquid water. Moreover, we demonstrate that thermal gradients polarize supercooled water with a thermal polarization coefficient of several mV/K. We explain the thermal polarization effect using a theoretical formulation introduced recently that connects the thermal polarization effect to the isobaric thermal expansion coefficient.

Kinetics of the bcc to hcp phase transition in [formula omitted]-oriented single crystal iron under ramp compression

Non-Equilibrium Molecular Dynamics simulations have been used to investigate the kinetics of the bcc to hcp phase transition in ramp-compressed [001]-oriented single crystal iron in the presence of defects consisting in micro-voids. Depending on the ramp maximum pressure, the microstructure of the daughter hcp phase is observed to be strongly affected by the plastic activity in the parent bcc phase. The hcp phase is first nucleated around the pre-existing microvoids, then it grows along preferential crystal directions and the subsequent volume reduction induces a local relaxation of the pressure. The velocity of the bcc-hcp phase boundary is shown to increase gradually with maximum pressure up to supersonic values. For high over-compression beyond the transition onset pressure, the phase transformation becomes almost instantaneous and essentially independent of the initial defects. Meanwhile, the ramp wave steepens with propagation distance and evolves into a shock front, up to the over-driven regime classically observed in both experiments and hydrodynamic, macro-scale simulations.

Tailoring thermal transport and confinement effect of GaN/Si3N4 nanowires utilizing core–shell architecture

The thermal transport properties of nanowires (NWs) can be significantly influenced by the implementation of a core-shell structure, which introduces interface scattering and phonon localization effects, opening avenues for novel applications. In this paper, we use the method of non-equilibrium molecular dynamics to simulate the effects of system temperature, cross-sectional width, and nanopillar interface on the thermal transport of GaN/Si3N4core-shell NWs. The thermal transport process of phonons in core-shell NWs is studied by calculating the vibrational density of states, phonon participation rate, and dispersion curve. The results show that the core-shell NWs characterized by smooth interfaces exhibit a 17.4% decrease in thermal conductivity (TC) at room temperature when contrasted with pristine GaN NWs. Furthermore, the TC of GaN/Si3N4core-shell NWs can be further reduced by adding nanopillars at the interface. Due to resonance effect, thus effectively regulating the thermal transport. The presence of nanopillars increases phonon-surface scattering intensity at low-frequency and modifies phonon dispersion to decrease the group velocity. In addition, the hybridization of phonon modes between those of the nanopillars and the Si3N4shell gives rise to numerous dispersionless resonance phonon modes that span the entire phonon spectrum. This research delves into the effects of nanopillars and interfaces on thermal transport, providing important guidance for understanding confinement effects and establishing a robust theoretical basis for the regulation of thermal transport.

Viscosity of Asphalt Binder through Equilibrium and Non-Equilibrium Molecular Dynamics Simulations

Viscosity is a curial indicator for evaluating asphalt performance, representing its ability to resist deformation under external forces. The Green–Kubo integral in equilibrium molecular dynamics simulations and the Muller-Plathe algorithm in reverse non-equilibrium molecular dynamics simulations were used to calculate the asphalt viscosity. Meanwhile, the key parameters of both methods were rationalized. The results show that in equilibrium calculations, using a 1/t weighting for the viscosity integral curve results in a well-fitted curve that closely matches the original data. The isotropy of the asphalt model improves for atomic counts exceeding 260,000, rendering viscosity calculations more reasonable. When the viscosity did not converge, it increased linearly with the number of atoms. In non-equilibrium calculations, the number of region divisions had almost no effect on the viscosity value. A momentum exchange period of 20 timesteps exhibits a favorable linear trend in velocity gradients, and an ideal momentum exchange period was found to be between 10 and 20 timesteps. As the model size increased, the linear relationship with the shear rate became more pronounced, and the isotropy of the asphalt system improved. Using an orthogonal simulation box with a side length of 75 Å effectively meets the computational requirements.

Chemical short-range order increases the phonon heat conductivity in a refractory high-entropy alloy

We study the effects of the chemical short-range order (SRO) on the thermal conductivity of the refractory high-entropy alloy HfNbTaTiZr using atomistic simulation. Samples with different degrees of chemical SRO are prepared by a Monte Carlo scheme. With increasing SRO, a tendency of forming HfTi and TiZr clusters is found. The phonon density of states is determined from the velocity auto-correlation function and chemical SRO modifies the high-frequency part of the phonon density of states. Lattice heat conductivity is calculated by non-equilibrium molecular dynamics simulations. The heat conductivity of the random alloy is lower than that of the segregated binary alloys. Phonon scattering by SRO precipitates might be expected to reduce scattering times and, therefore, decrease thermal conductivity. We find that, in contrast, due to the increase of the conductivity alongside SRO cluster percolation pathways, SRO increases the lattice heat conductivity by around 12 %. This is expected to be a general result, extending to other HEAs.

Nonequilibrium dynamics of coupled oscillators under the shear-velocity boundary condition.

Deterministic and stochastic coupled oscillators with inertia are studied on the rectangular lattice under the shear-velocity boundary condition. Our coupled oscillator model exhibits various nontrivial phenomena and there are various relationships with wide research areas such as the coupled limit-cycle oscillators, the dislocation theory, a block-spring model of earthquakes, and the nonequilibrium molecular dynamics. We show numerically several unique nonequilibrium properties of the coupled oscillators. We find that the spatial profiles of the average value and variance of the velocity become nonuniform when the dissipation rate is large. The probability distribution of the velocity sometimes deviates from the Gaussian distribution. The time evolution of kinetic energy becomes intermittent when the shear rate is small and the temperature is small but not zero. The intermittent jumps of the kinetic energy cause a long tail in the velocity distribution.

CO2 and NO2 formation on amorphous solid water

Context. The dynamics of molecule formation, relaxation, diffusion, and desorption on amorphous solid water (ASW) is studied in a quantitative fashion. Aims. The formation probability, stabilization, energy relaxation, and diffusion dynamics of CO2 and NO2 on cold ASW following atom+diatom recombination reactions are characterized quantitatively. Methods. Accurate machine-learned energy functions combined with fluctuating charge models were used to investigate the diffusion, interactions, and recombination dynamics of atomic oxygen with CO and NO on ASW. Energy relaxation to the ASW and into water internal degrees of freedom were determined from the analysis of the vibrational density of states. The surface diffusion and desorption energetics were investigated with extended and nonequilibrium MD simulations. Results. The reaction probability is determined quantitatively and it is demonstrated that surface diffusion of the reactants on the nanosecond time scale leads to recombination for initial separations of up to 20 Å. After recombination, both CO2 and NO2 stabilize by energy transfer to water internal and surface phonon modes on the picosecond timescale. The average diffusion barriers and desorption energies agree with those reported from experiments, which validates the energy functions. After recombination, the triatomic products diffuse easily, which contrasts with the equilibrium situation, in which both CO2 and NO2 are stationary on the multinanosecond timescale.

Mechanism of interfacial thermal resistance variation in diamond/Cu/CNT tri-layer during thermal cycles

The multilayer materials in aerospace applications generally subject to thermal shock, which may significantly affect the interfacial thermal resistance (ITR) and cause serious overheating issues. Therefore, understanding the interfacial thermal transport under temperature shock is essential for effective thermal management in aerospace devices. In the present study, the variation in ITR within the diamond/Cu/carbon nanotube tri-layer after thermal cycles was investigated through non-equilibrium molecular dynamics (NEMD) simulations. We conducted four groups of investigations, exploring the impact of maximum cycle temperature, heating/cooling rates, number of cycles, and interfacial structure on ITR. The results demonstrate that while the matching degree of phonon density of states remains almost constant, the ITR varies significantly among all groups. The structural deformations and changes in lattice type can be observed in the Cu layer. Based on these findings, we proposed an “atomic distribution method” to elucidate the mechanism behind ITR variation, which was verified as applicable and precise. Additionally, the thermal rectification results demonstrate the significant effect of interfacial structure on the rectification coefficient, indicating that interfacial transport and lattice thermal conduction deserve more in-depth study in the future. This work provides a novel perspective on understanding thermal transport across the irregular and complex interfaces, which is significant for the thermal management in aerospace field.

Orientation of Linear Liquid Alkane on Solid Surface of Face Centered Cubic Lattice of (100), (110) and (111)

Solid-liquid (S-L) interfaces are found in numerous engineering applications, such as lubrication and coating systems as well as thermal interface materials. Understanding the interactions between (S-L) is crucial for optimizing various engineering applications. The main objective of this study is to provide insight regarding liquid adsorption on solid surfaces using non-equilibrium molecular dynamics simulations. To achieve this goal, a liquid confined between solid surfaces will be modeled to match the real state of contact interfaces using a constant temperature as a baseline. The results highlight a significant relationship between the peak height values of the liquid adsorption layer, density profile, and radius of gyration. Specifically, the peak height density at the S-L interfaces for the crystal plane (110) is 784.756 kg/m3, followed by (100) at 801.786 kg/m3, and finally the highest is for the crystal plane (111), at 966.940 kg/m3. Whereas the radius of gyration at the S-L interfaces for crystal planes (100) and (111) is approximately 7.45 × 10-21 m2, but for crystal plane (110) it is less and measures approximately 7.06 × 10-21 m2. Conclusion, the adsorption layer of solid density near solid-liquid interfaces is significantly influenced by the peak height of the solid's density. Higher density results in a higher adsorption layer of liquid near solid-liquid interfaces. The number of solid density layers does not affect the height of adsorption layers for liquids.

A quantitative theory and atomistic simulation study on the soft-sphere crystal-melt interfacial properties. I. Kinetic coefficients.

By employing non-equilibrium molecular dynamics (NEMD) simulations and time-dependent Ginzburg-Landau (TDGL) theory for solidification kinetics [Cryst. Growth Des. 20, 7862 (2020)], we predict the kinetic coefficients of FCC(100) crystal-melt interface (CMI) of soft-spheres modeled with an inverse-sixth-power repulsive potential. The collective dynamics of the local interfacial liquid phase at the equilibrium FCC(100) CMIs are calculated based on a recently proposed algorithm [J. Chem. Phys. 157, 084 709 (2022)] and are employed as the resulting parameter that eliminates the discrepancy between the predictions of the kinetic coefficient using the NEMD simulations and the TDGL solidification theory. A speedup of the two modes of the interfacial liquid collective dynamics (at wavenumbers equal to the principal and the secondary reciprocal lattice vector of the grown crystal) is observed. With the insights provided by the quantitative predictive theory, the variation of the solidification kinetic coefficient along the crystal-melt coexistence boundary is discussed. The combined methodology (simulation and theory) presented in this study could be further applied to investigate the role of the inter-atomic potential (e.g., softness parameter s = 1/n of the inverse-power repulsive potential) in the kinetic coefficient.