Cryogenic liquid nitrogen fracturing is expected to provide an effective stimulation method for hot dry rock reservoirs to increase heat production. This paper establishes a three-dimensional model to calculate the distributions of temperature and stress of the reservoir rock when liquid nitrogen is injected into the wellbore. The sensitivity of different parameters and water fracturing to the stress state is studied. The results indicate that when liquid nitrogen is injected into the bottom of well, a huge heat exchange occurs on the rock surface, which generates great thermal stress on the fluid-solid interface, and the value of thermal stress exceeds the tensile strength of rock. For the effect of parameters, the primitive temperature of the rock has a significant impact on the value of maximum principal stress. The pressure drop and ambient pressure affect the thermal stress slightly. At the same time, a series of experiments are conducted to validate the effect of thermal stress induced by liquid nitrogen injection on the rock fracture. As the temperature rises, the shale samples are broken more severely at the action of thermal stress. Thus, the study of liquid nitrogen fracturing provides a scientific and effective method for geothermal exploitation.

1. Introduction

The hot dry rock (HDR) is green, low-carbon renewable energy in the reservoirs, which is regarded as an important alternative for conventional energy [1, 2]. HDR is defined as the high-temperature rock between 150 and 650°C in the deep subsurface. The common HDR includes granite, diorite, and gneiss [3]. The progress of geothermal mining from HDR is injecting cold water into reservoir and then returns to the ground through the production well. The cold water is heated by contacting the hot rock to achieve the purpose of extracting thermal energy. However, it is difficult to rely on the conductivity of natural fractures and the seepage ability of rocks to obtain economical heat flow due to the extremely low permeability of HDR reservoirs [4]. Because of extracting thermal energy, it not only needs to increase the permeability of the reservoir, but also needs to increase the surface area of the fracture as much as possible. Therefore, stimulation methods should be adapted to increase the permeability [57].

At present, the main stimulation method of HDR resources is hydraulic fracturing, which has been widely used in increasing production of unconventional oil and gas [811]. However, the traditional hydraulic fracturing has some problems, such as the high fracture pressure, pore blockage caused by water sensitivity damage, and large water consumption [1214]. In order to solve these problems, waterless stimulation methods are proposed. Compared to traditional hydraulic fracturing techniques, waterless fracturing is mainly unconstrained by water resources, which not only increases the fracturing fluid backdraft rate and reduces the damage of fracturing fluid to high-temperature reservoirs, but also reduces the waste of water resources.

As a method of waterless fracturing, cryogenic liquid nitrogen (LN) fracturing has its unique advantages in the process of fracturing HDR because of its extremely cryogenic property [15]. Cryogenic fracturing is defined as follows: huge temperature gradient causes local tensile stress to break rock, which is generated by the contact between cryogenic fluid and high-temperature rock [16]. The rock breakage due to thermal stress results from two aspects: thermal stress generated by the temperature gradient and adjacent minerals with different thermal expansion [17]. Thermal stress will cause microscopic damage inside the rock and increase the permeability, which is beneficial to extract heat from reservoir [18]. At the twentieth century, LN fracturing stimulation technology was used to increase production in the low-permeability Devonian shale, and favorable stimulation effects have been obtained [19]. For this technology, LN is injected into the target reservoir as a fracturing fluid instead of water-based fracturing fluid. In this case, thermal stress is generated to fracture rock due to cryogenic LN contact with high-temperature rock [20, 21]. During the LN fracturing, not only the macrofracture occurs in the rock, but also the intergranular breakage is generated in the microstructure, which promotes the appearance of complex fracture networks [22, 23].

