Numerical simulation analysis method for road snow melting and ice removal

Abstract

This invention belongs to the field of computer numerical simulation technology and relates to a numerical simulation analysis method for road snow melting and icing based on a CFD-DEM coupled model. The method includes: 1) generating a corresponding discrete element road model with corresponding particle composition, and assigning different micromechanical and thermodynamic parameters to the model for concrete and asphalt; 2) establishing an ice model and setting ice parameters, establishing a non-deterministic heat transfer equation for the control volume in the ice model, and using the liquid fraction β to represent the proportion of liquid in the control volume during heat conduction; 3) updating the liquid fraction β, and using the updated value of the liquid fraction β to determine the change in the ice model volume, providing a basis for road snow melting and icing analysis based on this change. This invention provides a numerical simulation analysis method for road snow melting and icing that can formulate optimal control and solution strategies for road snow melting and icing under different conditions.

Description

Technical Field

[0001] This invention belongs to the field of computer numerical simulation technology, and relates to a numerical simulation analysis method for road snow melting and ice removal, and more particularly to a numerical simulation analysis method for road snow melting and ice removal based on a CFD-DEM coupled model. Background Technology

[0002] In winter, especially in frigid regions, icy and snow-covered roads severely impact a country's transportation, economy, and normal outdoor activities and work. The severe snow and ice accumulation often causes traffic disruptions or accidents. According to incomplete statistics, approximately 30% of traffic accidents are related to icy and snowy roads, seriously affecting transportation and economic development. Road icing and snow problems have long plagued transportation departments worldwide. To address this challenge, countries globally primarily employ manual or mechanical methods such as spreading de-icing agents and gravel, and utilize external heating technologies, energy conversion de-icing technologies, and conductive pavement, among other de-icing and snow-melting techniques.

[0003] Carbon fiber electrothermal de-icing is a novel active snow and ice melting technology. Its working principle involves laying carbon fiber heating cables within the road concrete, with two different types of asphalt layers on top. The heating cables convert electrical energy into heat energy, which is then transferred to the ice layer on the road surface through thermal conduction, raising the road surface temperature below 0°C and melting the ice and snow. As a type of de-icing system, the electrothermal de-icing system uses an electric heat source to melt the ice accumulated on the outer surface of icy components. After the electrothermal de-icing system operates for a period of time, the ice will melt. The melting efficiency depends on two factors: the electrothermal power and the thermophysical properties of the ice. Studying the influence of thermophysical properties on the de-icing process helps in designing more suitable de-icing systems. The de-icing process involves many issues, such as temperature distribution in the melting zone, melting rate, melting area, and melting volume. Currently, research on the heat transfer and melting processes of ice has become a hot topic.

[0004] The core of road snow and ice melting is heat exchange, primarily the transfer of heat between layers, involving a three-phase system and anisotropic flow, accompanied by complex heat and mass transfer chemical reactions and strong interphase interactions. Based on their location and function, road surface layers mainly consist of the surface layer, base layer, and subbase layer. The surface layer area is an indispensable and crucial energy reaction zone for simulating the road snow and ice melting process. However, the complexity of the problem arises from the heat transfer within the road and the complex chemical reactions involving the three phases, making it difficult to accurately understand the heat transfer at each layer that leads to the melting of the ice layer.

[0005] To accurately and reliably describe the process of ice melting caused by heat transfer within roads, and to formulate optimal control and solution strategies for road snow and ice melting under different conditions, it is necessary to establish a numerical simulation analysis method for road snow and ice melting based on a CFD-DEM coupled model. Summary of the Invention

[0006] To address the aforementioned technical problems in the background art, this invention provides a numerical simulation analysis method for road snow melting and ice removal that can formulate optimal control and solution strategies for road snow melting and ice removal under different conditions.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] A numerical simulation analysis method for road snow melting and ice removal, characterized in that: the numerical simulation analysis method for road snow melting and ice removal includes the following steps:

