In this study, the molecular dynamics (MD) simulation method was adopted to simulate the kinetic characteristics of ice under different temperature conditions, analyze the variation law of the intermolecular structure of ice with temperature, and evaluate the variation law of hardness of ice with temperature from the perspective of chemical hardness. The TIP4P/Ice model was selected for the numerical simulations, which is a typical four-water point water model. The TIP4P model is compared with the conventional three-point water model, which adds a virtual atomic point M near the angle bisector O of the original three-point H-O-H model. Espinosa et al. (2016) used the TIP4P/Ice model to simulate homogeneous water nucleation and found that the model exhibited the best simulated water nucleation rate. The TIP4P/Ice model can emulate the overall phase diagram of ice with different densities, which can reconstruct the solid phase of ice. It is a water model that is specially designed to describe a range of properties of ice (Abascal et al., 2005).

Under different conditions, liquid water and water vapor form ice crystals with various structures. In nature, there are two forms of ice polycrystals: ordinary hexagonal ice (ice Ih) and cubic ice (ice Ic), both of which are molecular crystals with intermolecular forces including hydrogen bonds and van der Waals forces; and the transition from ice Ic to ice Ih occurs at approximately 185–190 K (Celli et al., 2020). As a result, the ice structure selected for the numerical simulation in this study was limited to 188 K. The ice Ih-type was used for the numerical simulation at temperatures greater than 188 K, while the ice Ic-type was used for the simulation at temperatures below 188 K. The ice model used for the simulation is illustrated in Figure 1.

The MD simulations were performed using the GROMACS 2020.4 software package (Hess et al., 2008; Habasaki et al., 2017), and a 4 nm × 4 nm × 4 nm cubic simulation box was first established. The TIP4P-ICE model was selected to fill the entire simulation box, and the total number of water molecules in the simulation process was 2,250. The control equations used in the simulations are as follows: m i r = i F i = ∑ j ≠ i f i j (1) F i = − Δ r i U r 1 , r 2 , k r N (2) U r 1 , r 2 , k r N = ∑ i = 1 N ∑ j = 1 ϕ r i j (3) where N denotes the number of atoms in the system, m i denotes the mass of the ith atom, r i denotes the coordinate mass of the atomic position, F i is the force on the atom, f i j is the force of atom j on atom i, ϕ r i j denotes the potential energy of interaction between atom i and atom j, and U r 1 , r 2 , k r N is the total potential energy of the system.

The OPLS-AA force field designed by Doherty et al. (2017) was used to deal with intermolecular interaction forces, and the system used the NVT system synthesis, i.e., the number of particles N, the volume V and the temperature T of the system are kept constant (Fraternali, 1990). Furthermore, the system temperature is controlled by the nose-hoover method, which is based on statistical mechanics and is able to simulate real physical effects, the system pressure is controlled by the Parrinello-Rahman homogeneous pressure control method, which simulates the variable size and shape of a single cell (Saito et al., 2003; Podio-Guidugli, 2010), PME algorithm for electrostatic interaction force calculation, Coulomb cut-off radius of 1 nm (Kholmurodov et al., 2000), and van der Waals cut-off radius of 1.2 nm. The total calculation time was 20 ns, and the time step of the simulation was 2 fs. The coordinate data were saved every 5 ps, and the data from the last 10 ns were used for the analysis to ensure the equilibrium of the system.

In this study, kinetic simulations of ice were performed at nine temperature conditions (153, 163, 173, 183, 193, 203, 213, 233, and 253 K). The hardness reflects the strength of the interparticle bonds of a crystal and it depends on the internal structure of the substance. The strength of the chemical bonding of a substance affects the hardness of the material with the same structure. Consequently, this study analyzed the law of the temperature dependence of ice hardness by analyzing the law of the temperature dependence of hydrogen bonding and van der Waals interactions between the ice molecules.

MD simulations can provide output parameters such as the hydrogen bond angle, van der Waals force, electrostatic force, and O-H radial distribution. The O-H radial and hydrogen bond angle distributions reflect the strength of hydrogen bonding. The data of the O-H radial distribution function, hydrogen bond angle, van der Waals force, and electrostatic force after 10 ns were used to analyze the variation law of intermolecular hydrogen bonding and van der Waals force with temperature, in turn evaluating the law of ice hardness variation with temperature.

The Mohs and Richter Hardness Scales were used to measure the hardness of ice at different temperatures. The Mohs hardness scale measurement principle uses the hardness plan engraving and scratch hardness test from low to high to observe if the measured sample plane has a scratch. If there is no scratch, the sample hardness is greater than that of the hardness tester, thereafter, the higher level of the hardness tester is tested, until the tester can scratch the plane of the sample and there is no difference between the hardness levels of the sample and the tester. Then the hardness of the sample is comparable to the value of the hardness tester.