Many previous investigations have been performed to study the rock damage due to LN cooling. Ren et al. [24] carried out experiments of LN thermal shock on coal sample and found that both the wave velocity and amplitude reduced sharply after LN treatment. The experimental results indicated that the internal structure damage occurred in the coal, which is beneficial for generating a fracture network to increase the permeability of rock. Cha et al. [16, 25] conducted cryogenic fracturing studies by injecting LN into concrete and sandstone samples. In this study, numerous microcracks were generated around the borehole due to the effect of thermal stress. This indicated that the LN entering the borehole can change the rock structure and induce fracture creation. Cai et al. [26, 27] investigated the change in the mechanical properties of rock samples after LN treatment. According to the experimental study, the obvious thermal fracture appeared on the rock surface, and the tensile strength and compressive strength of the rock were obviously weakened. This indicated that the treatment with LN has a great influence on the strength of rock samples. Yang et al. [28] conducted granite breakage experiment under different rock temperature and confining pressure. The experimental results showed that the fracture initiation pressure of granite was reduced and the microscopic pore structure was expanded after LN treatment. For the study of numerical simulation method, Yao et al. [29] simulated the LN stimulation process by the TOUGH2-EGS. The simulation results were also compared with experimental data. They found that the increase in reservoir seepage can be achieved by the injection of cryogenic fluid with high pressure. Cai et al. [30] simulated the transient fluid flow and heat transfer in the process of LN being injected into the downhole. The influence factors were also analyzed in this study. The results indicated that cryogenic LN can be generated in the downhole and has a great impact on reservoir temperature. Kim et al. [21] studied the effect of rapid cooling on granite samples and found that the crack growth occurred in granite samples. Furthermore, the thermal stress distribution of rapid cooling was illustrated by a transient thermodynamic model. Zhang et al. [31] studied the granite damage under the impact of cryogenic LN by applying SEM methods. They found that the number of cracks and complexity of fracture network increased with the increase of temperature gradient. This indicated that the rock structure was damaged by LN and the temperature gradient was an important factor. Wu et al. [32, 33] conducted a set of experiments to investigate variations of physical and mechanical properties of granite, including density, wave velocity, strength, and elastic modulus. According to the experimental results, they pointed out that the destruction of granite is aggravated with the number of heating and cooling cycles increasing. Xi et al. [34] performed an experimental study into the effect of water cooling on mechanical properties of granite samples. The influence factor of rock initial temperature was considered in this study. They found that the thermal cracking occurred inside the rock, leading to the degradation of mechanical properties of granite samples. Huang et al. [22] conducted experiments to investigate the characteristics and mechanism of granite samples with LN jet. The results showed the damage of granite samples and the corresponding mechanical properties deteriorated obviously with the growth of primitive temperature. This suggested that the primitive rock temperature is an important factor determining the performances of rock breaking and cracking with cryogenic LN.

According to the literature review, most of studies focused on the mechanical properties deterioration of granite caused by thermal stress after LN treatment. The fracture characteristic of granite attributed to LN was discussed in laboratory investigation. However, the effect after injecting LN into the bottom of reservoir is still not clear. Hence, before the LN can be extensively implemented as a fracturing fluid instead of water-based fracturing fluid, many fundamental problems should be addressed to provide guidance on HDR stimulation. These problems included the heat transfer and the distribution of stress at the bottom of well during injecting LN into wellbore. The heat transfer contributes to the change of fracturing conditions, and the distribution of stress is directly related to the effect of fracturing.

This paper conducted a numerical simulation study on the heat transfer and the stress distributions of reservoir when LN was injected into the HDR reservoir. During the process of LN flowing into the wellbore, the tensile stress was generated due to the thermal stress and fluid pressure. The maximum principal stress was selected for analyzing the effect of thermal stress on rock failure, because the maximum principal stress can effectively determine the rock fracture when the tensile failure was generated in rock. The computational fluid dynamics (CFD) was used in the fluid region with standard k-ɛ turbulence model. The fluid-solid interface was adapted to conjugate heat transfer method. In the simulation, granite was adopted in solid region, which is common in HDR. Due to the large elastic modulus and small deformation in granite, the thermoelastic mechanical model was adapted to calculate the stress and strain in solid region.

2. Geometric Model and Parameters

2.1. Physical Model

A three-dimensional model is established to study the heat transfer between LN and high-temperature rock and the stress distributions of rock during injecting LN into the reservoirs. The schematic of model is shown in Figure 1. This model is carried out by the thermal-hydraulic-mechanical coupling method in transient state. The geometric size of the model is shown in Table 1.

The model mainly includes two regions, the fluid region and the solid region. In the fluid region, the high-pressure LN is injected into a pipe and then flows through the annulus between the drill hole and pipe. The inlet and outlet boundary conditions of fluid region are set as pressure inlet and pressure outlet boundary condition, respectively. During LN flowing through the hole of rock, the heat of the rock will be transferred to the cryogenic fluid quickly, which results in rock temperature decreasing. Consequently, the thermal stress is generated in reservoir rock, leading to the variation of original stress distribution and even rock damage.

