Molecular Dynamics Study of Anti-Wear Erosion and Corrosion Protection of PTFE/Al2O3 (010) Coating Composite in Water Hydraulic Valves

Cavitation erosion and corrosion commonly occur on the surface of fluid dynamic system components, mostly water hydraulic valves, causing the failure of metal parts. Coating of polytetrafluoroethylene (PTFE) on Al2O3 (010) was created by varying the chain length of polytetrafluoroethylene. Calculations were conducted by molecular dynamic (MD) simulations. This study shows that the K10 and K20 chain lengths’ mechanical properties possess negative elastic, shear, and bulk modulus values. We have found that the K10 chain length composition shows the high results of binding energy and negative bulk modulus of 6267.16 kJ/mol and −3709.54 GPa, respectively. The K10 chain length was observed to possess a higher cohesive energy density (CED) and solubility parameter of (6.885 ± 0.00076) × 109 J/m3 and (82.974 ± 0.005) (J/cm3)0.5, respectively. It was also found that increasing the chain length contributes to decreasing the binding energy and solubility parameter of PTFE/Al2O3 (010) composition. These results are vital for overcoming the repetitive regime of high compressive strength of water microjets on the valves’ material surface. Improved values of the cohesive energy density and solubility parameters imply the water’s superior hydrophobic effect.


Introduction
Polymeric materials are prevalent in nature and daily life. It has been demonstrated that composite materials with polymer matrices cover surfaces pressured with cavitation flows and abrasive pattern suspension [1]. Polymers are most widely used in applications and structures that work under such hydrodynamic fluid conditions due to the characteristics of mechanical resilience, chemical degradation, and abrasive erosion resistance. Protecting the surfaces of hydro mechanical parts that work under cavitation wear conditions, such as rotors and valves in hydraulic machines, is among such applications [1]. Molecular dynamic simulations have contributed a lot in recent years to amorphous alloy development and improvement. Several investigations have been documented into amorphous alloy forming capacity and mechanical properties using molecular dynamic simulations [2]. Xu et al. [3] uses Monte Carlo dynamic simulations to investigate the chain length effects in polymer compound systems of specific chain length on stereocomplex crystallites (SCs). Studies have shown that decreasing the chain's length increases the concluding stereocomplex crystallites material. In the short-chain systems, three aspects lead to improved SC creation, comprising miscibility between dissimilar polymer types, chain flexibility, and nuclear mode. The structures with shorter chains are more miscible, have better compartment mobility, and tend to crystallize in a parallel intermolecular alignment; thus, providing a more significant SC material than long chains.
Meanwhile, in ref. [4], Chen et al. conducted a study investigating the influence of glucose polymer chain length on the heat and physical characteristics of milk-protein isolate (MPI) sugar-nutrient drinks containing 8.5% w/w total protein and 5% w/w carbohydrate. The findings show that the chain length of the glucose polymer heavily influences MPI carbohydrate supplement beverages. Zhao et al. [5] also investigated chain-length and temperature-based tensile and shear failure activity using molecular dynamic simulations in amorphous polymers. They concluded that the breakup intensity under shear and tension rises with chain length and the temperature increases. The breakup rates under shear are often higher for a defined chain length and temperature than those below the same high strain. The tensile and shear stresses decrease with rising temperature for a given chain length. The effects of polymer chain length on formless drug recrystallization hindrance were investigated by Pacult et al. [6].
To contemplate this issue, they arranged a drug-polymer mix in the proportion of 10:1 (drug to polymer) comprising bicalutamide (BIC) and polyvinylpyrrolidone (PVP) with peculiar chain lengths K10, K30, K90. They believe that the polymer chain may have different configurations in space that contribute to better or worse stability in research, based on its length. Mlela et al. [7] aimed to match Cr 2 O 3 , Al 2 O 3 , and Ti 2 O 3 with the most appropriate polymer coating. Through using the Preference Ranking Organization Method for Enrichment Evaluations (PROMETHEE), the most robust obtained polymer in a watery environment valve with an evaluated efficiency attribute is polytetrafluoroethylene with a greater outranking (0.5932052). Their molecular dynamics studied the robust bonding concerning the superior cleave plane between polytetrafluoroethylene and the selected substrate. PTFE/Al 2 O 3 cleaved in (010) crystallographic plane was found to have the strongest bond in terms of binding energy (3188 kJ/mol).
In this paper, we carried out molecular dynamic (MD) simulations of PTFE/Al 2 O 3 (010) by varying the PTFE chain length to calculate the binding energy, mechanical properties, cohesive energy density (CED), and solubility parameter. The greater the mechanical properties' parameters of PTFE/Al 2 O 3 (010), the greater the potential to resist stress without microjet deformation and breakage. Binding energy and coherent energy density (CED) reflect their component molecules' bonding strength. The solubility parameter holds great promise for developing a poorly water-soluble composition and, hence, smaller corrosion vulnerability. This work's results significantly impact the protection of water hydraulics' valves from erosion and corrosion and thus improve the water hydraulic system's performance.

