Abstract
The long-span ice composite shell structure is a new type of ice and snow structure developed in recent years. The engineering practice of ice composite shell shows that sublimation is one of the important reasons for its damage and even collapse. In this paper, we firstly supplemented the existing H-K equation and obtained the revised ice sublimation equation through indoor evaporative plate experiment considering the influence of admixtures and wind speed. Afterwards, combining the simulations of solar radiation and CFD, the numerical simulation of sublimation distribution on the surface of were realized by programming in Grasshopper platform. During sublimation, the thickness of the ice composite shell decreases by 0.38 mm every 10 days and the sublimation rate on the sunny side was 1.7 times that on the shady side. Finally, the static performance and stability of the sublimated ice composite spherical shell were analyzed. After 70 days of sublimation, the thickness of the ice composite shell structure becomes thinner and uneven, which leads its sensitivity to external load increases.
Keywords
Introduction
As a new type of ice and snow structure, long-span ice composite shell structure has attracted the attention of researchers in recent years.1,2 The fiber reinforced ice (FRI) 3 is widely used in ice composite shell in cold regions. The ice composite structure is formed by mixing a certain proportion of pulp fibers in water to form a homogeneous composite aqueous solution, and then spraying and freezing in low temperature environment through aerated membrane template 4 and spray construction technology. As a kind of thin shell structure, the long-span composite ice shell is highly sensitive to the change of its section thickness, and the small change of thickness may have great influence on its mechanical performance. As shown in Figures 1 and 2,4 with the continuous sublimation, the ice loss on the surface leads to the gradual thinning of the thickness of the structure and the change of its shape. In severe cases, it may cause the perforation damage of the structure during the service period.

Ice shell thickness curve.

Ice shell sublimation perforation failure.
The uneven thickness of thin shell structure caused by sublimation, can seriously affect the structure safety and long-term mechanical properties. It not only brings new problems to the design and analysis of the long-span ice shell structure, but also increases the difficulty in the maintenance of the ice shell structure during its service. Therefore, ice sublimation should not be ignored in the research of large span ice shell structure.
The essence of sublimation is the diffusion of molecules at the microscopic level from a high concentration of solid to a low concentration of air. Scholars5–9 have investigated the effects of vapor pressure, temperature, light source, and impurities on the sublimation rate of pure ice. At the beginning of the 20th century, Miyamoto 10 proposed a theoretical calculation equation for the sublimation rate of pure ice plane in a vacuum environment based on the relevant theories of thermodynamics and gas dynamics. In this equation, the sublimation rate was expressed as the difference between the sublimation process and the condensation process, namely the net sublimation rate. This equation was also known as Hertz-Knudsen-Langmuir equation (HK equation for short). Since then, scholars11–15 have carried out a series of related researches on ice sublimation based on this equation.
The HK equation is applicable to calculate the sublimation rate of pure ice materials in the vacuum environment, and only considers the influence of ice surface temperature and water vapor pressure difference. But it is not applicable to calculate the sublimation rate of pure ice and composite ice materials mixed with other substances in the actual environment.
In this paper, the sublimation calculation theory of pure ice and FRI in the actual environment will be studied. Various factors affecting the ice sublimation rate in the real environment were revealed through different experiments, and the ice sublimation model considering different influencing factors was modified. We also established a numerical simulation method for the cumulative sublimation distribution on the surface of long-span composite ice shell structure.
The paper is organized as follows: Firstly, the sublimation properties of composite ice materials were tested and the revised H-K equation is proposed in Chapter 2. Secondly, the numerical simulation method for surface sublimation of ice shell structure is presented in Chapter 3. Afterwards, the static performance and stability of an ice composite spherical shell with the consideration of sublimation are analyzed in Chapter 4. Finally, the main conclusions are presented in Chapter 5.
Sublimation Experiment of ice materials
Theory and experimental methods of ice sublimation
Suggesting that the molecular escape rate distribution on the ice surface satisfies the Maxwell Boltzmann rate distribution function 16 and ideal gas, the H-K equation were used to calculate the sublimation rate on the pure ice surface under the vacuum environment can be expressed as shown in equation (1) 17 :
For ease of use in practical engineering, Relative humidity of air, which is closely related to saturated vapor pressure and vapor pressure, is introduced in this paper, as shown in equation (2):
The saturated vapor pressure is dependent only on the surface temperature of pure ice, and we choose a simple form, 13 as shown in equation (3):
Referring to the relevant theories, 22 when calculating the sublimation of outdoor ice sheets or snow fields, we assumed that wind were not coupled with other meteorological factors. So we introduced the influence function g(v) of wind speed on sublimation into equation. For admixtures, we consider by the proportional coefficient k. Therefore, after considering the influence of wind speed and admixture and relationship between vapor pressure and humidity, we obtained the sublimation rate calculation model of pure ice and FRI on the basis of the equations (1) and (2) and equation (3), as shown in equation (4):
The meanings of the related symbols are shown in table below.
The surface size of the specimen is 90 mm × 90 mm and it is made of distilled water or pulp fiber frozen at −15°C. For the measurement of sublimation of outdoor ice cover and snow field, the most commonly adopted methods are evaporation disk method19,20 and snow pile method,21,22 both of which obtain the corresponding sublimation situation through regular thickness measurement or weighing. We used evaporative plate method for sublimation experiment of pure ice and FRI. In the pure ice and FRI sublimation experiment, the mass loss of the specimen was recorded every 3 hours, and the corresponding surface height decline rate was converted. The experimental work conditions are shown in Table 1:
Experimental work condition.
In order to get the influence function g(v) of wind speed, all experiments were carried out in a low-temperature experiment box of −15°C and 70% RH with a size of 1500 mm × 550 mm × 100 mm, as shown in Figure 3, whose wind speed can be adjusted between 0 and 10 m/s. The prepared specimen is placed at the exit of the wind passage so that the wind flows parallel to the upper surface of the specimen.