Based on the present analysis, the following assumptions are proposed: (1) the influence of seepage flow in rock is ignored in the calculation. (2) No deformation of rock is assumed in flow field computation. (3) The rock is assumed as a homogeneous, isotropic, and linearly elastic material. (4) The calculating process does not involve rock failure.

2.2. Numerical Model

For the progress of injecting fluid into the wellbore, the standard k-ɛ model is adapted for the simulation of turbulent flow because of the larger Reynolds number. Besides, the intense heat transfer will also occur between the cryogenic fluid and warm rock. Thus, the equations of mass conservation, momentum, and energy should be solved in the flow field.Continuity equation:where ρl is fluid density, kg/m3; t is time, s; and v is velocity vector, m/s.Momentum equation:where p is pressure, Pa; μ is dynamic viscosity, N·s/m2; and g is gravitational acceleration, m2/s.Energy equation of fluid:where Tl is fluid temperature, °C; λl is thermal conductivity of fluid, W/m·K; and cp is specific heat at constant pressure, J/kg·K.Turbulent kinetic energy equation:where k is turbulent kinetic energy, m2/s2; is the velocity, m/s; μt is turbulent viscosity, N·s/m2; σk is the turbulent Prandtl number for k; Gk is generation of turbulence kinetic energy; ɛ is specific dissipation rate, J/kg·s; and YM is the contribution of the fluctuating dilatation.Specific dissipation rate of turbulent kinetic energy:where σɛ is the turbulent Prandtl number for ɛ, and C and C are turbulent constant.

Besides, μt is defined aswhere Cμ is turbulent constant. The empirical constants appearing in the above equations are given by the following values: C = 1.44, C = 1.92, Cμ = 0.09, σk = 1.0, σɛ = 1.3.

In the solid region, the heat conduction equation will be solved to calculate the temperature distribution. Furthermore, the stress state of rock owing to thermal stress is also considered in our model. The governing equations of the thermoelastic model are considered in the calculation.

The heat conduction of equation for rock:where λs is thermal conductivity of rock, W/m·K; c is specific heat of rock, J/kg·K; ρs is rock density, kg/m3; and Ts is rock temperature, °C.Physical equation:where τ is the strain of rock; G is shear modulus; σ is the total stress; is Poisson’s ratio; Θ is the sum of the normal stresses; I is the Kronecker delta; α is the coefficient thermal expansion, W/m·K; and ΔTs is rock temperature difference, °C.Equilibrium equation:where ρs is rock density.

The strain-displacement relation iswhere u is the displacement.

2.3. Boundary Conditions

To analyze the temperature and stress distributions of the rock during the process of injecting LN into the wellbore, the simulation parameters are set as in Table 2. Besides, it is assumed that the fluid region is initially filled with gaseous nitrogen (GN). During the simulation, the GN and rock are set to the same initial temperature. The LN is injected from the pipe inlet (i.e., inlet boundary condition), and its temperature is regarded as the injection temperature. In the solid region, the geological conditions of wellbore Newberry are referenced. It is assumed that the wellbore is located at the subsurface of 1.5 km depths, whose temperature is 150°C. The maximum principal stress is vertical with a gradient of 24.1 MPa/km, and the minimum principal stress is 14.9 MPa/km. Therefore, the front, left, and bottom surface of the model are set as displacement boundary conditions, in which the normal displacement is zero. Both the right and backside surfaces are set as stress boundary conditions, and their value is 22.35 MPa. The upper surface is set as stress boundary condition with the value of 36.15 MPa [35, 36].

In addition, the surface contacted the fluid region, and the solid region is set as the fluid-solid interface, which causes a large amount of heat exchange. The conjugate heat transfer method is used to simulate the heat transfer in the fluid-solid interface. Under conjugated boundary conditions, the temperature and the heat flux of the fluid and solid are equal, respectively. The corresponding mathematical expressions are as follows:where the Tl and Ts are the fluid temperature and the solid temperature at the solid-fluid interface, respectively; λl is the thermal conductivity of fluid; λs is the solid’s thermal conductivity; and n is the common normal direction of the solid-fluid interface.

2.4. Material Parameters

