Numerical Simulation Research on Precooling Characteristics of LNG Discharge Pipeline

The low-temperature fluid flowing into the LNG unloading pipeline will cause violent heat exchange phenomenon, resulting in damage to the pipeline and affecting normal operation. Therefore, pre-cooling treatment is required in advance. Considering the fluid-solid heat transfer and fluid flow characteristics comprehensively, multi-physics coupling simulation study on the pre-cooling problem of the LNG discharge pipe with a π-shaped elbow is carried out. The analysis from BOG and LNG pre-cooling processes respectively shows that: the fluid with low temperature and high density is mainly concentrated at the bottom of the pipe and outside the bend, the wall temperature has a sudden drop at the bend, and there is a large temperature difference between the gas side and the liquid sidewall surface. Mixed convection heat exchange mainly occurs in the tube. There will be a small area of natural convection-dominated heat exchange state in the elbow and the downstream pipe section.


INTRODUCTION
In order to meet the increasing energy demand, look for fuels that can replace hydrocarbons such as diesel and gasoline [1], which produce little or no greenhouse gases. while natural gas meets the current needs, it is clean, renewable and mass production, etc. The advantage is that it generates less greenhouse gases, can exist in a variety of physical forms and is more widely used in various environments. When natural gas is cooled to -162℃, it will be converted into liquefied natural gas [2]. When converted to liquid, its volume was reduced by more than 600 times [3]. LNG is not only the cleanest form of natural gas; it also provides a safer and more economical option for transportation and also increases its storage capacity [4][5][6]. Therefore, in the current energy market, LNG has been widely used around the world, and now is the right time to review the development status of liquefied natural gas (LNG) [7][8][9].
Domestic scholar Lu Chao analyzed the BOG+LNG pre-cooling method from the BOG pre-cooling process and the LNG pre-cooling process, and established a one-dimensional flow heat transfer model to analyze the changes of various parameters of the pre-cooling process [10]; Ji Junyi et al. studied the nitrogen pre-cooling scheme, and simulated the precooling process from theory, numerical simulation and example verification [11]; Domestic scholar Chen Taoqiang established gas-liquid stratified flow model and dynamic simulation operation platform for the LNG discharge pipe pre-cooling process, and not only optimized the BOG+LNG pre-cooling method, but also proposed three typical pre-cooling schemes and used to establish Platform's simulation calculation and optimization [12].
The earliest experimental research was carried out by foreign scholars Burke [13], Graham [14], Bronson [15], Chi and Vetere [16], Steward [17] and others. Burke et al. based on the assumption of one-dimensional heat transfer through the tube wall and the infinite heat transfer rate from the cryogenic fluid to the tube wall, and established a rough cooling model, but the model ignores the effect of the flow method on the heat transfer rate [13]. Graham et al. used their experimental data to correlate the heat transfer coefficient and pressure drop with the Martinelli number [14]. Chi and Vetere studied the flow and heat transfer conditions in the horizontal tube during low temperature cooling. According to the temperature and volume fraction of the fluid in the cooling process, information about the flow pattern is obtained [16]. Steward established a uniform flow model for cryogenic cooling. The model treats cryogenic fluids as a homogeneous mixture. Solve the continuity, momentum, and energy equations of the mixture to obtain the density, pressure, and temperature of the mixture. Later [17].Jun Liao established a quasi-steady-state model to predict the temperature change process of the tube wall in the low-temperature fluid horizontal pipeline, and combined with the layered lowtemperature cooling flow structure to develop a new film boiling heat transfer model [18]. Reid Shaeffer explored the cooling method of low-temperature pipelines and studied pulse-flow cooling pipes. A comparison of heat transfer characteristics between continuous flow and pulsed flow modes is proposed, and the influence of flow oscillation on heat transfer is revealed [19]. Yuan [20] scholar used experiments to conduct experiments on cryogenic fluid cooling equipment, and found that the data curves of the cryogenic cooling experiment and the boiling experiment have certain similarities, but the difference is that the test part did not supplement the external energy to cause the temperature to drop. The final test data shows that cryogenic liquids will have three boiling regions during the pipeline cooling period: thin film, transition and nucleate boiling. Steward [21] has also carried out related experiments, using liquid nitrogen and liquid hydrogen as working fluids. The study found that the liquid in the pipeline inlet is too cold, the transmission pipeline is long, the liquid density is high, and the inlet valve is opened quickly, which can cause or aggravate the surge phenomenon in the pipeline cooling process. Klimenko [22] scholars studied the effect of pipe orientation and geometry on the heat transfer of nitrogen in a two-phase forced flow. Using the change of Froude number, it was concluded that the effect of gravity on heat transfer in vertical pipes and horizontal pipes is in the horizontal pipe section. It is mainly manifested in the increase of heat transfer in the upper part of the channel as the Froude number decreases. The reason is that this leads to a significant increase in the average cross-section of the heat transfer coefficient in the nucleated boiling region (20-30%).
Liquefied natural gas (LNG) is a low-temperature cryogenic medium of -162℃, and the design temperature range of its unloading pipeline is -170℃~60℃. In order to prevent the safety problems caused by the rapid thermal changes of the pipe wall temperature, it is particularly important to choose a reasonable pre-cooling medium and control the temperature drop rate. In order to ensure the safe operation of the LNG receiving station, the numerical simulation method is used to analyze the BOG pre-cooling method of the short horizontal section and the expansion elbow, and the process of LNG cooling pipeline is simulated. LNG and BOG are complex mixtures, however, their methane content is above 90%, the main properties and physical properties are similar to methane, so in the simulation of pre-cooling research, methane is used instead of LNG for simulation; The material properties of the insulation layer are simplified to fixed parameters that do not change with temperature.