Schematic diagram of wind speed experimental device.
Results of sublimation experiments
The sublimation rates of pure ice under different wind speeds can be obtained, and the experiment results are plotted as scatter plots and fitted, as shown in Figure 4. It can be seen that the sublimation rate increases linearly with the increase of wind speed, and the sublimation rate of pure ice is obviously greater than that of FRI.

Value of sublimation rates at different wind speeds.
For pure ice and FRI, the linear function of g(v) can be obtained by fitting wind speed and sublimation rate for normalization, as shown in equations (5) and (6):
For the atmosphere where the ice composite shell is located, the air pressure is usually about 1 atm. When the air pressure is about 1.0 atm, the sublimation rate of the pure ice basically remains unchanged, which is about 0.077 mm/d. According U = 70%, Ts = 258.15 K, g(v) = 1.0 and k = 1.0 at this time, we can obatian A0 = 15.321. Due to limitations of experiment results in laboratory, the actual sublimation rate of FRI was corrected by result of outdoor natural environment experiment and k is 0.557 when the material is 2% FRI. 18
Therefore, we can get the sublimation rate equation of ice material considering the wind speed and admixture, as shown in equations (7) and (8):
Numerical simulation of ice composite shell sublimation
Numerical simulation of the effect of solar radiation on sublimation of ice composite shell
In the previous section, we did not study the effect of solar radiation on thesublimation of ice. As a form of energy transfer, solar radiation has an indirect effect on sublimation of FRI. The heat transfer between solar radiation and the ice composite shell increases the surface temperature and exceeds the ambient temperature, thus increasing the sublimation rate. In equations (7) and (8), temperature Ts refers to surface temperature, instead of air temperature. Since values of the two can be regarded as the same at night, Ts can be directly replaced by air temperature. However, during the daytime when the solar radiation is intense, the values of the two are often different, and the difference is even more than 5°C. Therefore, in the actual study of sublimation of ice composite shell, we must consider the temperature rises △T of FRI caused by solar radiation and wind speed, the surface temperature of FRI can be expressed by equation (9):
For the actual ice composite shell structure, the heat transfer process in this area can be simplified as a one-dimensional heat transfer problem along the direction of plate thickness.Using the heat transfer analysis module in ABAQUS to simulate the radiative and convective heat transfer process of 2% FRI in the natural environment, the temperature rises △T can be obtained.
The plate finite element model was established.The plane size of each plate was 1 m. Five kinds of plate thickness were set, namely h = 5, 10, 15, 20, and 25 cm. In Harbin, the instantaneous total horizontal solar radiation intensity usually does not exceed 1400 W/m2. Under the condition that the radiation absorption coefficient of composite ice is 0.25, the actual solar radiation absorbed by FRI does not exceed 350 W/m2 at most. In addition, environmental wind speed usually does not exceed 10 m/s. Eight kinds of solar radiation intensity were set, namely Q = 0, 50, 100, 150, 200, 250, 300, and 350 W/m2. Six kinds wind speed also were set, namely v = 0, 2, 4, 6, and 8 m/s.
The initial temperature field in the analysis was set as −15°C, and the transient temperature field heat transfer simulation was carried out without considering the change of actual temperature with time. In addition, considering that the storage frequency of meteorological data in the meteorological file is usually 1 h, the simulation duration is 1 h. As shown in Figure 5, the temperature distribution of the plate with different thickness after 1 h is given by taking the working condition with the solar radiation intensity Q = 350 W/m2 and wind speed v = 0 m/s as an example. Temperature distributions of different plate thicknesses along section thicknesses are shown in Figure 6. When h ⩾ 10 cm, the difference of surface temperature between different plates is very small, and the trend of temperature distribution along the section plate thickness also tends to be consistent. Therefore, When h ⩾ 10 cm, the temperature rise can be treated as if h = 10 cm.