In the progress of the simulation, the detailed modeling parameters used are shown in Table 3. It is assumed that the parameters, including density, specific heat, thermal conductivity, and viscosity of LN, are unchanged in simulation progress. For injecting LN into the wellbore, the injection temperature is set as -173°C. In addition, the granite properties in calculation are referred to in the experimental data by Park [37]. Regarding the thermal-physical properties of granite, the fitting equations of specific heat, thermal conductivity, and thermal expansion with temperature are adopted in the calculation based on previous studies [37, 38].

The specific heat of granite:where Cs is rock specific heat, J/kg·K.

The thermal conductivity coefficient of granite:where λs is thermal conductivity coefficient of rock, W/m·K.

The thermal expansion coefficient of granite:where αs is thermal expansion coefficient of rock.

3. Results and Discussion

3.1. Mesh Independence and Time Step Selection

To ensure the accuracy of the simulation, mesh independence is validated with four different cell number cases (see Table 4), in which the primitive rock temperature is set as 150°C, the time step is set as 0.2 s, and other parameters remain unchanged. The point below the bottom surface of wellbore 1 mm is employed as the monitoring point. Variation of temperature against grid number is plotted in the logarithmic coordinate shown in Figure 2(a). Temperature just changes by 0.25% as the grid number increases from 2071302 to 2250028; thus, Case 3 is adopted in the following calculation.

In the meantime, time step is also an important influencing parameter in transient calculation. Different time steps, including 0.05 s, 0.1 s, 0.2 s, 0.5 s, and 1 s, are computed to verify that the time step is independent of simulation results. As shown in Figure 2(b), when time step is less than 0.1 s, the temperature of the monitoring point changes slightly; however, it increases sharply when time step is beyond 0.1 s. Considering computational precision and efficiency, 0.1 s is selected as the optimal time step for the computation in later simulation.

3.2. Flow Field Analysis

The fluid flow numerical solution is performed based on the CFD code of Fluent with coupled algorithm. It can simulate flow field and heat transfer by coupling mass conservation, momentum, energy, and heat conductivity equations. The computation results indicate that the parameters of pressure and velocity change slightly with time during the process of computation. Thus, the steady state flow is analyzed in this section. Figure 3 shows the contours of velocity and static pressure in flow field.

As shown in Figure 3(a), the velocity of LN increases near the outlet of tube and then decreases sharply at the bottom of the wellbore. Besides, at the center point of the bottom of well, the velocity of the LN drops to zero due to the boundary effect, and the hydrostatic pressure increases greatly at the same time as shown in Figure 3(b). Then, the LN flows diffusively around the stagnation point at the action of hydrostatic pressure. Moreover, there is a low velocity area near the pipe outlet due to the effect of fluid turbulence, which results in a lower pressure in the near area.

3.3. Rock Temperature Distribution

The temperature contours in the symmetry plane of model at different times are shown in Figure 4. The primitive rock temperature is set as 150°C in this case. At the action of pressure, the LN flows into the bottom of well through the pipe, which results in a lower temperature in fluid-solid interface at 0.1 s. The cooled-region of rock constantly increases with time going by. Because the cooled-region is circular in rock, the radius of cooled-region is chosen to compare the effect of heat transfer at different conditions. Specifically, it is expressed by the distance along path 1 and path 2 (see Figure 1). Path 1 is the line along the axial direction of rock center at the bottom of well, and path 2 is 100 mm above the bottom of well along the radial direction in rock. The temperature distribution at different times along path 1 and path 2 is shown in Figure 5. It is noted that the initial point temperature of path 1 and path 2 quickly drops to −156.3°C and −136.6°C, respectively, as soon as LN connects to rock for 0.1 s. However, both the cooling distances of path 1 and path 2 are about 24 mm at 20 s. During the process of injecting LN into wellbore, the fluid comes into contact with the rock of wellbore bottom first, which results in a low temperature region at the bottom rather than the side area at 0.1 s. However, with the high-speed flow of the LN, the temperature of solid-fluid interface will remain the same, and the distance of cooled-region perpendicular to the solid-fluid interface will also be similar.

For the inner region of rock, its cooled-region is constantly increasing, but the growth trend in cooled-region becomes slower and slower. Based on (10) and (11), we could conclude that, with rock temperature rising, the thermal conductivity drops, while the specific heat increases. In this case, more heat is required to be carried away to cool the rock. As a result, the growth of cooled-region in inner rock will be slowed down.