[0009] 1) The dimensional parameters of the road model are set using discrete element method (DEM) software. Based on the dimensional parameters and the particle size distribution of the sample, a corresponding discrete element road model with particle composition is generated. Different micromechanical and thermodynamic parameters of concrete and asphalt are assigned to the discrete element road model. Using the thermal module program embedded in the discrete element software, based on the set micromechanical and thermodynamic parameters, each particle in the discrete element road model is regarded as a heat storage device in the heat conduction process, and the temperature change and heat transfer of the particle system in the discrete element road model are calculated.

[0010] 2) Using the OpenFoam software for computational fluid dynamics, an ice model was built on the upper surface of the discrete element road model and ice parameters were set. At the same time, the ice model was meshed. The moving boundary that satisfies the heat and mass balance conditions was solved according to the enthalpy method. The non-deterministic heat transfer equation of the control volume in the ice model was established. The liquid fraction β was used to represent the proportion of liquid in the control volume during the heat conduction process.

[0011] 3) Based on the temperature changes and heat transfer of the particle system within the discrete element road model obtained in step 1), record the temperature changes and accumulated heat of the particles at the upper asphalt boundary of the discrete element road model. Through unidirectional coupling calculation, input the temperature changes and accumulated heat of the particles at the upper asphalt boundary of the discrete element road model as boundary input conditions into the OpenFoam software. Using the non-constant heat transfer equation of the control volume obtained in step 2), calculate the heat transfer within a given load time, update the liquid fraction β, and use the updated liquid fraction β value to determine the change in the ice model volume. Based on this change, provide a basis for the analysis of road snow melting and ice removal.

[0012] Preferably, the discrete element method (DEM) software used in this invention is Yade DEM software. The DEM road model includes an upper asphalt layer, a lower asphalt layer, and a concrete layer arranged sequentially from top to bottom. The thickness of the upper asphalt layer is h1, the thickness of the lower asphalt layer is h2, and the thickness of the concrete layer is h3, where h1 < h2 < h3. In the DEM road model, based on the contact characteristics between particles in the concrete and particles in the asphalt, the contact between particles is determined by the parallel bonding model connectivity.

[0013] Preferably, in step 2) of the present invention, the mesh division is to divide the ice model into computational units of uniform size, calculate the temperature change and heat accumulation inside each computational unit, and divide the mesh size into the size of the smallest diameter particle of the upper asphalt interface of the discrete element road model.

[0014] Preferably, in step 2) of this invention, the specific implementation method for solving the moving boundary that satisfies the heat and mass balance conditions using the enthalpy method is as follows:

[0015] Using enthalpy as a dependent variable, and assuming that enthalpy is a function of temperature, the relationship between enthalpy and temperature is used to determine the temperature of ice during the melting process;

[0016] Temperature ;

[0017] in:

[0018] The T melt It is the temperature at which ice begins to reach its melting point during the melting process;

[0019] The H sm The enthalpy of ice at its melting point;

[0020] The H sm =C ps T melt The H lm The enthalpy of water at its melting point;

[0021] The H lm =C ps T melt +L;

[0022] H represents the total enthalpy of ice;

[0023] The C ps This represents the specific heat capacity of ice under constant pressure.

[0024] The C pl This refers to the specific heat capacity of water under constant pressure.

[0025] When H <H smAt that time, the ice had not yet reached its melting point and had not yet begun to melt;

[0026] When H sm <H<H lm At that time, the ice reaches its melting point but has not yet melted completely. When H = H lm When it reaches the melting critical point;

[0027] When H>H lm At that moment, the ice began to melt.

[0028] Preferably, the expression for the non-constant heat transfer equation of the control volume in the ice model used in step 2) of this invention is:

[0029]

[0030] in:

[0031] The j represents the layer of the ice model built in OpenFOAM;

[0032] The x, y, and z directions are the X, Y, and Z directions in spatial coordinates, respectively.

[0033] The The density of ice;

[0034] The C p This represents the specific heat capacity of ice under constant pressure.