Model assumption
In the actual LNG discharge pipe pre-cooling process, in order to ensure the safe and accurate operation of pre-cooling, a variety of tasks are required. At the same time, heat leakage will occur in the pipeline pre-cooling process, such as valves and connections. Since this article only conducts numerical simulation studies and considers it to be an ideal state, the following assumptions are made about the model: (1) The pipeline has been purged; (2) No cooling loss during the pre-cooling process; (3) The physical properties of the pre-cooling medium only change with temperature.

Governing equation
According to the assumptions, in the pre-cooling numerical simulation, the model needs to satisfy the control equations necessary for fluid flow and heat transfer. Following the basic conservation laws, the corresponding three-dimensional control equations are established as follows: (1) Continuity equation Given in Equation (1), where, is density, is time, and u is velocity vector, , and are the components of velocity vector in the , , and directions, respectively.
Since the LNG pre-cooling process is a variable-density flow process, parameters such as density viscosity change due to temperature changes, so equation (1) is used.
(2) Momentum equation Given in Equations (2) and (3), where, is static pressure, is the viscous stress acting on the surface of the microelement caused by the molecular viscosity, and is gravity volume force in i direction, is external volume force in the i direction.( , is research dimension, In the two-dimensional model, , = 1, 2 and in the three-dimensional model, , = 1, 2, 3). (4) and (5),

(3) Energy equation Given in Equations
where is viscous force tensor, is strain rate tensor, and is heat source.
In the COMSOL software simulation, ignoring the viscous effect caused by the fluid flow will cause the fluid to heat, and the fluid density will change according to temperature. Therefore, the pre-cooling model belongs to the convection heat transfer model, where the energy equation can be expressed by temperature: In this paper, the COMSOL finite element software is used to solve the two-phase flow equations, fluid and solid heat transfer equations. The pre-cooling method of the LNG discharge main is BOG+LNG pre-cooling method, and there is a coexistence state of gas and liquid between the two stages. In order to accurately simulate the flow heat transfer process at this stage, it is necessary to distinguish between the gas phase and the liquid phase. For the pre-cooling simulation of the LNG discharge pipe, the LNG filling process is a gas-liquid twophase medium, and there is obvious gas-liquid stratification Interfaces, fluids are not miscible, the choice of interface tracking method is more suitable for this simulation. The phase field method with higher applicable precision is used to trace the interface of the gas-liquid twophase precooling process.

PHYSICAL MODEL AND MESHING
3.1. Three-dimensional discharge pipeline model In this paper, the numerical simulation analysis of the pre-cooling of the discharge main pipe in the pre-cooling system is carried out. The physical model of the discharge elbow is shown in Figure 1. Most of the actual LNG discharge pipelines on the site are of the 100-meter level or even the kilometer level, and the pipe wall thickness is relatively few orders of magnitude. In the discharge main pipe, the expansion elbow will be installed in the pipeline every 100m to reduce the deformation caused by thermal stress, so this article only simulates the unloading main pipe containing the π-shaped elbow part, the model is established and the mesh is divided. The horizontal pipe section in front of the bend pipe entrance is the distance from the outlet of the previous compensator to the entrance of the next compensator, and the horizontal pipe section at the outlet of the bend pipe is 5m. The specific parameters of the model are shown in Table 1.