3.4. Maximum Principal Stress Analysis

The thermal stress will be generated due to the rock temperature decreasing. Therefore, the stress state in rock would be changed with the effects of thermal stress and fluid pressure. A vertical plane 1 which is 100 mm above the downhole is chosen to show the maximum principal stress contour of rock at different times (see Figure 6). Here, the positive and negative values represent the tensile and compressive stresses, respectively. The thermal stress is generated by the temperature gradient inside the rock with cryogenic LN flowing along the surface of the rock. In addition, it can be found that tensile stress zone increases with time going by due to the inside rock cooling down. This shows that injecting LN can change the original stress state and the fracturing environment.

Figure 7 illustrates the maximum principal stress distributions at different times along path 2. There are small changes for maximum principal stress values at different times except 0.1 s. This is because the surface temperature of rock is in unsteady state at 0.1 s. This also indicates that the maximum principal stress is gradually increasing during the period of LN contact with rock surface. The value of tensile stress generated at 100 s is much greater than the tensile strength of granite, which shows that it is feasible to break rock by injecting LN. Besides, the tensile stress zone gradually increases with time going by. This is favorable for the rock to form complex fracture networks [18].

3.5. Parameter Analysis

Figure 8 shows the maximum principal stress along path 2 with different primitive rock temperature at 20 s. The results indicate that, with primitive rock temperature rising, the value of maximum principal stress also increases along path 2. The maximum principal value of stress is about 154.55 MPa for the rock with the primitive temperature of 300°C, which is much greater than the value with the primitive temperature of 150°C. The reason is that larger temperature difference is generated on rock surface with the primitive rock temperature rising, and it also leads to much larger deformation correspondingly. In addition, when the primitive temperature increases from 150°C to 300°C, the zone of tensile stress zone increases due to the larger maximum principal stress being generated on the rock surface, which is beneficial to the cracking and breaking of rock.

The difference between the inlet and ambient pressure is defined as the pressure drop. Figure 9 shows the maximum principal stress along path 2 with different pressure drop at 20 s. The results indicate that the distributions of maximum principal stress are hardly affected by the inlet pressure with the same ambient pressure. The reason is that the thermal conductivity and the heat transfer rate of rock are limited; the rock cooling performance cannot be enhanced efficiently by improving the inlet pressure. When the ambient pressure of reservoir has not been changed, the thermal stress induced by rock cooling is the main factor on the rock stress distribution.

The ambient pressure is also a significant factor on rock fracturing. The maximum principal stress along path 2 with different ambient pressure is shown in Figure 10. Similarly, there is a small change in the distribution of stress along path 2. Because the inlet pressure also hardly affects the rock cooling performance during the LN flowing along wellbore, the thermal stress of rock almost remains unchanged. The value of maximum principal stress is slightly increased due to the effect of fluid pressure.

To study the influence of different stress boundary on stress distribution, we calculate the maximum principal stress with different stress boundary as shown in Figure 11. It indicates that tensile stress will appear on the coupled surface when stresses on both sides of the plane are equal, and the maximum principal stresses along the circumference are similar. When the stresses on both sides are unequal, the maximum principal stress on the coupled surface of wellbore will change. With the increasing of σH, the maximum principal stress area in the axial direction will decrease gradually, and the maximum principal stress near wellbore along the radial direction also decreases at the same time. The greater the difference is, the more obvious the trend will be. Figure 12 illustrates the value of maximum principal stress in rock along path 3 under different stress boundary. It could be clearly seen that maximum principal stress near coupled interface decreases with σH improving. Thus, it could be reasonably predicted that the rock will be fractured toward the direction of larger in situ stress when LN flows into wellbore.

To analyze the effect of thermal stress induced by LN cooling, the water cooling is also considered under the same conditions. The physical properties of water are set as follows: density is 998.2 kg/m3, specific heat is 4182 J/(kg·K), thermal conductivity is 0.6 W/(m·K), and viscosity is 0.001003 kg/(m·s). When the stress caused by injecting water into the reservoirs is calculated, the inlet temperature is set as 25°C. Figure 13 shows the maximum principal stress along path 2 with LN and water flowing into wellbore at 100 s, respectively. In this case, the maximum principal stress in the side surface of rock is about 21.02 MPa at 100 s with the water flowing into wellbore, which is about 39.90% caused by LN under the same conditions. This indicated that, due to the extremely cryogenic characteristic, LN will generate greater tensile stress on the rock surface during the progress of flowing into wellbore compared with water. It is beneficial to reduce the injection pressure.