The measurement range of the Mohs hardness is wide and the measured value is non-absolute. In this study, the Richter hardness scale was selected for the secondary measurement of hardness of ice based on the Mohs hardness classification. The Richter hardness measurement principle involves the use of an impact body with a certain mass to impact the surface of the ice specimen under the action of fixed impact/test force. The impact and rebound velocities at 1 mm from the impact body to the surface of the specimen were measured. Using electromagnetic principles, induction, and speed proportional to the voltage, the hardness value is calculated according to Eq. 4. H L = 1000 × V B V A (4) where HL is the hardness value on the Richter hardness scale, V _B is the rebounding velocity of the impactor, and V _A is the impact velocity of the impactor body. In this study, the Richter hardness measurement experiment was performed using a D-type universal impact device.

Hardness is one of the critical mechanical properties of solids, and it uses certain substances as standards; specifying their hardness from 1 to 10 constitutes the Mohs hardness. Using these standard substances to scratch the surface of the material, the degree to which the object resists external mechanical action to produce scratching is called scratching hardness. When a solid is scratched, the lattice is broken, and the greater the lattice energy, the greater the resistance required to scratch and the harder the crystal. In this study, a Mohs hardness scale was used to measure the hardness of ice by scratching ice specimens. The other equipment required for the experiment included a Richter hardness scale, a temperature sensor, and a mold for preparing ice cubes (Figure 2). Temperature sensors were used with 1.6 and .254 mm diameter thermocouple beads head respectively, and all thermocouples were sampled at 10 Hz using a data logger, with an error range of ±.5°C, all of which can meet the requirements for water freezing time and maintain the temperature inside the ice-forming cavity surface to determine the surface temperature of the ice. The temperature of the ice particles is monitored by a thermocouple, and the mold temperature is controlled by the liquid nitrogen flow rate. Liquid nitrogen flow rate controlled by solenoid flow meter and temperature sensor linkage, the error range of the solenoid flow meter is ±.5%, measurement range is 1.19–119.4 m³/h.

Based on the experimental requirements, the ice-making mold was placed in a liquid nitrogen environment, the temperature sensor was placed in the mold, and distilled water (prepared in advance and cooled to 20°C) was injected into the mold. The amount of liquid nitrogen was adjusted according to the temperature of the prepared ice. The size of the ice prepared by the mold was 4 cm × 4 cm × 2 cm (L × W × H), and the as-prepared ice cube is shown in Figure 3. The experimental scheme for hardness testing of the as-prepared ice cubes at different temperatures shown in Table 1.

Measuring the hardness of ice at different temperatures requires the preparation of ice at different temperatures. There are four existing methods of ice preparation: 1) the broken ice preparation, 2) the vacuum rapid cooling method, 3) direct contact refrigeration, and 4) subcooling dynamic icing methods. Ice prepared by the ice-breaking method is irregular in shape and the ice temperature is not easy to control, while that prepared by the vacuum rapid cooling method has serious bonding issues and poor practicality; ice prepared by direct contact refrigeration is prone to cracking and affects the hardness test results and that by the subcooling method has a low generation rate and high equipment requirements.

Thus, in this study, based on these methods, liquid nitrogen was used as the refrigerant, the method of indirect contact with the refrigerant was selected, and a square ice making mold with an expansion effect was used to prepare the ice (Figure 4). The liquid nitrogen flow rate controls the temperature of the formed ice, while the temperature sensor displays the temperature of the ice. The size of the square mold was 4 cm × 4 cm and the bottom was made of a steel plate with a size of 7 cm × 7 cm as the ice-forming substrate, having good thermal conductivity and a significant ice-forming effect. The physical parameters of the substrate used in this experiment are listed in Table 2. The ρ, λ and C p in the table represent the density, thermal conductivity and constant pressure specific heat of the material used to make the ice cavity substrate, respectively.

The hydrogen bond angle of a water molecule refers to the angle between the OH and the H linkages in one water molecule and the O linkage of neighboring water molecule. The structural diagram corresponding to the geometric criterion for hydrogen bond determination is shown in Figure 5, where β is the hydrogen bond angle, which can be calculated using Eq. 5. β = arcsin R s i n α R 2 + d 2 − 2 R d c o s α . (5)

The hydrogen bonding angle reflects the strength of the hydrogen bond (Hakala et al., 2006). The smaller the β value, the stronger the hydrogen bond. The probability distribution of the bond angles is approximately a normal distribution, and the sharpness of the distribution peak shape indicates the degree of concentration of the bond angle distribution. The probability distributions of the hydrogen bond angles at different temperatures are shown in Figure 6. It can be observed from the figure that the largest and sharpest peak occurred at a bond angle of 4.5°at T = 153 K with a probability of .128637, indicating that the hydrogen bond angle is mostly concentrated at approximately 4.5°. As the temperature increased, the peaks of the bond angle distribution decreased, peak position shifted to the right, probability distribution curve widened, and peak state tended to flatten, indicating that the angle of the bond angle concentration distribution increases with increasing temperature, and the angle of intermolecular hydrogen bonding becomes greater. The hydrogen bonding was weakened when β increased, which is consistent with previous literature (Yang et al., 2019).