Molecular Dynamics (MD) Simulation and Properties Calculation
By utilizing Material Studio ® software version 7.0 created by Accelrys Inc. (San Diego, CA, USA), molecular dynamic simulations were applied to examine the binding energies, mechanical and dissolvability properties of the various PTFE/Al 2 O 3 (010). Figure 1 shows the simulation and calculation models of different phases of this work. In this research, as seen in Figure 2a,b, tasks were to build the cleaved crystal surface of Al 2 O 3 (010) of a thickness of 5.4539 nm, with the lattice parameters of U = 1.2991 nm, V = 0.4759 nm, and θ = 90 • . The next task was to optimize the geometry using the Forcite module obtained in materials studio (MS) version 7.0 of Accelrys Inc. and a smart algorithm [8] with the maximum iteration of 1000 steps, with ultra-fine quality that employs 8.368 × 10 −5 kJ/mol of convergence energy and a force tolerance of 0.04184 kJ/mol/nm with the displacement of 1.0 × 10 −6 nm. The surface of the crystals of Al 2 O 3 (010) was minimized further by the discover module using a smarter ultra-fine minimizer with the maximum iteration of 2000 to achieve the stable structure of the lowest energy.   The Al2O3 (010) surface area was increased to obtain the 3D periodicity of the lattice vector of U × V= 6 × 6 by keeping a vacuum thickness of 0.0 nm. The 3D triclinic crystal lattice had the parameters of a = 7.7946 nm, b = 2.8554 nm, c = 2.8554 nm, and α = β = γ = 90° was formed, see Figure 2c. The monomer repeat unit, tetrafluoroethylene, was imported from the materials library in the materials studio (MS) version 7.0 of Accelrys Inc. to build a polymer and amorphous cell. As usual, the job optimized by Forcite uses ultra-fine quality and a smart algorithm with the 1.0 × 10 3 maximum number of geometry optimization cycles (maximum iteration). By specifying the temperature 298 K, 2.20 g/cm 3 targeted density, the periodic amorphous cells of the final configuration were created in the amorphous cell tool, as Figure 3 shows. Table 1 gives the designated chain number, chain length, and polymer configuration of a different model. We must bring together the Al2O3 (010) surface layer and the other polytetrafluoroethylene layer for interaction. We specified a vacuum of 7.0 nm spacing for polytetrafluoroethylene in c-direction to secure enough space for physicochemical phenomena. The molecular dynamics simulation is essential for obtaining the actual density and conformation of wear-resistant coating to reach the equilibration state. We established the simulation at room temperature (298 K), with the Andersen temperature control system with a collision frequency of every 325 steps. The velocity Verlet integrator algorithm was employed considering the canonical ensemble of the constant number of particles, constant volume, and constant temperature (i.e., NVT)    The Al2O3 (010) surface area was increased to obtain the 3D periodicity of the lattice vector of U × V= 6 × 6 by keeping a vacuum thickness of 0.0 nm. The 3D triclinic crystal lattice had the parameters of a = 7.7946 nm, b = 2.8554 nm, c = 2.8554 nm, and α = β = γ = 90° was formed, see Figure 2c. The monomer repeat unit, tetrafluoroethylene, was imported from the materials library in the materials studio (MS) version 7.0 of Accelrys Inc. to build a polymer and amorphous cell. As usual, the job optimized by Forcite uses ultra-fine quality and a smart algorithm with the 1.0 × 10 3 maximum number of geometry optimization cycles (maximum iteration). By specifying the temperature 298 K, 2.20 g/cm 3 targeted density, the periodic amorphous cells of the final configuration were created in the amorphous cell tool, as Figure 3 shows. Table 1 gives the designated chain number, chain length, and polymer configuration of a different model. We must bring together the Al2O3 (010) surface layer and the other polytetrafluoroethylene layer for interaction. We specified a vacuum of 7.0 nm spacing for polytetrafluoroethylene in c-direction to secure enough space for physicochemical phenomena. The molecular dynamics simulation is essential for obtaining the actual density and conformation of wear-resistant coating to reach the equilibration state. We established the simulation at room temperature (298 K), with the Andersen temperature control system with a collision frequency of every 325 steps. The velocity Verlet integrator algorithm was employed considering the canonical ensemble of the constant number of particles, constant volume, and constant temperature (i.e., NVT) The Al 2 O 3 (010) surface area was increased to obtain the 3D periodicity of the lattice vector of U × V = 6 × 6 by keeping a vacuum thickness of 0.0 nm. The 3D triclinic crystal lattice had the parameters of a = 7.7946 nm, b = 2.8554 nm, c = 2.8554 nm, and α = β = γ = 90 • was formed, see Figure 2c. The monomer repeat unit, tetrafluoroethylene, was imported from the materials library in the materials studio (MS) version 7.0 of Accelrys Inc. to build a polymer and amorphous cell. As usual, the job optimized by Forcite uses ultra-fine quality and a smart algorithm with the 1.0 × 10 3 maximum number of geometry optimization cycles (maximum iteration). By specifying the temperature 298 K, 2.20 g/cm 3 targeted density, the periodic amorphous cells of the final configuration were created in the amorphous cell tool, as Figure 3 shows. Table 1 gives the designated chain number, chain length, and polymer configuration of a different model. We must bring together the Al 2 O 3 (010) surface layer and the other polytetrafluoroethylene layer for interaction. We specified a vacuum of 7.0 nm spacing for polytetrafluoroethylene in c-direction to secure enough space for physicochemical phenomena. The molecular dynamics simulation is essential for obtaining the actual density and conformation of wear-resistant coating to reach the equilibration state. We established the simulation at room temperature (298 K), with the Andersen temperature control system with a collision frequency of every 325 steps. The velocity Verlet integrator algorithm was employed considering the canonical ensemble of the constant number of particles, constant volume, and constant temperature (i.e., NVT) Coatings 2020, 10, x; doi: FOR PEER REVIEW 4 of 14  The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) ab initio forcefield was used for all these applications. The simulation was performed with an interval of 1 femtosecond (fs) in each scenario, and the simulation steps were 2.0 × 10 5 with a cumulative dynamic time of 200.0 picoseconds (ps). In the meantime, atom configurations were saved as trajectory data in every 5000 frames in each case and dynamic. Figures 4-7 show the coating configuration of PTFE's amorphous cell on Al2O3 (010) before and after MD simulation. Equation (1) was designed to model the interaction of covalent bonds between atoms and the interaction of van der Waals (vdW) in the wear-resistant layers [9]. Additionally, Equation (1) describes all bonded and non-bonded interactions in twelve forms. The first ten terms in Equation (1) are used to model cross-linked diagonals and off-diagonals, and the last two terms model non-bonded interactions (electrostatic (coulombic) and van der Waals). The first four terms reflect, in particular, the energies produced by the bonds stretching (b), bonds bending (θ), bond torsion (φ), and Wilson out-of-plane (χ) internal coordinates. The cross-coupling terms include combinations of two or three internal  The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) ab initio forcefield was used for all these applications. The simulation was performed with an interval of 1 femtosecond (fs) in each scenario, and the simulation steps were 2.0 × 10 5 with a cumulative dynamic time of 200.0 picoseconds (ps). In the meantime, atom configurations were saved as trajectory data in every 5000 frames in each case and dynamic. Figures 4-7 show the coating configuration of PTFE's amorphous cell on Al 2 O 3 (010) before and after MD simulation. Equation (1) was designed to model the interaction of covalent bonds between atoms and the interaction of van der Waals (vdW) in the wear-resistant layers [9]. Additionally, Equation (1) describes all bonded and non-bonded interactions in twelve forms. The first ten terms in Equation (1) are used to model cross-linked diagonals and off-diagonals, and the last two terms model non-bonded interactions (electrostatic (coulombic) and van der Waals). The first four terms reflect, in particular, the energies produced by the bonds stretching (b), bonds bending (θ), bond torsion (ϕ), and Wilson out-of-plane (χ) internal coordinates. The cross-coupling terms include combinations of two or three internal     The static mechanical properties of the equilibrated system of PTFE/Al 2 O 3 (010) need to be calculated. The balance of temperature and energy (non-bond and potential) simultaneously determines the system's equilibrium, significantly, when energy and temperature fluctuate within a maximum of ±5%-10% [10]. The scenario is considered to reach the equilibrium state when the temperature is 298 K. The Forcite module in materials studio (MS) version 7.0 of Accelrys Inc. was used to implement all molecular dynamic calculations. Mechanical properties, cohesive energy density, and binding energy were calculated using the universal force field (UFF) based on the quantum mechanical method. The last eleven frames of the stable systems were determined. The stiffness matrix and flexibility matrix of the composite, shear modulus, and bulk modulus were determined using Reuss, Voigt, and Hill approximation [11]. The macro-level enhanced anti-erosive wear is close to the isotropic material, and comparable to isotropic materials.
The relationship between stress (σ) and strain (ε), with stiffness or elastic constants (C ij ) of the material, follows the generalized Hooke's law [12,13], as shown in Equation (2). Lame constants, lambda (λ), and Mu (µ) describe the stress-strain relationships of isotropic material, as demonstrated in Equation (3). Meanwhile, Equation (4) is used to evaluate standard elastic parameters such as elastic modulus (E), bulk modulus (B), shear modulus (G), and Poisson ratio (v). Further, the binding energy (E b ) is the mechanical (binding) energy required to disassemble a whole composite into separate parts. Equation (5) is used to quantify the binding energy where E t is the single-point total energy of the equilibrated system after the interaction, E p is the single point energy of the polymer, and E s is the single point energy of the Al 2 O 3 (010) surface before interaction with the polymer. The Hildebrand solubility parameter and cohesive energy density are defined in Equation (6) and calculated through the Forcite tool in each erosive wear-resistant composite of PTFE/Al 2 O 3 (010).

