Icing is a complex phase change process that is widespread in nature and industry and may have a number of negative effects. During the freezing of water into ice, air bubbles are often trapped in ice and affect the physical properties of the ice. To control the icing process, it is necessary to study these air bubbles in ice. Here, an experimental setup is built to study the growth and distribution characteristics of trapped air bubbles. The results show that the critical freezing rates for the transitions from the egg-shaped bubble region to the egg-/needle-shaped bubble region and from the egg-/needle-shaped bubble region to needle-shaped region are 22.45 ± 3.24 and 12.64 ± 1.65 μm/s, respectively. A mathematical model that can predict bubble growth is obtained by coupling the gas diffusion equation, Henry's law, and the Young–Laplace equation. The model shows that both the maximum width of the bubble and the distance between adjacent bubbles mainly depend on the freezing rate and are proportional to the inverse of the second power of the freezing rate, meaning that the maximum width and the distance gradually increase as the freezing rate decreases. These results contribute to a better understanding of icing mechanisms and inform the optimization of anti-icing and deicing methods.
NOMENCLATURE
Variables
- A
-
Area of water–air boundary (m3)
- AR
-
Aspect ratio (/)
- c
-
Air concentration (mol/m3)
- H
-
Henry law constant (Pa m3/mol)
- k
-
Diffusion constant of air in water (m2/s)
- Li
-
Ice length (m)
- Le
-
Egg-shaped bubble length (m)
- Ln
-
Needle-shaped bubble length (m)
- P
-
Pressure (Pa)
- t
-
Time (s)
- V
-
Volume of bubble (m3)
- r
-
Radius of the part of the bubble in water (m)
- We
-
Egg-shaped bubble width (m)
- Wn
-
Needle-shaped bubble length (m)
Greek symbols
Abbreviations
I. INTRODUCTION
Ice is widely present in everyday life and industry and to some extent influences social development.1,2 Ice can be used as a preservative to extend the lifetime of food.3,4 It is also used as a cold compress to treat sprains.5 While ice brings us convenience, it can also have some negative effects, especially when it is in inappropriate places. For example, ice and frost on the outdoor heat exchanger of an air source heat pump can reduce the operating efficiency of the system.6,7 When it is present on the wind turbine blades, the operational load on the system increases and the efficiency of power generation decreases.8,9 Once the aircraft is covered in ice, the sensors will not work properly and the aircraft will be less maneuverable.10 To eliminate and reduce the negative effects of ice cover on these devices, numerous anti-icing methods are used, including mechanical deicing,11 thermally melt ice,12,13 and superhydrophobic surface.14,15 However, these deicing methods have not been able to achieve perfect results for various reasons. Understanding icing processes and properties is the first prerequisite for icing control. Therefore, to optimize and improve anti-icing methods, it is necessary to study the icing.
The typical solidification process of a water droplet can be divided into four stages, including supercooling, recalescence, freezing, and cooling, based on the changes in temperature and morphology it undergoes during the process.16,17 Nucleation of ice crystals occurs spontaneously in the supercooled droplet, at which point it enters the recalescence stage.18 This stage usually lasts only a few tens of milliseconds before entering the freezing stage. At this point, the droplet transforms into a mixture of ice and water, and the temperature remains at the freezing temperature.19 Compared to the other three stages, the freezing stage is longer in duration and more distinctive in character and has a somewhat greater impact on the ice formed by freezing.20,21 Bubbles trapped in ice are a distinctive feature of the freezing process.
Air is less soluble in ice than it is in water. As water freezes into ice, large amounts of air are squeezed out of the ice and converge at the freezing front.22,23 As the freezing front moves, the concentration of air at the freezing front increases. When the air concentration reaches a critical value, bubble nucleation occurs.24,25 Freezing fronts are jagged structures with high gas concentrations in their depressions, so these depressions tend to be sites of bubble nucleation. Due to the higher gas concentration at the freezing front, the bubbles grow gradually. Although the bubbles receive buoyancy in the water, most of them do not escape from the freezing fronts due to adhesion at the interface. As the bubbles grow, the gas concentration at the freezing front gradually decreases, resulting in slower bubble growth. When the freezing front covers the bubbles, the bubbles no longer grow.26,27
Trapped bubbles not only reduce the clarity of the ice,28 but also affect the thermal properties of the ice.29 When bubble ice is thermally melted, the rate of ice melting over the same period of time increases significantly as the porosity in the ice increases.30 In addition, trapped air bubbles significantly reduce the mechanical properties of ice, with the compressive strength of ice with 10% bubble content being only 25% of the compressive strength of clear ice.31,32 Given the impact of bubbles on the physical properties of ice, it is necessary to conduct research. Most of the studies on the growth and distribution of trapped bubbles are based on experiments, while there are fewer models of bubble growth. Compared to experimental research, mathematical modeling can predict and control bubbles. In terms of mathematical modeling, Chu et al. focused on frozen droplets and built a model that can predict the formation and influencing factors of bubbles.33 A model for predicting the growth and distribution of bubbles in ice has not yet been developed.
To optimize and innovate anti-icing and deicing techniques, it is necessary to study the freezing properties of water. In view of the important influence of trapped bubbles on the physical properties of ice, an experimental and mathematical modeling approach has been used to investigate the growth and distribution characteristics of trapped bubbles. In this study, it is proposed to investigate them using a combination of experiments and mathematical modeling. To better observe these properties, a transparent Hele–Shaw cell is chosen as the research equipment for this experiment. The growth and distribution characteristics of trapped bubbles include bubble shape and aspect ratio. This study can not only scientifically improve the understanding of ice, but also provide a reference for accurate prediction and efficient use of ice.
II. METHODOLOGY
To quantify the process of trapped air bubbles in ice during freezing, an experimental study will be carried out. As shown in Fig. 1, a horizontal Hele–Shaw cell is placed on the stand and a camera for recording the experimental procedure is placed above it. A cold light source is placed on the lower side of the Hele–Shaw cell to increase the ambient brightness. The temperature of the Hele–Shaw cell is regulated by varying the power supply through a digital thermostat. A low-temperature glycol solution flows from the constant temperature water bath into a micro heat exchanger to carry away heat generated during the freezing. Temperature and humidity sensors are connected to the data acquisition instrument to record ambient conditions.
Schematic diagram of the experimental setup. (a) Experimental equipment connection diagram and (b) schematic diagram of the freezing device. The red and blue arrows indicate the directions of camera movement and freezing, respectively. (c) Layout diagram of six thermocouples on a copper sheet.
Schematic diagram of the experimental setup. (a) Experimental equipment connection diagram and (b) schematic diagram of the freezing device. The red and blue arrows indicate the directions of camera movement and freezing, respectively. (c) Layout diagram of six thermocouples on a copper sheet.
To facilitate the observation of the growth and distribution of trapped air bubbles in ice, the presence of a single layer of bubbles within the Hele–Shaw cell is ideal. To determine the thickness of the Hele–Shaw cell, the effect of cell thickness on bubble width is experimentally explored. The results showed that the maximum width of the bubbles is less than 0.25 mm when the cell thickness is varied in the range of 0.1–2 mm. Therefore, the thickness of the Hele–Shaw cell is determined to be 0.5 mm. The Hele–Shaw cell used in this experiment consisted of two transparent acrylic sheets with a length of 50 mm, a width of 50 mm, and a thickness of 5 mm and a 0.5 mm thick acrylic spacer. The acrylic spacer is placed in the middle of the two acrylic plates and coated with glue where it touched. To freeze the water inside the Hele–Shaw cell, a copper sheet with a length of 50 mm, a width of 40 mm, and a thickness of 0.5 mm was inserted into one side of the Hele–Shaw cell.
The top of the copper sheet is covered with insulation cotton to minimize heat exchange with the environment. The thermoelectric cooler is in contact with the lower side of the copper sheets to control its temperature. To better control the freezing temperature of the experiment and to monitor whether the temperature distribution on the copper sheet is uniform or not, six thermocouples were placed on the copper sheet close to the Hele–Shaw cell, and their arrangement is shown in Fig. 1(c). The area where the two are in contact is sealed with thermally conductive silicone to improve the heat transfer coefficient. The lower side of the thermoelectric cooler is in contact with the micro heat exchanger to carry away the heat generated during operation. Also, the area where the two are in contact is sealed with thermally conductive silicone. The low-temperature glycol solution passes through the micro heat exchanger to take away the heat generated during the operation of the thermoelectric cooler. Pictures from the experiment are recorded by a camera carrying a macro lens. The microscope objective will move with the freezing front.
All materials are treated before being assembled into the Hele–Shaw cell. The copper sheet is first polished with 2000 mesh sandpaper, then sequentially ultrasonically cleaned with de-ionized water and ethylene glycol reagent for 15 min, and then placed in a blast oven at 80 °C for 5 min. The acrylic sheets are then cleaned with de-ionized water and ethylene glycol reagent for 15 min. The acrylic sheets are ultrasonically cleaned with de-ionized water and glycol reagent for 15 min and then placed in a blast oven at 80 °C for 5 min. The contact angles of the copper and acrylic sheets are 74 ± 3° and 85 ± 3°, respectively. Table I lists the description of all the equipment used in this experimental study. To ensure accuracy and reproducibility of the experiments, saturated de-ionized water that is left to stand for 24 h at room temperature of 20.2 ± 0.2 °C is employed as the material. To prevent the surface of the cell from fogging and affecting the observation of the experiment, the relative humidity of the environment should not be greater than 30.5 ± 3%. The constant −25 °C freezing temperature in the experiment is controlled by a thermoelectric cooler. During the experiment, the copper sheet was lowered from the room temperature of 20.2 °C to the freezing temperature of −25 °C in 45 s, with an average cooling rate of 1 °C/s. At this time, the height of the ice was about 2 mm, and the bubbles in this region were not selected for analysis in order to ensure the accuracy of the experimental results.
Descriptions of the experimental equipment used in this study.
Item . | Equipment . | Brand . | Unit type . |
---|---|---|---|
1 | Blast drying oven | SUPO | 101-0AB |
2 | Camera | Canon | EOS 90D |
3 | Cold light | BETICAL | ULP-L20-SL |
4 | Contact angle meter | AISRY | ASR-709B |
5 | Constant temperature water bath | HMKXYQ | HMBD-4030 |
6 | Date acquisition instrument | Agilent | 34971A |
7 | Digital thermostat | OMRON | E5CWL-Q1TC |
8 | Power supply | MRINS | MSP6700 |
9 | Temperature and humidity senor | HUAHUANWEI | TH10S-B |
10 | Ultrasonic cleaning machine | Yu clean | AK-009SD |
Item . | Equipment . | Brand . | Unit type . |
---|---|---|---|
1 | Blast drying oven | SUPO | 101-0AB |
2 | Camera | Canon | EOS 90D |
3 | Cold light | BETICAL | ULP-L20-SL |
4 | Contact angle meter | AISRY | ASR-709B |
5 | Constant temperature water bath | HMKXYQ | HMBD-4030 |
6 | Date acquisition instrument | Agilent | 34971A |
7 | Digital thermostat | OMRON | E5CWL-Q1TC |
8 | Power supply | MRINS | MSP6700 |
9 | Temperature and humidity senor | HUAHUANWEI | TH10S-B |
10 | Ultrasonic cleaning machine | Yu clean | AK-009SD |
III. RESULTS AND DISCUSSION
A. Bubble formation and distribution characteristics
Trapped air bubbles in ice can be classified into egg-shaped (ESB) and needle-shaped (NSB) bubbles when the aspect ratio of bubble length to width is greater and smaller than 5.24 The growth process is the same for both types of bubbles. To facilitate the analysis of the formation process of trapped bubbles, the egg-shaped bubble (ESB) is a suitable object of study due to its small size, short growth time, and distinctive change characteristics. Normalized dynamic width and length of egg-shaped bubbles with normalized growth time are analyzed, and the results are displayed in Fig. 2. It is seen that one part of the bubble is in the ice and the other in the water before it is trapped in the ice. The dynamic width of the bubble is the size at the freezing front, while the length is the size trapped in the ice. As the normalized bubble growth time increases, the bubbles width shows a trend of first increasing and then decreasing, while the bubbles length shows a linear increase trend. This may be due to the fact that the growth rate in the length direction remains almost identical to the freezing rate until the bubbles are completely trapped in the ice, thus presenting little apparent change over a short period of time. The portion of the bubble in the water grows rapidly due to the higher surrounding gas concentration. The growth rate of the bubble in the width direction slows down due to a decrease in the gas concentration in front of the freezing front caused by the growth of its own and surrounding bubbles. Therefore, the bubble growth process after nucleation can be divided into an expansion stage and a contraction stage based on the dimension change in bubble width.
Variations of normalized dynamic bubble width and dynamic bubble length with normalized bubble growth time. The width increases significantly and then decreases, while the length almost increases linearly. According to the variation of the width, the growth process can be divided into expansion and contraction stages.
Variations of normalized dynamic bubble width and dynamic bubble length with normalized bubble growth time. The width increases significantly and then decreases, while the length almost increases linearly. According to the variation of the width, the growth process can be divided into expansion and contraction stages.
An interesting phenomenon can be observed when the bubble is completely closed. That is that water seems to be present from the place where the bubble edges the ice. These phenomena may be related to a phenomenon known as interfacial pre-melting.34,35 At constant pressure, the melting point of a substance is generally fixed. Instead, when the boundary of a substance is deformed, the melting point of the surrounding substance decreases and is proportional to the local curvature of the particle and the surface energy of the solidification front. This localized decrease in melting temperature at the expense of interfacial capacity is known as the Gibbs–Thomson effect.36,37 The growth of trapped air bubbles leads to a change in the shape of the boundary of the ice, so that there is a water film on the boundary between the bubble and the ice. However, the occlusion of substances, such as ice crystals, makes it difficult to observe this phenomenon in the experiment. When the bubble is closed, the pre-melting phenomenon is more pronounced due to the large change in the radius of curvature of the bubble boundary (Multimedia view). As demonstrated in Fig. 3, after the bubble is completely captured by the freezing front, the curvature of the in-ice portion of the bubble increases, leading to a decrease in the melting point of the ice in the local region. Subsequently, as the temperature decreases, the film of water that eventually surrounds the top of the bubble also freezes.
(a) An experimental image in which air bubble is captured by freezing front. (b)–(d) Schematic diagram of the process of a bubble being trapped in ice. A water film between the bubble and ice can be observed during the process of the bubble being frozen in ice. Multimedia available online.
(a) An experimental image in which air bubble is captured by freezing front. (b)–(d) Schematic diagram of the process of a bubble being trapped in ice. A water film between the bubble and ice can be observed during the process of the bubble being frozen in ice. Multimedia available online.
Similarly, in the early stage of trapped bubble formation, a film of water exists where the bubble is in contact with the ice due to the greater curvature of the bubble. Although as early as 1967, Maneno demonstrated the presence of water at their boundary when a moving front comes into contact with an air bubble by pre-positioning a bubble in front of a frozen front.38 However, there may still be some problems with this experiment and the results. The temperature inside the bubble is higher than the freezing temperature, causing the ice around the bubble to melt. The second is that the experiment used implanted bubbles rather than bubbles that were subjected to being trapped in ice. In addition, microflow in front of the freezing front may also be responsible for bubble escape.39 Nakaya also demonstrated the existence of a film of water on the surface of a hockey puck through a simple hanging hockey puck experiment.40 The experiment appears to be affected by steam in the air. Therefore, to clarify whether a water film exists between the bubble and the ice, direct experimental observation is a better means of verification.
The distribution of individual bubbles in the ice after formation is an important factor affecting the physical properties of the ice. The study of experimental photographs reveals that there is a certain pattern in the distribution of the two types of bubbles in ice. Egg-shaped bubbles and needle-shaped bubbles will appear in the ice sequentially. As shown in Fig. 4, according to their distribution in the ice, the ice can be divided into three different regions, i.e., an egg-shaped bubble region (ESBR), an egg-/needle-shaped bubble region (E&NSBR), and a needle-shaped bubble region (NSBR). As the ice height increases, the thermal resistance increases gradually. In addition, since the thermal conductivity of air bubbles in ice is relatively small, an increase in air bubbles will also increase the thermal resistance as well. Therefore, the freezing rate decreases as the ice height increases. A further study found that the dimensions of the three bubble regions are 2782.89 ± 301, 4140.2 ± 324, and 3, 792.96 ± 59 μm. Thus, the percentages of the three bubble regions are 25.92 ± 1.16%, 38.61 ± 0.56%, and 35.47 ± 2.18%. It is also found that a correlation between the behavior of the bubbles and the rate of freezing independent of the freezing temperature could be observed. The freezing rates for the transitions from ESBR to E&NSBR and from E&NSBR to NSBR are 22.45 ± 3.24 and 12.64 ± 1.65 μm/s, respectively.
Morphology of trapped bubbles at different locations in ice. Egg-shaped bubbles and needle-shaped bubbles appear sequentially in ice. Based on their distribution in the ice, the ice can be divided into an egg-shaped bubble region (ESBR), an egg-/needle-shaped bubble region (E&NSBR), and a needle-shaped bubble region (NSBR).
Morphology of trapped bubbles at different locations in ice. Egg-shaped bubbles and needle-shaped bubbles appear sequentially in ice. Based on their distribution in the ice, the ice can be divided into an egg-shaped bubble region (ESBR), an egg-/needle-shaped bubble region (E&NSBR), and a needle-shaped bubble region (NSBR).
B. ESB and NSB distribution in ice
To study the distribution characteristics of bubbles in ice, it is necessary to analyze the size of two types of bubbles in each of the three bubble regions. The width, length, and aspect ratio of two types of bubbles are selected, and their changes along the ice length direction are shown in Figs. 5 and 6. As demonstrated in Fig. 5(a), the ESBs widths have the same trend in both ESBR to E&NSBR and gradually increases along the ice length direction. The slope of the fitted straight line for ESBs from all experimental points in the ESBR is 21.39, while the slope of the fitted straight line for ESBs from all experimental points in the E&NSBR is 17.54. Due to the high number of egg-shaped bubbles in the ESBR, the growth process of adjacent bubbles is similar, resulting in a smaller difference in the size. In E&NSBR, the growth of egg-shaped bubbles is influenced by needle-shaped bubbles, which are smaller in size. In addition, due to the slowing down of the freezing process, the growth of the bubbles takes a longer time, and thus, its size is larger. This results in a more dispersed size distribution of egg-shaped bubbles in the E&NSBR. The freezing rate decreases as the ice length increases, which leads to a prolongation of the bubble growth time, while the air concentration in front of the freezing front decreases. Especially in E&NSBR, the appearance of NSBs further reduces the air concentration in front of the freezing front. This leads to a gradual increase in the rate of the bubble width.
Distribution of egg-shaped trapped air bubbles trapped in ice. Bubble (a) width, (b) length, and (c) aspect ratio (AR) of bubbles in different ice regions. The blue and brown lines are fitted from all experimental points.
Distribution of egg-shaped trapped air bubbles trapped in ice. Bubble (a) width, (b) length, and (c) aspect ratio (AR) of bubbles in different ice regions. The blue and brown lines are fitted from all experimental points.
Distribution of needle-shaped air bubbles trapped in ice. (a) Bubble width, (b) length, and (c) aspect ratio (AR) of bubbles in ice. The blue and brown lines are fitted from all experimental points.
Distribution of needle-shaped air bubbles trapped in ice. (a) Bubble width, (b) length, and (c) aspect ratio (AR) of bubbles in ice. The blue and brown lines are fitted from all experimental points.
Similar to the variation trend of bubble width, the bubble length gradually increases along the ice length direction. As seen in Fig. 5(b), the slopes of the fitted lines in the two bubble regions are 58.21 and 67.52, respectively. On the one hand, the gradual decrease in the freezing rate results in a longer growth time for bubbles. On the other hand, it may also be due to the decrease in the number of bubbles in E&NSBR, and so more air is absorbed and the bubbles continuously grow along the length direction. In addition, the aspect ratio of bubbles length to width also shows a gradual increase along the ice length direction, although there is still a significant difference in the aspect ratio (AR) between the two bubble regions. The slopes of the fitted lines in the two bubble regions are 0.44 and 0.24, respectively. The AR value can comprehensively reflect the difference between bubbles in two bubble regions. The experimental results indicate that the expansion of bubble width is greater than the elongation of bubble length in E&NSBR. This indicates that the decrease in the freezing rate leads to an extension of bubble growth time, which has a greater impact on width than on length.
The variations of the NSB sizes along the ice length direction are also analyzed and shown in Fig. 6. The NSB sizes also gradually increase along the ice length direction. The slopes of the fitted lines for bubble width in the two bubble regions are 9.07 and 16.30, respectively. The slopes of fitted lines within the two bubble regions for bubble length are 0.25 and 0.16, respectively. The variation trend of these two sizes along the ice length direction is the same as that of ESBs. In E&NSBR, the reduction in the freezing rate allows more time for the bubbles to grow in the width direction. The decrease in the freezing rate leads to a decrease in the air concentration in front of the freezing front; thus, the growth of the bubble length gradually slows down. The AR values of NSBs can verify this very well. As shown in Fig. 6(c), the slopes of the fitted lines for the AR values are 1.82 and −1.47 in the two bubble regions, respectively. Although the freezing rate decreases significantly with an increasing ice length. The reduction in the freezing rate results in lower gas concentrations at the freezing front. As a result, the growth of gas bubbles slows down or even stops. At the same time, bubble nucleation may not occur because the gas concentration at the freezing front does not reach the bubble nucleation concentration.
C. Effect of freezing rate on bubble formation and distribution
The results of the experimental studies described above show that the growth and distribution of trapped bubbles in ice are related to the freezing rate. To further investigate the effect of freezing rate on bubble shape, the egg-shaped bubbles are chosen for the study due to their short growth time and significant variation in dimension. The aspect ratios of egg-shaped bubbles in the range of freezing rate from 15 to 75 μm/s are counted. As demonstrated in Fig. 7, the aspect ratio shows a gradual increase as the freezing rate decreases. After nucleation, the tiny bubbles initially assume a perfect spherical shape, and during growth, the bubbles are deformed by the compression of the ice. The larger the freezing rate, the shorter the growth time of the bubbles, and thus, their deformation is smaller and the aspect ratio is smaller. This indicates that the freezing rate has an important effect on both the size and shape of the bubbles.
Aspect ratio of the bubble as a function of the freezing rate. As the freezing rate decreases, the size and aspect ratio of the trapped air bubbles show a gradual increase.
Aspect ratio of the bubble as a function of the freezing rate. As the freezing rate decreases, the size and aspect ratio of the trapped air bubbles show a gradual increase.
Comparison of experimental and modeling values for the relationship between bubble width and freezing rate.
Comparison of experimental and modeling values for the relationship between bubble width and freezing rate.
Variation of normalized bubble width with normalized growth time. The third power of the normalized bubble width is proportional to the normalized growth time. The insets show in detail the size change of a single bubble during the expansion stage. Multimedia available online.
Variation of normalized bubble width with normalized growth time. The third power of the normalized bubble width is proportional to the normalized growth time. The insets show in detail the size change of a single bubble during the expansion stage. Multimedia available online.
Comparison of experimental and modeling values for the relationship between bubble distance and freezing rate. The experimental values suggest that the bubble width is a function of 1/v0.5.
Comparison of experimental and modeling values for the relationship between bubble distance and freezing rate. The experimental values suggest that the bubble width is a function of 1/v0.5.
IV. CONCLUSIONS
To quantify and visualize the bubbles trapped in the ice during the freezing process, an experimental setup is constructed, and a series of experiments are carried out. (1) According to the change in the bubble width, bubble growth process can be divided into two stages of expansion and contraction. (2) The critical freezing rates for the transitions from the egg-shaped bubble region to the egg-/needle-shaped bubble region and from the egg-/needle-shaped bubble region to needle-shaped region are 22.45 ± 3.24 and 12.64 ± 1.65 μm/s, respectively. (3) The characteristic dimensions of egg-shaped and needle-shaped bubbles follow the same trend within the egg-shaped bubble region and the egg-/needle-shaped bubble region. The bubble width and bubble length gradually increase along the ice length direction, while the aspect ratio gradually decreases. (4) Both the bubble width of bubbles and the distance between two bubbles increase with a decreasing freezing rate. The mathematical modeling results show that they are both inversely proportional to the second power of the freezing rate.
ACKNOWLEDGMENTS
The corresponding author acknowledges the financial supports from the National Natural Science Foundation of China (No. 52076013), Beijing Municipal Science & Technology Commission (Nos. 3212024 and No. 3232032), China Postdoctoral Science Foundation (No. 2022M720423), Young Elite Scientist Sponsorship Program by BAST (No. BYESS2023352), and Key Laboratory of Icing and Anti/De-icing of CARDC (No. IADL20220304).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Keke Shao: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Mengjie Song: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Xuan Zhang: Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Long Zhang: Funding acquisition (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.