The radial distribution function (RDF) is a physical quantity that reflects the microstructural characteristics of the fluid. Each peak of the RDF represents a coordination shell of the reference particle, and the r value corresponding to the peak indicates the range of the coordination particles at the nearest distance from the target particle. The RDF is a description of the spherically localized distribution of other atoms around an atom and is expressed as shown in Eq. 6: g r = < 1 ρ N ∑ i = 1 N ∑ j = 1 N δ r i j − r > (6) where N is the number of water molecules in the system, ρ N is the average density of water, r i j is the distance between molecules i and j, and the term within the angular brackets <> denotes the average system synthesis.

Figure 7 shows the RDF of O-H at different temperatures. According to the definition of RDF, the r values corresponding to the first, second, and third peaks of g _O-H (r) in Figure 7 characterize the distance of hydrogen and oxygen atoms between water molecules and the water molecules of their first, second, and third coordination circles. It can be observed from Figure 7 that there are differences in the peak and trough values of g _O-H (r) under different temperature conditions. The g _O-H (r) peak size represents the strength of hydrogen bonding; the larger the peak, the stronger the hydrogen bonding (Chialvo and Cummings, 1994).

At 153 K, the OH radial distribution function exhibited the first, second, and third peaks at r values of .308, .363, and .45 nm, indicating that the distances of O-H between the central and the nearest neighboring water molecule, and water molecules of second and third coordination circles at 153 K were .308, .363, and .45 nm, respectively. The first peak of the RDF of O-H appeared at r = .308 nm with peak values of 2.538, 2.465, 2.399, 2.321, 2.263, 2.228 at T = 163, 173, 183, 193, 203, and 213 K, respectively; when the temperature was increased to 233 and 253 K, the first peaks of the RDF of O-H were observed at r = .31 nm with peak values of 2.119 and 2.044 at 233 and 253 K, respectively. Experimental evidence of neutron diffraction showed that as the temperature increases, the first peak of the RDF of O-H decreases, the peak position shifted to the right, and the hydrogen bonding between the water molecules weakens (Bruni et al., 1996). From Figures 8A, B, it can be observed that the first and second peaks of the RDF of O-H gradually decreased and the peak position shifted to the right as the temperature increased, which is consistent with the previously reported neutron diffraction experimental results, indicating that the hydrogen bonding weakens when the temperature increases. Similarly, as shown in Figure 8C, the third peak changes in the same pattern as the first and second peaks.

At 153 K, the RDF of O-H showed trough values at r = .337 and .396 nm, is .505 and .366, where the first trough value corresponds to the r value indicating the maximum distance of intermolecular oxygen-hydrogen formation of hydrogen bonds, and this distance is an important criterion to define and calculate the number of hydrogen bonds by geometric methods (Kalinichev and Bass, 1994). Figure 8D shows that the first valleys of the RDF of O-H at 163,173, 183,193, 203, 213, 233, and 253 K were .539, .574, .607, .65, .677, .684, .736, and .76, respectively. As the temperature increased, the value of O-H radial distribution valley increased and the maximum distance of intermolecular hydrogen-oxygen formation of hydrogen bonds also increased; the increase in O-H distance leads to the weakening of intermolecular hydrogen bonding and the decrease in the number of hydrogen bonds formed. Comprehensive analysis of the O-H radial and bond angle distributions under different temperature conditions showed that as the temperature increased, the movement between the water molecules intensified, the hydrogen-oxygen distance between the water molecules and the water molecules of their coordination circles increased, the network structure of hydrogen bonds was destroyed, and the hydrogen bonding between the water molecules weakened.

The intermolecular interactions of ice include hydrogen bonding and van der Waals forces, where the latter are the effects that exist between water molecules, and are different levels of electron correlation. The existence of van der Waals forces can lead to the coalescence of molecules under certain conditions, which is also the reason for water to exists in a multimeric state. The van der Waals forces arise from electrostatic interactions between molecules or atoms. To clarify the variation of intermolecular van der Waals interactions with temperature, this study analyzed the van der Waals forces and electrostatic forces under different temperature conditions. Figure 9 shows the variation of van der Waals and electrostatic forces at different temperatures.

As shown in Figure 9, the van der Waals and electrostatic forces decreased with increasing temperature; at 153 K, the maximum of van der Waals and electrostatic forces were 56,402.34032 and 244,229.2604 kJ/mol, respectively. The consistent pattern of van der Waals and hydrogen-bonding interactions between ice molecules with temperature variation indicates that the intermolecular forces between ice molecules weaken as the temperature increases. The intermolecular forces affect the lattice energy of the molecules; the smaller the forces, the smaller the lattice energy, lesser the energy required to destroy the crystal, and lower the hardness of the crystal.