In order to study the effect of inlet pressure on stress distribution around wellbore as water flows into the bottom of the wellbore, the maximum principal stress of wellbore surface at 100 s under the same pressure drop is calculated. Figure 14 shows the maximum principal stress in rock with different inlet pressure at 100 s. It is noted that there is an obvious linear relationship between the water inlet pressure and the maximum principal stress. The fitting expression isWhere P is inlet pressure and σ is maximum principal stress.

In addition, the value of maximum principal stress is about 50.65 MPa with the inlet pressure of 55 MPa, which is similar to the value caused by LN flowing at 25 MPa under the same conditions. This shows that LN injection into the wellbore for fracturing operations can effectively reduce the inlet pressure due to its cryogenic properties compared with fracturing by water.

The variation of stress values due to water injection under different in situ stress conditions was analyzed. Figure 15 shows the values of maximum principal stress at the surface of wellbore with the relationship of σH. It indicated that as σH increases, the value of maximum principal stress decreases at the surface of wellbore in the direction of minimum horizontal in situ stress, which is the same trend as the effect of LN injection on the values of maximum principal stress. However, when σH is 32.35 MPa, the maximum principal stress in the direction of maximum horizontal in situ stress is about 31.04 MPa generated by the water injection, which is lower than the stress value of 35.04 MPa at the direction of minimum horizontal in situ stress by the injected LN. It can be concluded that the tensile damage could be generated in the direction of minimum horizontal in situ stress by the LN injection during the fracture process by injecting water and LN at the same inlet pressure. This is mainly due to the fact that the LN injection process generates a huge thermal stress in the surface of wellbore. The value of thermal stress is not affected by the in situ stress distribution, so the thermal cracks are formed under the thermal stress, which is beneficial to the cracks expansion. Based on this, a more complex fracture network is formed in the reservoir. It has been demonstrated that the main fracture direction of rock samples after LN fracturing does not extend exactly in the direction of the maximum principal stress [39].

4. Conclusion

In this paper, the heat transfer and stress state during injecting LN into HDR reservoirs are analyzed based on a 3D thermal-hydraulic-mechanical coupling numerical model. The parameters’ sensitivity is analyzed in detail, and the injection of water under the same conditions is also considered. Finally, a set of experiments is conducted to validate the effect of thermal stress on the rock. The main conclusions of this paper are as follows:(1)The LN injected into the wellbore causes the rocks around the wellbore to cool down. As time passes by, the increasing amplitude of cooled-region is weakened. During this progress, the great tensile stress that exceeds the tensile strength of granite is generated around the wellbore due to the action of thermal stress and fluid pressure.(2)The primitive temperature of the reservoir has a significant impact on the stress distributions around the wellbore during LN injecting. With the growth of the primitive temperature, the thermal stress value around wellbore becomes larger. The inlet pressure and the ambient pressure have little effect on the thermal stress generated by LN cooling.(3)The in situ stresses affect the stress distribution during LN fracturing. Under the unequal in situ stress in different directions, the value of maximum principal stress is larger in the direction with larger in situ stress; i.e., the rock has a tendency to break in the direction of the larger in situ stress value.(4)Compared with the injection of LN, the value of thermal stress around the wellbore caused by water injection is reduced under the same conditions. The thermal stress value around the wellbore during LN injecting with injection pressure of 25 MPa is similar to that during water injecting with injection pressure of 55 MPa. This means that LN fracturing can reduce the injection pressure effectively and achieve a better fracturing effect.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

Chengzheng Cai contributed to conceptualization. Keda Ren participated in data curation and reviewed and edited the manuscript. Both authors were responsible for methodology, software, supervision, and original draft preparation.


This research was funded by the National Key R&D Program of China (2020YFA0711800), the Assistance Program for Future Outstanding Talents of China University of Mining and Technology (2020WLJCRCZL010), the Postgraduate Research and Practice Innovation Program of Jiangsu Province (KYCX20_2042), and the National Natural Science Foundation of China (51604263).