[0035] k is the thermal conductivity of ice;

[0036] T represents the temperature of the ice model during the melting process;

[0037] Q represents the amount of heat source;

[0038] The non-deterministic heat transfer equation for the control volume in the ice model represents the energy accumulation rate of the sum of heat conduction and local heat generation in the X, Y, and Z directions.

[0039] Preferably, the functional expression of the non-constant heat transfer equation of the control volume in the ice model used in this invention within the ice layer is:

[0040]

[0041] in:

[0042] The H ice It is the enthalpy of ice;

[0043] T represents the temperature of the ice model during the melting process;

[0044] The K xice K yice K ziceThese are the thermal conductivity coefficients in the x, y, and z directions of the ice model's spatial coordinates, respectively.

[0045] t is the calculation time.

[0046] Preferably, the liquid fraction β used in this invention is defined as follows:

[0047]

[0048] in:

[0049] The T < T melt At that time, the ice was in an unmelted state;

[0050] The T=T melt At that time, the ice was in an ice-water paste state;

[0051] The T>T melt At that time, the ice was in a completely melted state;

[0052] T represents the temperature of the ice model during the melting process;

[0053] The T melt It is the melting point of ice.

[0054] As a preferred embodiment, step 3) of the present invention is specifically implemented as follows:

[0055] 3.1) Based on the temperature change and accumulated heat of the upper asphalt boundary particles in the discrete element road model obtained in step 1), these are used as boundary input conditions and substituted into the OpenFoam software as initial conditions to determine the temperature field.

[0056] 3.2) Based on the temperature field determined in step 3.1), using the enthalpy method and the non-indeterminate heat transfer equation of the control volume established in the ice model, the temperature change and heat accumulation within the computational cells formed after meshing in the ice model are calculated within a given load time to simulate ice melting; wherein the load step in CFD is 10. -6 s;

[0057] 3.3) In the calculation of ice melting process using the enthalpy method, under the iterative load step, the computational cells formed after the meshing in the ice model continuously absorb heat, update the temperature field and material field, update the liquid fraction β, and finally determine the change in ice melting volume by the updated value of the liquid fraction β.

[0058] The numerical simulation analysis method for road snow melting and ice removal proposed in this invention has the following beneficial effects:

[0059] This invention employs a simplified model for simulation, which shortens computation time, reduces computational costs, achieves high computational accuracy, and has strong universal applicability. It can more easily obtain the distribution patterns of certain microscopic mechanical properties that are difficult to obtain through experimental methods. By recording temperature changes and heat accumulation between particles in the DEM during the simulation process, it can reveal the laws governing temperature and heat transfer between particles. By reading the boundary condition information input from the DEM using OpenFoam, it simulates the heat transfer process in the ice model, further revealing the laws governing ice melting. This numerical simulation method can provide a basis for analyzing road snow melting and ice removal under real-world conditions. Attached Figure Description

[0060] Figure 1 This is a schematic diagram of the road's geometric model;

[0061] Figure 2 This is a diagram showing the distribution of carbon fiber heating cables in the road geometry model.

[0062] Figure 3 This serves as a numerical sample for a road snow melting and ice removal simulation experiment. Detailed Implementation

[0063] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Preferred embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete.

[0064] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0065] To accurately and reliably describe the process of ice melting caused by heat transfer within roads, and to formulate optimal control and solution strategies for road snow and ice melting under different conditions, it is necessary to establish a numerical simulation analysis method for road snow and ice melting.

[0066] Please see Figures 1 to 3 This invention proposes a numerical simulation analysis method for road snow melting and ice removal, which includes the following steps:

[0067] 1. In the DEM, according to the set size parameters and particle size distribution of the road model, generate the corresponding particle composition road model. The size parameters are a length of 1m and a thickness of 0.2m, wherein the thickness from top to bottom is a 0.04m asphalt layer, a 0.06m asphalt layer and a 0.1m concrete layer.