Result after 1 h heating. (Q = 350 W/m2, v = 0 m/s).

Curve of temperature along section thickness (Q = 350 W/m2, v = 0 m/s).
The temperature rise value △T obtained by numerical simulation was sorted out and fitted, and the temperature rise value △T fitted surface with plate thickness of 5 and 10 cm was finally obtained, as shown in Figure 7.

Fitting surface of △T, Q, and v: (a) h = 5 cm and (b) h = 10 cm.
Extreme Cum surface model is used for fitting, and when the thickness of the plate is 5–10 cm, the surface formula after fitting was simplified as equation (10):
Where k = (0.0510 – 0.00202*h) – (0.0340 – 0.00158*h)*exp(–exp(B)), B = (1.2203+0.0135*h) – (0.6545+0.0008*h)*v
Combined with the equations (8) and (9) and equation (10), the sublimation of FRI can be calculated according to the meteorological data such as the ambient temperature, relative humidity, wind speed and solar radiation intensity. In engineering practice, the calculation equation of 2% FRI sublimation rate on ice composite shell surface is shown in equation (11):
Numerical simulation method of sublimation ice composite shell
Although the sublimation rate of FRI calculation equation S (T, U, Q, v) has been obtained, there are still some difficulties in the practical application of this equation.Using meteorological monitoring equipment such as the weather station, we can get meteorological conditions of a point on the ice composite shell in the environment and the real-time changes.Because the complexity of ice composite shell structure form and structure on the surface of the meteorological conditions of different regions is different, especially the distribution of surface wind speed and solar radiation intensity, the weather stations measurement result can not reflect the actual meteorological environment condition of the structure surface.We also cannot directly calculate the sublimation rate on the surface of the structure and the distribution of accumulation of sublimation by meteorological data from weather stations.
In this section, the sublimation rate of FRI calculation equation S (T, U, Q, v) will be used to establish the numerical simulation method for the sublimation of the surface of the ice composite shell structure. The basic flow of the numerical simulation of the sublimation is shown in Figure 8. According to the meteorological data provided by the weather station, the numerical simulations for solar radiation and CFD flow field were carried out on the surface of the discretized ice composite shell structure respectively to obtain the solar radiation intensity value and wind speed value at each node. Combined with the temperature value and relative humidity value of each time obtained from the weather station, the upward sublimation rate and cumulative sublimation value of each node can be obtained after substituting into S (T, U, Q, v ).