Results and Discussion
The results obtained in the MD simulations and the elastic constants of each PTFE/Al 2 O 3 (010) are shown in Tables 2 and 3 below. In the case of the isotropic materials, Poisson's ratio ν may be expressed in terms of the bulk modulus B and the shear modulus G, which relate, respectively, to changes in size and shape. The Poisson ratio of all stable isotropic materials is found between numerical limits 1 2 and −1 [14]. If the material were not isotropic, the Poisson ratio would be greater than one (see Figure 8). The plastic range is typically between 0.2 and 0.4; natural rubber is about 0.5 and thermoplastic starch (TPs) is between 0.1 and 0.4 [15]. Nevertheless, according to the polymer properties database in [16], the rubber Poisson ratio is never exactly 0.5; this means that the elastic bulk modulus would reach infinity, which is not possible. However, very close to 1 2 values have been reported (0.4999) [16]. The composite attributed to K10 chain lengths of PTFE is typically a natural rubber, as seen in Figure 9e and Table 3. In comparison, the composite assigned to the chain length of K20 belongs to plastic or thermoplastic starch. The chain length of K30 shows the unstable isotropic, and the K40 chain length seems to possess the anisotropic range of Poisson's ratio, unstable by volume change but stable if constrained. The idea that the Milton map of bulk isothermal Coatings 2020, 10, 1214 7 of 12 modulus B versus shear module G, when B/G >> 1 and ν→ 1 /2 (vertical axis) means materials such as rubber are extremely incompressible. This idea was discussed and clarified in [14] by Greaves et al. Bulk module's negative compressibility is usually forbidden in classical thermodynamics, but it could be possible due to a recent study [17]. A stable negative bulk modulus of material has been observed experimentally in the constrained form [18]. Structures can consist of bi-material components that can have negative properties, significantly negative compressibility (negative bulk modulus, i.e., increased size when the external pressure is increased and decreased size when external pressure diminishes) [17]. The degree of compression typically represents bond strengths and, thus, softer materials [19]. For this reason, the bulk modulus of the K10 chain length of PTFE/Al 2 O 3 (010) in Figure 9c exhibits the negative value, which demonstrates the resistance of the PTFE/Al 2 O 3 (010) composite to erosion and corrosion compared to other simulated arrangements of the PTFE/Al 2 O 3 (010) composite. This effect is strongly supported by the obtained calculated binding (mechanical) energy and bulk modulus of 6267.16 kJ/mol and −3709.54 GPa, respectively (see Figure 9c,d). The mass density, bulk modulus, and Young's modulus work in the same way to ensure their auto-correlation functionality is similar [20]. Mass density can also be related to the corrosion degree of the material in the sense that, if the surfaces corrode, oxidizing occurs and eventually creates rust, which takes up more room than the initial material. This expansion and reduced density can lead to a variety of defects, such as cracking [20]. The higher the density of the material, the closer the particles are squeezed into the material, and the denser the PTFE/Al 2 O 3 (010) composite is, which contributes to hardness [21]. Negative effective shear modulus depicts axisymmetrically deformed composites with opposite axisymmetric loading conditions [22].