[0068] 2. Based on the road geometry model, different micromechanical and thermodynamic parameters are assigned to the asphalt concrete particles to facilitate the calculation of particle temperature changes and heat transfer in the subsequent thermal module. The density among these micromechanical parameters is 2463 kg / m³. 3 The porosity is 0.11, the damping is 0.7, the effective deformation modulus is 1.0e9Pa, and the time step is 10. -8 The specific heat capacities of concrete, lower asphalt, and upper asphalt in the thermodynamic parameters are 970 J / kg·℃, 921.1 J / kg·℃, and 1168.0 J / kg·℃, respectively, and their thermal conductivityes are 1.28 W / m·℃, 1.163 W / m·℃, and 2.057 W / m·℃, respectively.

[0069] 3. Considering the contact characteristics between particles, which transfer temperature and heat, and taking into account the contact characteristics between asphalt concrete particles, the contact type between particles should be selected from the parallel bonding model embedded in the discrete element software to achieve bonding between particles and provide shear and tensile strength between particles.

[0070] 4. In the embedded thermal module program of the Discrete Element Method (DEM) software, based on the generated particles and the assigned micromechanical and thermodynamic parameters, the temperature and heat of the particles are calculated within each set time step. Each particle in the road model is considered as a heat storage device in the heat conduction process, and the temperature change and heat transfer of the particles are iterated at each time step. The time step in the DEM is 10. -8 s.

[0071] 5. With the iterative accumulation of time steps, the temperature and heat of the particles continuously rise. According to the contact characteristics between particles, the temperature and heat are continuously transferred upwards until they reach the particles on the surface of the upper asphalt layer. The temperature and heat of the particles at this interface are recorded and used as a boundary input condition in the CFD as an initial condition to determine the temperature field and calculate the melting process of ice.

[0072] 6. Create an ice model in OpenFoam, set the appropriate ice parameters, and mesh the ice model, dividing the mesh size to the size of the smallest diameter particles at the upper asphalt interface. The ice parameters are: thickness 0.02 m, density 600 kg / m³, porosity 35%, thermal conductivity 1.45 W / mK, specific heat capacity 1725.65 J / kg·K, thermal diffusivity 1.41 m² / s, and initial temperature -6.5℃.

[0073] 7. In OpenFoam, initial conditions are set based on the temperature changes and heat accumulation of particles on the surface of the upper asphalt layer recorded in the DEM. These are used as boundary input conditions and input into OpenFoam as initial conditions to determine the temperature field and calculate the ice melting process. The time step for calculating ice melting in OpenFoam is 10. -6 s.

[0074] 8. When ice melts into water, two states exist, with ice and water separated by a moving interface. The difficulty in simulating phase transition processes lies in the existence of a moving boundary or region that must satisfy heat and mass equilibrium conditions. The enthalpy method is generally used to solve moving boundary problems. In the solution process, the enthalpy method treats enthalpy as a dependent variable, eliminating the need to determine the interface location. It assumes enthalpy is a function of temperature, using the relationship between enthalpy and temperature to determine the temperature of the ice. The temperature of the ice is:

[0075]

[0076] H sm H is the enthalpy of ice at its melting point. lm Let C be the enthalpy of water at its melting point, H be the total enthalpy of ice, and C be the total enthalpy of ice. ps C is the specific heat capacity of ice under constant pressure. pl T is the specific heat capacity of water under constant pressure. melt It is the temperature at which ice begins to reach its melting point during the melting process.

[0077] 9. Solve the non-steady-state heat transfer equations for the control volume in the ice model using OpenFoam. The non-steady-state heat transfer equations for the control volume in the ice model are:

[0078]

[0079] The subscript j indicates the determined level. Let Cp be the density of ice, Cp be the isobaric specific heat capacity of ice, K be the thermal conductivity of ice, T be the temperature of the ice model during melting, Q be the amount of heat source, and t be the calculation time. This equation represents the energy accumulation rate of the sum of heat conduction in three directions and local heat generation.