Three-dimensional model of conjugate heat transfer grid
The mesh division of the LNG discharge pipe will adopt different division methods according to the calculation domain. The pipe wall adopts structured mesh with a faster calculation speed and easy convergence because it only involves temperature field changes. For the fluid calculation domain in the tube, there are changes in speed and material properties at the same time as temperature changes are involved, so tetrahedral mesh with higher accuracy is used for segmentation. The meshing of the 3D model is shown in Figure 2.
In this simulation, since the solid domain is the pipe wall and the insulation layer, it has the characteristics of regular shape. At the same time, only the temperature field is calculated in the solid domain. The grid of the solid domain adopts the calibration method of the general physical model, and the structure of four iterations is used. The grid divides the model, the maximum grid division is 0.03cm, and the maximum processing depth is 4. For the construction of fluid basin grids, unstructured grids are used to ensure the continuous and effective flow of fluids due to changes in density, viscosity, heat transfer parameters, and flow parameters between fluids. At the same time, the boundary layer is divided in the fluid domain. Between the boundary layer and the tetrahedral mesh in the fluid domain, a pentahedral cone mesh is used for transition. The solid and fluid domain meshing is shown in Figure 3.

BOG AND LNG PRE-COOLING CHARACTERISTICS ANALYSIS
In this paper, BOG+LNG is used to simulate the pre-cooling process of the discharge short horizontal section and the expansion elbow. First, BOG gas is used to pre-cool the pipeline to reduce the pipeline temperature to -120 ℃, then slowly access LNG liquid to deep cool the pipeline, and finally the temperature of the pipeline drops to -150℃, and the pre-cooling ends.
BOG pre-cooling is divided into 5 stages. Pre-cooling due to the longer pre-cooling time, the tube wall temperature is characterized every few hours. The initial temperature of the LNG stage has reached -120 ℃, so the pre-cooling method in sections is no longer used. At the same time, the same speed (0.6m/s) precooling is adopted for the two-stage precooling. This article uses segmented pre-cooling parameters as shown in Table 2:

Wall temperature analysis of discharge pipe
Pre-cool the LNG discharge elbows connected by four 90° elbows in sequence, and perform stepped pipeline pre-cooling on the BOG pre-cooling stage. The total calculation time is 28 hours. Because the pre-cooling medium directly contacts the pipe wall pipe material, the outer wall insulation layer is hidden in the analysis of the pipe wall temperature result. In order to facilitate the analysis of the pipeline, the straight pipe sections are defined as Z1, Z2, Z3, Z4 and Z5 along the fluid flow direction; the 90° elbows are defined as W1, W2, W3 and W4, as shown in Figure 4. As shown in Figure 5, the wall temperature rises along the direction of fluid flow, and the temperature of the pipe wall is unevenly distributed. The main performance is that the temperature of the upper wall surface of the straight pipe section is higher than that of the lower wall surface. The main reason is that the thermodynamic properties of the BOG gas change with the change of temperature, resulting in fluid stratification. The temperature of the pipe wall is unevenly distributed in the 90° elbow and the straight pipe section near the downstream. The temperature of the inner wall surface of the bend pipe is higher than that of the outside, and in the downstream straight pipe section, the wall temperature in the same direction as the inside of the upstream curve is higher than the wall temperature in the same direction on the outside of the upstream curve, and the low temperature area of the outer wall surface temperature is significantly longer than that of the inside. According to the figure, the temperature distribution of the LNG wall surface is the same as that of the BOG pre-cooling stage, and the temperature of the pipe wall gradually increases with the direction of fluid flow, and there is a sudden temperature drop at the bend, however, compared with the BOG precooling stage, the overall temperature is more uniform.  Figure 6 shows the wall temperature curve of BOG+LNG at each pre-cooling stage .With the progress of different BOG pre-cooling stages, it was found that the temperature distribution in each stage was the same, with the lowest temperature at the pipe inlet and the highest temperature at the pipe outlet. Lateral comparison of the various pre-cooling stages reveals that as the pre-cooling time progresses, the overall temperature of the pipeline exhibits a temperature drop at different rates. The wall temperature decreases rapidly during the beginning of each stage of pre-cooling, and approaches the end of each stage of pre-cooling, the temperature drops gradually. There are obvious fluctuations in the temperature drop rate curve of the pipe wall at the bend, and a sudden change in temperature drop occurs at or near the entrance of the bend. When the pre-cooling time is close to the end of each stage, the temperature at each position of the pipe wall is basically maintained consistent, infinitely close to the temperature of the incoming pre-cooling fluid to meet the stage pre-cooling requirements.
In the LNG pre-cooling stage, the temperature distribution on the wall surface is more uniform, and the temperature difference between the inlet and outlet is small, But its temperature drop rate is too fast, due to the small fluctuation of LNG temperature. In order to slow down the rate, the fluid flow should be reduced to achieve the purpose of safe precooling. In the practical application of the LNG receiving station project, it is considered that the pre-cooling ends when the temperature at the outlet of the pipeline reaches the basic requirements. According to the data curve compiled by this article, the method can be judged to be effective.  Figure 7 shows the temperature cross-sectional diagram of the LNG discharge pipeline in the x direction at different pre-cooling stages. Looking at the overall temperature change of the pipeline from the figure, the low-temperature pre-cooling medium enters the discharge pipe from the left Z1 pipe and is heated, and it can be seen that as the time course increases, the pre-cooling medium flowing through the Z1 pipe follows the direction of fluid flow. The temperature gradually decreases. At the same time, analysis of W1 shows that the pre-cooling medium once enters the elbow, the temperature distribution of the fluid will be reconstructed. And it can be clearly seen that no matter how the orientation of the elbow changes, the stratification of the straight pipe sections in and downstream of the elbow still exists. The fluid temperature is stratified according to the bending center of the elbow to the outside of the elbow, and high-temperature fluid and low-temperature fluid are distributed from the inside of the elbow to the outside of the elbow.

Pre-cooling medium temperature analysis
This stratification phenomenon will affect the downstream pipeline. By observing the W2 pipeline, it is found that the cryogenic fluid is mainly concentrated on the inside of the pipeline. This is because the upstream fluid is stratified, and the pipeline Z3 downstream of W2 can be seen that the cryogenic fluid still flows outward. The effect of fluid stratification caused by the bend pipe continues to the Z5 pipeline, but after the pre-cooling medium passes through multiple bend pipes with different bending centers, the fluid temperature stratification phenomenon is gradually weakened, mainly due to the fast flow velocity of the fluid and the short distance of the bent pipe section, the temperature distribution of the fluid is gradually disturbed after being repeatedly disturbed.  Figure 8 shows the temperature distribution of the fluid medium after LNG pre-cooling for 2 hours. According to the figure, the temperature of the fluid in the tube gradually increases with the direction of the fluid flow, but due to the larger specific heat capacity of the fluid itself, the temperature of the fluid in the tube is retained, the cooling loss is smaller, and the overall temperature difference is smaller than the BOG stage.  Figure 9 shows the longitudinal temperature cross-sectional distribution along the direction of fluid flow. The temperature cross-section is obtained for each pipe section. Due to the difference between the wall surface and the fluid temperature, the internal fluid temperature gradually increases with the flow direction. In the Z1 pipeline, the fluid temperature mainly shows a gradual decrease in the temperature in the direction of gravity. This is due to the phenomenon caused by the dual action of gravity and buoyancy. After the fluid exchanges heat with the tube wall, the temperature rises and the density decreases, and the stratification phenomenon appears due to the influence of gravity. The pre-cooling medium with higher temperature is concentrated on the top of the pipe, and the pre-cooling medium with lower temperature shows the opposite trend, concentrated on the bottom of the pipe.
The overall temperature gradient of the pipe section shows an increasing trend along the way, and the temperature gradient at the outlet reaches the maximum. This is due to the temperature and velocity gradually increase as the fluid moves forward and the existence of the buoyancy force, which make the fluid on the lower wall absorb heat. The velocity is obviously greater than the fluid on the upper wall surface, which enhances the mixing of cold and hot fluids in the flow field. When the fluid flows to the bend section or even the downstream section, the temperature distribution gradually shifts. Before the W3 pipeline, it was greatly affected by gravity. When the fluid completely passed through the W3, the temperature gradient of the downstream fluid gradually distributed mainly in the direction of centripetal force. At the same time, due to the multi-section elbow, the temperature mixing is more intense, and the temperature distribution at section 10 reaches the most uniform state. Observation of sections 3, 5, 7, 9 shows that the secondary flow at the section disturbs the cold fluid from its initial position to deviate from the flow in the bend, The heated fluid near the pipe wall is brought to the inside of the pipe following the circulation, and at the same time, it is extrapolated to the outside of the elbow due to the appearance of the secondary flow, and the low-temperature fluid is brought to the inside of the pipe by the secondary flow. The repeated flow phenomenon caused the thermal boundary layer at the bend to gradually decrease, and this caused the sudden temperature drop of the bend and the nearby straight pipe.  Figure 10 shows the temperature distribution of the Z1 section of a straight pipe section. According to the figure, the stratification trend still exists in the temperature of the LNG stage, and the analysis of the BOG pre-cooling stage found that the stratification of the fluid in the pipe cannot be ignored, and the fluid property changes caused by the temperature are more obvious. The formation of vortices accelerates the cooling of the pipeline. At the same time, the change of the attribute also causes the formation of vortex in the straight pipe section, and accelerates the cooling of the pipe.