Flow chart of numerical simulation of sublimation of ice composite shell.
The numerical simulation in this section will be implemented in a unified platform Rhino-Grasshopper. On one hand, the spatial inhomogeneity of wind speed and solar radiation intensity on the surface of the ice shell is considered, on the other, it is assumed that air temperature and relative humidity are spatially uniform.
One of the research objectives of this paper is to simulate the real sublimation of the ice composite shell structure as far as possible. Therefore, in the numerical simulation of solar radiation, more attention is paid to the heat gain on the surface of the structure under the real-time sky radiation distribution. So we will adopt the hourly sky brightness distribution model 23 commonly used in the field of building energy consumption simulation. The Ladybug plug-in supported on Grasshopper was used to simulate the distribution of solar radiation on the surface of ice composite shell structure. By reading the horizontal radiation intensity values provided in the weather file and calling RADIANCE, a lighting simulation software with the Perez 24 model was built in and the annual sky dome was generated firstly. Then the sky dome at a specified time period was extracted, and finally the solar radiation projected on the surface of the ice composite shell structure was calculated.
Since 2016, Grasshopper platform began to support the access of Butterfly, a fluid mechanics plug-in, which provided a possibility for CFD simulation by calling computational fluid mechanics software openFoam. The basic flow of the numerical simulation of CFD is shown in Figure 9.

Surface wind speed distribution simulation “battery” connection diagram.
The distribution of solar radiation and wind speed on the surface of the ice composite shell structure at a certain time is simulated. However, the sublimation on the surface of the ice composite shell structure itself is a very slow process, and the sublimation rate varies with the constant change of meteorological conditions. So it is of little significance to obtain only the instantaneous sublimation rate distribution on the surface of the ice composite shell structure. The ultimate goal should be to obtain the distribution of accumulated sublimation on the surface of the ice shell structure in a specified period of time. Using Python programming language, new batteries were developed on Grasshopper platform to calculate accumulated sublimation.
Effects of sublimation on an ice shell
Numerical simulation results of sublimation of ice composite shell
Numerical simulation of Sublimation of a ice composite spherical shell structure with 10 m span, 5 m height and an intial thickness of 5 cm was performed. The typical annual meteorological data obtained from the actual meteorological data in Harbin area from 1957 to 2017 were used for the meteorological data of numerical simulation. According to the service time of the ice composite shell in the actual engineering, the sublimation of ice composite spherical shell from January 1st to March 12th for 70 days is simulated by using meteorological data. Considering that a CFD simulation consumes a lot of computation time, this paper adopts some simplified processing in CFD simulation. We consider the change of wind speed without wind direction angle,and CFD simulation was only carried out for the main wind Angle in the above research period.
The ice composite spherical shell is simulated with solar radiation and wind speed distribution and the final cumulative sublimation distribution, as shown in Figure 10.

(a) Distribution of solar radiation intensity on the surface of the ice composite spherical shell (12:00 to 13:00 on January 1st), (b) distribution of wind speed on the surface of the ice composite spherical shell (12:00 to 13:00 on January 1st), and (c) distribution of cumulative sublimation.
As shown in Figure 11, the initial weight of the spherical shell was 7052.72 kg, and the rate of weight decline in the first four 10-days was not significantly different, and the rate of weight decline in the last three 10-days was significantly accelerated. The final weight of the spherical shell was 3577.14 kg, and further obtained that the total weight loss rate caused by sublimation was 49.3%, with an average daily weight loss of 0.68%. During sublimation, the thickness of the ice composite shell decreases by 0.38 mm every ten days.