Stable
Stiff    It seems in Table 3 and Figure 9b that the K10 and K20 chain lengths of PTFE/Al2O3 (010) show the negative shear modulus of −25.56 and −97.99 GPa, respectively. Combining their binding energies of 6267.16 and 4271.05 kJ/mol, respectively, forms the most desired strengths to overcome the surge imploding of pressure from the water cavity, causing water valves' erosion and corrosion. The estimated implosive pressure was of the order of 10 3 -10 4 MPa and its period 2-3 μs [23]. The simulated composite of K10 chain length exhibits stronger incompressibility of −3709.54 GPa and abundant binding energy of 6267.16 kJ/mol and unique moderate negative elastic and shear modulus (Figure 9a,b). Negative stiffness is possible in the systems, such as prestrained objects, including post-buckled elements, single-cell models, and post-buckled flexible tubes, containing stored energy [24,25]. The negative stiffness principle has been studied theoretically and experimentally in composite materials [25]. Stable negative structural stiffness is known, both theoretically and experimentally. The negative material stiffness, namely the bulk modulus (inverse compressibility) of a constrained solid elastic body of any dimension, is stable within the elasticity principle [18]. Double-negative metamaterials can be obtained by combining negative effective mass and negative effective modulus [26]. Similar to electromagnetic (EM) metamaterials, the integration of resonant structures into blocks of an elastic metamaterial achieves negative successful elastic parameters [27]. Negative parameters, especially compressibility, will resist the compressive force generated by the microjets and prevent the material surface from undergoing plastic deformation (plastic flow), 10 20 It seems in Table 3 and Figure 9b that the K10 and K20 chain lengths of PTFE/Al 2 O 3 (010) show the negative shear modulus of −25.56 and −97.99 GPa, respectively. Combining their binding energies of 6267.16 and 4271.05 kJ/mol, respectively, forms the most desired strengths to overcome the surge imploding of pressure from the water cavity, causing water valves' erosion and corrosion. The estimated implosive pressure was of the order of 10 3 -10 4 MPa and its period 2-3 µs [23]. The simulated composite of K10 chain length exhibits stronger incompressibility of −3709.54 GPa and abundant binding energy of 6267.16 kJ/mol and unique moderate negative elastic and shear modulus (Figure 9a,b). Negative stiffness is possible in the systems, such as prestrained objects, including post-buckled elements, single-cell models, and post-buckled flexible tubes, containing stored energy [24,25]. The negative stiffness principle has been studied theoretically and experimentally in composite materials [25]. Stable negative structural stiffness is known, both theoretically and experimentally. The negative material stiffness, namely the bulk modulus (inverse compressibility) of a constrained solid elastic body of any dimension, is stable within the elasticity principle [18]. Double-negative metamaterials can be obtained by combining negative effective mass and negative effective modulus [26]. Similar to electromagnetic (EM) metamaterials, the integration of resonant structures into blocks of an elastic metamaterial achieves negative successful elastic parameters [27]. Negative parameters, especially compressibility, will resist the compressive force generated by the microjets and prevent the material surface from undergoing plastic deformation (plastic flow), ultimately enhancing erosion and corrosion protective performance. According to Lee SH and Lee SB in [28], the Hildebrand solubility parameter (δ H ) is usually expressed as the square root for a chemical compound's cohesive energy density (CED), as Equation (6) depicts.
where ∆U, ∆H vap and V are internal molar energy, vaporization enthalpy (298 K), and molar volume, respectively. The solubility parameter (δ H ) has been widely used for predicting the solubility of various chemicals in organic solvents. Water as a polar solvent is likely to cause a hydrophobic effect in a valve chamber. The solubility parameter as a square root of cohesive energy density is a notable characteristic that distinguishes the nonpolar compound's hydrophobic nature. The higher value of the solubility or cohesive energy density of the nonpolar compound, the higher the hydrophobicity and bond strength, and the higher the corrosion protective performance. Molecules with stronger interactions have higher values of the cohesive energy density and solubility parameters, promoting the hydrophobic effect and improving the valve chamber's anti-corrosion and wear properties. In this work, it has been identified that solubility parameter (δ H ) and cohesive energy density values tended to decrease as the chain length increased (see Table 4). Based on Table 4, the wear-resistant coating of PTFE/Al 2 O 3 (010) with K10 chain length shows promising results for valve protection. These results can resist water hydraulic cavitation erosion and corrosion due to its higher cohesive energy density (CED) and solubility parameter of (6.885 ± 0.00076) × 10 9 J/m 3 and (82.974 ± 0.005) (J/cm 3 ) 0.5 , respectively, compared to other simulated results of K (20, 30 and 40) of PTFE/Al 2 O 3 (010).

Conclusions
The mechanical properties of the K10 and K20 chain lengths of PTFE/Al 2 O 3 (010) possess negative values of elastic modulus, shear modulus, and bulk modulus. In most elastic structures, we expect positive stiffness, such that the deformed material experiences force in the same direction as the deformation. Negative stiffness is possible in structures, such as prestrained objects and post-buckled components, that contain stored energy essential to rebound cavitation jet impacts. The same applies to bulk modulus; negative bulk modules are shown to be feasible in a prestrained lattice and other crystalline materials. Furthermore, the K10 chain length composition shows the higher value of binding energy and negative bulk modulus of 6267.16 kJ/mol and −3709.54 GPa, respectively. Likewise, the K10 chain length possesses a higher cohesive energy density and solubility parameter of (6.885 ± 0.00076) × 10 9 J/m 3 and (82.974 ± 0.005) (J/cm 3 ) 0.5 , respectively, which provides the best characteristics for high bond strength and hydrophobic effects and hence improved corrosion and wear-resistance of the surface. To conclude, this work is helpful and can be instructive for material design in water hydraulics to alleviate erosion and corrosion. Further research in the future should investigate the behavior of PTFE coatings of other substrates, such as tungsten carbide (WC), boron carbide (B4C), etc.