Analysis of heat transfer performance
According to the above analysis, due to the heat exchange in the tube, the change of the thermal physical properties of the fluid leads to the change of the flow speed and the secondary heat exchange caused by the buoyancy. In order to clarify the heat exchange situation inside the tube, the Reynolds number , Prandt number , Nusselt number and the Graffian number that caused the buoyancy will be studied along the pipeline. In order to investigate the influence of natural convection on the fluid flow and heat transfer in large-diameter pipes, the ratio of buoyancy to inertial force * (Equation (9)) is used, where the magnitude of this value can determine the intensity of natural convection in the heat exchange process.
* judges the effect of natural convection in the heat transfer process in a detailed table as shown in Table 3. The specific formula is as follows in Equation (6-9): Figure 11 shows the distribution of Reynolds numbers along the pipeline. It can be seen from the figure that the distribution of Reynolds number in the pipeline shows the same trend as the distribution of the fluid velocity field. As the temperature of the fluid increases in the Z1 pipe section, its dynamic viscosity and the flow velocity gradually increases, resulting the Reynolds number increases with the direction of fluid flow in the Z1 pipeline, which strengthens the turbulence intensity of the fluid and improves convective heat transfer.
In the W1 pipe section and its downstream pipelines, affected by the bend, the number is larger at the center of the bend. At the same time, the downstream pipeline is also greatly affected by its size distribution. Its basic trend is similar to the speed distribution, and will not be elaborated too much. Table 3. List of Gr*judging the intensity of natural convection in heat exchange Gr*≤0.1 Forced convection plays a leading role in the heat transfer process, it is considered that the effect of natural convection on heat transfer can be ignored