[0080] Within the ice layer, the governing equation can be expressed as:

[0081]

[0082] Where Hice is the total enthalpy of ice, Kxice, Kyice, and Kzice are the thermal conductivity in the x, y, and z directions of ice, T is the temperature of the ice model during melting, and t is the calculation time.

[0083] 10. In OpenFoam, the enthalpy method is used to calculate the ice melting process, with the liquid fraction (β) representing the proportion of liquid in the control volume. During iterative load steps, the elements continuously absorb heat and update the temperature and material fields (phase transitions), ultimately determining the change in ice-melting volume using the value of the liquid fraction (β). The ice-water paste region is defined by a liquid fraction between 0 and 1; when the ice is completely melted, the liquid fraction is 1. The liquid fraction β can be defined as:

[0084]

[0085] Where T is the temperature of the ice model during the melting process, T melt It is the melting point of ice.

[0086] Steps 1) to 5) above are calculated in Discrete Element Method (DEM) software; steps 6) to 10) are calculated in OpenFoam.

[0087] It should be noted that in the moving boundary method based on the enthalpy method to solve for the heat and mass balance conditions, a real ice-melting physical model is usually very complex, containing many uncertainties such as dust, density inhomogeneity, and irregular shape. Therefore, for numerical simulation studies, the following assumptions need to be made: The thermal properties of the internal structural materials are assumed to be independent of temperature; the ambient temperature and convective heat transfer coefficient are assumed to be constant and independent of time; perfect thermal contact exists between layers; ice and water are assumed to be homogeneous and isotropic; the effect of surface curvature is assumed to be neglected; and a phase transition occurs at the melting point temperature.

Claims

1. A numerical simulation analysis method for road snow melting and ice removal, characterized in that: The numerical simulation analysis method for road snow melting and ice removal includes the following steps: 1) The dimensional parameters of the road model are set using discrete element method (DEM) software. Based on the dimensional parameters and the particle size distribution of the sample, a corresponding discrete element road model with particle composition is generated. Different micromechanical and thermodynamic parameters of concrete and asphalt are assigned to the discrete element road model. Using the thermal module program embedded in the discrete element software, based on the set micromechanical and thermodynamic parameters, each particle in the discrete element road model is regarded as a heat storage device in the heat conduction process, and the temperature change and heat transfer of the particle system in the discrete element road model are calculated. 2) Using the OpenFoam software for computational fluid dynamics, an ice model was built on the upper surface of the discrete element road model and ice parameters were set. At the same time, the ice model was meshed. The moving boundary that satisfies the heat and mass balance conditions was solved according to the enthalpy method. The non-deterministic heat transfer equation of the control volume in the ice model was established. The liquid fraction β was used to represent the proportion of liquid in the control volume during the heat conduction process. 3) Based on the temperature changes and heat transfer of the particle system within the discrete element road model obtained in step 1), record the temperature changes and accumulated heat of the particles at the upper asphalt boundary of the discrete element road model. Through unidirectional coupling calculation, input the temperature changes and accumulated heat of the particles at the upper asphalt boundary of the discrete element road model as boundary input conditions into the OpenFoam software. Using the non-constant heat transfer equation of the control volume obtained in step 2), calculate the heat transfer within a given load time, update the liquid fraction β, and use the updated liquid fraction β value to determine the change in the ice model volume. Based on this change, provide a basis for the analysis of road snow melting and ice removal.

2. The numerical simulation analysis method for road snow melting and ice removal according to claim 1, characterized in that: In step 1): the discrete element software is Yade discrete element software, and the discrete element road model includes an upper asphalt layer, a lower asphalt layer, and a concrete layer arranged sequentially from top to bottom; the thickness of the upper asphalt layer is h1, the thickness of the lower asphalt layer is h2, and the thickness of the concrete layer is h3, wherein h1 < h2 < h3; in the discrete element road model, based on the contact characteristics between particles in concrete and particles in asphalt, the contact between particles is determined by the parallel bonding model connectivity.