Curve of structural weight change.
Along the east-west symmetry axis of the ice composite ice structure, the thickness variation curve of 5 points at the height of 2, 4, and 5 m is given successively from south to north, as shown in Figure 12. We can find that the sublimation rate of different locations of the ice shell is obviously different, and the sublimation rate on the sunny side was 1.7 times that on the shady side.The area with the most serious sublimation of the ice shell was located at the top of the structure, with the sublimation amount reaching 41.3 mm, the thickness of the structure decreased from the initial 5 cm to less than 0.9 cm.

Curve of structural thickness.
Changes of static behavior of ice composite shell during sublimation
Two finite element models, the initial spherical shell and the sublimated spherical shell, were respectively analyzed statically under normal load in ABAQUS. Based on the experience of practical projects, the structural analysis is carried out considering the dead weight of the structure. According to the existing 2% FRI material test results, 3 the material constitutive model used in this section is assumed to be linear elastic. Fixed constraints are adopted at the bottom of the ice composite shell.
The curves of the maximum tensile stress at the lower part of the waist, the maximum compressive stress at the top and the vertical displacement at the top over time is shown in Figure 13. It can be seen that the three curves decrease gradually. The stress and the vertical displacement of the ice composite shell decrease due to the decrease of the weight with the sublimation of the ice composite shell. However, the structure changed from regular shape to irregular shape, which may be increase the sensitivity of the structure to external loads.

Structural internal forces and deformation changes.
The ice composite shell structure is a kind of defect sensitive structure. In the actual construction, due to the deformation error of the formwork and the freezing speed, there will be some initial defects in the shape and interior of the structure, which will greatly reduce the stable bearing capacity of the structure. Therefore, a nonlinear buckling analysis is carried out considering the influence of defects.The uniform defect mode method 25 was adopted to impose geometric defects, and the first buckling mode was taken as the initial defect shape. The maximum defect value was 1/300 of the ice conposite shell structure span, about 33 mm.
According to the results of nonlinear buckling analysis, the nonlinear instability failure modes of ice composite spherical shell at the initial stage and 70 days after sublimation are shown in Figure 14. The failure mode of ice composite spherical shell before and after sublimation of FRI is the overall instability caused by local buckling. The initial failure mode of ice composite spherical shell is located at the waist of the northeast side of the structure, while the failure mode of sublimated ice composite spherical shell transfers to the area with the thinniest thickness at the top of the structure.

Failure modes of structural buckling: (a) initial ice composite shell and (b) sublimated ice composite shell.
As shown in Figure 15, it is the change of structural stability bearing capacity with the progress of sublimation. It can be seen that the stable bearing capacity of the structure decreases gradually due to sublimation of FRI. After 70 days of sublimation, the bearing capacity of the structure was reduced by a total of 95.8%. Therefore, the structural thickness thinning and uneven spatial distribution caused by sublimation of FRI have a serious effect on its stable bearing capacity.

Curve of bearing capacity with consideration of sublimation.
Conclusion
In this paper, based on the existing theory of pure ice sublimation rate calculation, the sublimation tests of pure ice and FRI were designed and carried out successively. The influence of air velocity and admixtures on sublimation of ice was studied by single factor control experiments in the laboratory and the calculation equation of different ice sublimation rates were given. Combining numerical simulations of the solar radiation and CFD, a numerical simulation method for the surface sublimation of composite ice shell structure was developed. A finite element method for reconstruction of composite ice shell structure with sublimation was presented, and a specific ice shell was analyzed. The following conclusions can be drawn: (1) The increase of wind speed generally increases the sublimation rate of ice. Under the same wind speed, the sublimation rate of pure ice is obviously higher than that of composite ice. (2) During sublimation, the thickness of the ice composite shell decreases by 0.38 mm every 10 days and the sublimation rate on the sunny side was 1.7 times that on the shady side. (3) After 70 days of sublimation, the thickness of the ice composite spherical shell becomes thinner and uneven, which increase the sensitivity of the structure to external loads.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors are grateful for the financial support from National Natural Science Foundation of China (NSFC) under Grant No. 51778182.