0.1< Gr*<10
Natural convection and forced convection are equally important in the heat exchange process, and the heat exchange is considered to be mixed convection at this time 10≤Gr* Natural convection plays a dominant role in the heat transfer process, and the effect of forced convection on heat transfer can be ignored Figure 11. Distribution Figure 12 shows the Prandtl number distribution along the pipeline. As shown in the figure, along the fluid flow direction, the Prandtl number tends to decrease gradually. The size of Prandtl number represents the ratio of the velocity boundary layer and the temperature boundary layer. When the value decreases, it is proved that the increase in the temperature boundary layer causes the degree of heat transfer. The Prandtl number reached the maximum at the entrance of the pipeline, which proves that the temperature boundary layer at this point has the strongest heat transfer relative to the minimum, which leads to a tendency for the wall temperature to increase in the direction of fluid flow. Through integral integration of the Prandtl number, it is found that the value greater than 0.7 conforms to the applicable range of Nusselt's formula.  Figure 13 shows the distribution diagram of the Grashof number along the pipeline. As shown in the figure, along the direction of fluid flow, the Grashof number tends to decrease gradually. The size of the Grashof number represents the ratio of buoyancy force to viscous force. The heat exchange intensity in the tube gradually decreases with the direction of fluid flow. The Grashof number at the bend is mainly distributed on the outside of the pipe, and is affected by the secondary flow at the bend, and the heat transfer is mainly affected by the bend. Figure 13. Distribution Figure 14 shows the distribution of Nusselt number along the pipeline. According to the formula of Nusselt number, Nusselt number and Reynolds number have a clear correlation, so they have the same trend. It can be seen that there is an obvious heat transfer enhancement phenomenon at the pipe pre-cooling elbow. The section 1 (at Z1) and section 3 (at W1) were averaged, the Nusselt number at section 1 was 106, and the Nusselt number at section 3 is 110, which is 1.03 times. In the elbow, the Nusselt number is mainly concentrated near the center of the elbow, but it can be clearly seen that after passing the elbow, there is an obvious low value area of the Nusselt number. The existence of this area will make the heat exchange of the pipe uneven, resulting in inconsistent temperature drop rate of the pipe wall, and thermal stress caused by the temperature difference between the upper and the lower. For the LNG full-flow state, the Nusselt number distribution law is affected by the tube type and its value is larger compared to the BOG stage. Regardless of the BOG or LNG pre-cooling stage, the Nusselt number distribution is the same.  Figure 15 shows the distribution diagram of * number. According to the calculation, the minimum value of * in the pipeline is 0.5, which is much greater than 0.1 in the discriminant, and the color is displayed as black when the * value is greater than 10 in the figure. On a non-negligible order of magnitude, the influence of buoyancy will dominate especially in some areas. According to the provisions of the discriminant, the heat exchange in the tube is mainly mixed convection heat exchange. In some areas, In some areas, due to the geometry of the elbow, the centripetal force causes the fluid to concentrate on the wall of the pipe, which reduces its natural convection effect. However, in the opposite area, because there is less fluid flowing in the black part of the figure, the Reynolds number and Nusselt number are lower and the heat transfer strength is reduced. Therefore, the heat transfer in the black area is mainly natural convection.

CONCLUSIONS
The numerical simulation method was used to analyze the flow and heat transfer of the whole process of LNG discharge pipe pre-cooling, and the following conclusions were drawn: (1) During the BOG pre-cooling process, the wall surface temperature rises along the direction of fluid flow, and there is a phenomenon of a sudden temperature drop of 3K at the bend. The inner wall surface temperature at the bend is higher than the outer wall surface temperature. The upstream pipeline has vigorous heat exchange and a rapid temperature drop rate. As the fluid flows downstream, the cooling capacity of the precooling medium decreases, the heat transfer slows down, and the temperature drop rate slows down; (2) During the LNG pre-cooling process, the parameters in the full-flow stage are the same as the BOG pre-cooling stage. For the non-full flow stage, due to the low liquid phase temperature and high density, the temperature of the liquid sidewall surface is much lower than the temperature of the gas sidewall surface (temperature difference is about 30K), and the temperature drop rate of the liquid sidewall surface is higher than that of the gas sidewall surface. The presence of the expansion elbow causes a significant disturbance of the internal fluid, which significantly reduces the temperature difference between the elbow and the downstream pipe wall; (3) During the overall pre-cooling process, the fluid in the tube is stratified. Under the action of gravity, the pre-cooling medium with low temperature and high density is concentrated at the bottom of the pipe; however, at the elbow, the pre-cooling medium is concentrated outside the elbow due to the presence of centrifugal force. Judgment and calculation of the internal fluid's buoyancy, and find that mixed convection heat exchange mainly occurs in the tube, and the buoyancy and its natural convection cannot be ignored; (4) The flow heat transfer in the tube is analyzed, and the Prandtl number is calculated to float up and down at 0.7, which is in line with the calculation range of the Nusselt number. Draw the cloud map of the Nusselt number distribution in the pipe. At the lower temperature of the pipe wall, there is a large Nusselt number phenomenon. At the same time, the geometric structure of the pipe has a great influence on the flow heat transfer. According to the analysis of buoyancy, it is found that mixed convection heat exchange mainly occurs in the tube, and there will be a small area of natural convection-dominated heat transfer state in the bend outlet and downstream pipe section of the pipe.
of SW-SAGD with between Heel and Toe Injection in heavy oil reservoirs and University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province (No. UNPYSCT-2017035).