3. The numerical simulation analysis method for road snow melting and ice removal according to claim 2, characterized in that: In step 2), the meshing is to divide the ice model into computational units of uniform size, calculate the temperature change and heat accumulation inside each computational unit, and divide the mesh size into the size of the smallest diameter particle at the upper asphalt interface of the discrete element road model.

4. The numerical simulation analysis method for road snow melting and ice removal according to claim 3, characterized in that: In step 2), the specific implementation method for solving the moving boundary that satisfies the heat and mass balance conditions using the enthalpy method is as follows: Using enthalpy as a dependent variable, and assuming that enthalpy is a function of temperature, the relationship between enthalpy and temperature is used to determine the temperature of ice during the melting process; Temperature ; in: The T melt It is the temperature at which ice begins to reach its melting point during the melting process; The H sm The enthalpy of ice at its melting point; The H sm =C ps T melt The H lm The enthalpy of water at its melting point; The H lm =C ps T melt +L; H represents the total enthalpy of ice; The C ps This represents the specific heat capacity of ice under constant pressure. The C pl This refers to the specific heat capacity of water under constant pressure. When H <H sm At that time, the ice had not yet reached its melting point and had not yet begun to melt; When H sm <H<H lm At that time, the ice reaches its melting point but has not yet melted completely. When H = H lm When it reaches the melting critical point; When H>H lm At that moment, the ice began to melt.

5. The numerical simulation analysis method for road snow melting and ice removal according to claim 4, characterized in that: The expression for the non-steady-state heat transfer equation of the control volume in the ice model described in step 2) is: ， in: The j represents the layer of the ice model built in OpenFOAM; The x, y, and z directions are the X, Y, and Z directions in spatial coordinates, respectively. The The density of ice; The C p This represents the specific heat capacity of ice under constant pressure. k is the thermal conductivity of ice; T represents the temperature of the ice model during the melting process; Q represents the amount of heat source; The non-deterministic heat transfer equation for the control volume in the ice model represents the energy accumulation rate of the sum of heat conduction and local heat generation in the X, Y, and Z directions.

6. The numerical simulation analysis method for road snow melting and ice removal according to claim 5, characterized in that: The functional expression of the non-indeterminate heat transfer equation of the control volume in the ice model within the ice layer is: ， in: The H ice It is the enthalpy of ice; T represents the temperature of the ice model during the melting process; The K xice K yice K zice These are the thermal conductivity coefficients in the x, y, and z directions of the ice model's spatial coordinates, respectively. t is the calculation time.

7. The numerical simulation analysis method for road snow melting and ice removal according to claim 6, characterized in that: The liquid fraction β is defined as follows: ,in: The T < T melt At that time, the ice was in an unmelted state; The T=T melt At that time, the ice was in an ice-water paste state; The T>T melt At that time, the ice was in a completely melted state; T represents the temperature of the ice model during the melting process; The T melt It is the melting point of ice.

8. The numerical simulation analysis method for road snow melting and ice removal according to claim 7, characterized in that: The specific implementation method of step 3) is as follows: 3.1) Based on the temperature change and accumulated heat of the upper asphalt boundary particles in the discrete element road model obtained in step 1), these are used as boundary input conditions and substituted into the OpenFoam software as initial conditions to determine the temperature field. 3.2) Based on the temperature field determined in step 3.1), using the enthalpy method and the non-indeterminate heat transfer equation of the control volume established in the ice model, the temperature change and heat accumulation within the computational cells formed after meshing in the ice model are calculated within a given load time to simulate ice melting; wherein the load step in CFD is 10. -6 s; 3.3) In the calculation of ice melting process using the enthalpy method, under the iterative load step, the computational cells formed after the meshing in the ice model continuously absorb heat, update the temperature field and material field, update the liquid fraction β, and finally determine the change in ice melting volume by the updated value of the liquid fraction β.

Images