Study on Temperature Distribution Law of Tunnel Portal Section in Cold Region Considering Fluid–Structure Interaction

Abstract

The design of tunnels in cold regions contributes greatly to the feasibility and sustainability of highways. Based on the heat transfer mechanism of the tunnel surrounding rock–lining–air, this paper uses FEPG software to carry out secondary excavation and development, then the air heat convection calculation model is established by using a three-dimensional extension of the characteristic-based operator-splitting (CBOS) finite-element method and the explicit characteristic–Galerkin method. By coupling with the heat conduction model of the tunnel lining and surrounding rock, the heat conduction-thermal convection fluid–structure interaction finite-element calculation model of tunnels in cold regions is established. Relying on the Qinghai Hekashan tunnel project, the temperature field of the tunnel portal section is calculated and studied by employing the fluid–structure interaction finite-element model and then compared with the field monitoring results. It is found that the calculated values are basically consistent with the measured values over time, which proves the reliability of the model. The calculation results are threefold: (1) The temperature of the air, lining, and surrounding rock in the tunnel changes sinusoidally with the ambient temperature. (2) The temperature of each layer gradually lags behind, and the temperature variation amplitude of the extreme value of the layer temperature gradually decreases with the increase in the radial distance of the lining. (3) In the vicinity of the tunnel entrance, the lining temperature of each layer remains unchanged, and the temperature gradually decreases or increases with the increase in the depth. The model can be used to study and analyze the temperature field distribution law of the lining and surrounding rock under different boundary conditions, and then provide a calculation model with both research and practical value for the study of the temperature distribution law of tunnels in cold regions in the future.

1. Introduction

With the rapid development of China’s economy and society, tunnel excavation, oil and gas development, coal mining, and other underground space development projects are also developing rapidly. The design of tunnels in cold regions contributes greatly to the feasibility and sustainability of highways, which is an inevitable requirement for the development of major national infrastructure. Among them, the rock mechanics of tunneling projects constructed in cold regions are outstanding, as tunnels in cold regions often have serious frost damage [1,2,3,4,5]. After the tunnel is connected, hot air will flow inside, gradually melting permafrost if the surrounding rock is permanently frozen. This will significantly reduce the strength of the rock and compromise the safety of the tunnel’s operation. If the surrounding rock of the tunnel is non-frozen soil, the low-temperature air inflow after the tunnel penetration causes seasonal freeze–thaw cycles in the non-frozen soil. The safety of the lining structure will be significantly impacted by this phenomena, which will further create frost heaving force and load it on the lining [6,7].

The frost damage of tunnels in cold regions will cause a series of lining diseases, such as cracking of the lining vault, serious weathering and spalling of the lining surface, extrusion of the side wall, leakage of the lining, and hanging ice [8]. Following the lining’s exposure to freezing damage, it will significantly decrease its structural safety and stability, threaten the safety of oncoming vehicles, and destroy auxiliary facilities like lighting and ventilation in the tunnel. This will make it more difficult for the tunnel to operate normally and significantly raise the cost of maintenance. There are many factors affecting the freezing damage of tunnels in cold regions, including the geographical location, topography, geological structure and hydrological conditions, thermal parameters of the surrounding rock and lining, local meteorological conditions, and even the traffic volume and type of tunnel operation. The operation of ancillary facilities includes lighting facilities, ventilation, and so on [9].

How to accurately grasp the temperature distribution law of the tunnel surrounding rock and lining is the premise and key to solving the problem of freezing damage in cold regions. At present, domestic and foreign scholars’ research on the temperature field of tunnels in cold regions is mainly divided into two aspects: theoretical analysis and numerical simulation. In terms of theoretical analysis, Li et al. studied the theoretical solution method of the temperature field of tunnels in cold regions, and deduced the heat transfer model of the airflow according to the basic methods of heat transfer and fluid mechanics [10]. Sun et al. [11,12] focused on the study of the temperature field of tunnels in cold regions under the coupling effect, obtained the temperature field distribution law of the longitudinal and transverse sections of the tunnel, and analyzed the influence of the initial ground temperature of the surrounding rock. At the same time, considering that the temperature field in the tunnel is affected by a certain degree of meteorological factors, combined with previous research on the initial ground temperature of the surrounding rock, the temperature field of the tunnel in the meteorological–geological system was quantitatively analyzed by taking the temperature and initial ground temperature in the tunnel and the action time as variables. Taking the Qilian Mountain Tunnel as the research object, Gao et al. used the control variable method to explore the distribution law of the air and surrounding rock temperature field in a tunnel under different time and space conditions [13] and obtained the relationship between the air temperature outside the tunnel, the surrounding rock temperature, and the tunnel temperature field. Zhou et al. deduced the control equation of a tunnel temperature field. Combined with a comparison with the measured temperature field, the numerical solution of the temperature field of the tunnel lining surrounding rock was obtained [14]. Zhang et al. collected field data and obtained the temporal and spatial distribution law of the tunnel temperature field by a regression fitting formula [15]. In terms of numerical simulation, Zhang et al. established a three-dimensional temperature field model of a tunnel on the basis of the above-mentioned field measurement and explored the temperature field distribution law under multi-field coupling [16]. Based on the mixture theory, considering the factors of water ice phase change and water migration in a tunnel, Zhang et al. established a hydrothermal coupling model to analyze the distribution law of the temperature field in a tunnel under the dynamic action of water state change [17]. Yu et al. focuses on water–heat coupling and considered random cracks as action factors to establish a coupling model and obtained the influence of the heat transfer coefficient and surrounding rock conditions on the tunnel temperature field [18]. Ma et al. combined numerical simulation data with information collected on site and concluded that the temperature distribution of each section is closely related to the temperature in the tunnel as the tunnel depth increases [19]. Ding et al. analyzed and calculated a finite-element model by the control variable method [20]. In addition, the thermal stress was found to cause damage inside the rock, resulting in variations in physical and mechanical properties [21,22,23]. It was concluded that the convective heat transfer of air is the main factor affecting the temperature field of a tunnel under the same initial ground temperature. From the above research, it can be seen that the current research carried out by practitioners in the field is mainly based on theoretical analysis, and the theoretical results of the research have been quite fruitful. According to the direction of numerical simulation, more research data are still needed to provide support for the development of theory. The main numerical simulation methods used for rock mechanics problems are the finite-element method and the discrete-element method [24,25,26,27,28,29,30]. At the same time, it is not difficult to see from the above research that the current mainstream research on the temperature field of tunnels in cold regions involves the study of the temperature distribution under multi-field coupling, which mainly explores the influence of the surrounding rock conditions on the temperature field distribution [11]. Furthermore, the influence of water–heat coupling in a tunnel was analyzed [17,18]. Ma et al. reflected the influence of air flow on the temperature field in a tunnel, and the effect of air flow was mainly reflected in the convective heat transfer between air and surrounding rock-lining system [19,20,31]. As mentioned above, under the premise of multi-field coupling, the influence of air convection heat transfer on tunnel temperature cannot be ignored. Some scholars have begun to pay attention to the influence of convection heat transfer. In addition, there is a certain heat conduction between the lining surface and the air contact part. The heat transfer process of the lining–air involves the process of heat conduction and heat convection. The systematic study of the tunnel temperature field must consider the effect of heat conduction while studying the air convection heat transfer.

In view of this, on the basis of the above research, this paper takes the air–lining–surrounding rock system in the tunnel as the research object and uses Finite-Element Program Generator (FEPG V8.0) software for secondary development. Considering the dynamic influence of air velocity and temperature on the lining surface in the tunnel, the N-S equation in the air thermal convection control equation is extended from one-dimensional to three-dimensional, and the algorithm of the tunnel air thermal convection equation is obtained. Then, the heat conduction effect is coupled with the thermal convection effect, and a comprehensive and accurate fluid–structure interaction model is proposed to provide a reference for subsequent research.

2. The Basic Theory of Heat Conduction and Heat Convection Calculation in Cold Tunnel

2.1. Control Equation of Heat Conduction in Tunnel in Cold Region

Due to the close contact between the primary lining of the tunnel and the surrounding rock in the project, the nature is close, and the thickness is thin. This paper regards the primary lining and the surrounding rock as a whole, mainly taking the secondary lining (hereinafter referred to as the lining) as the research object. The following assumptions are made in the finite-element model of heat transfer in a cold tunnel: (1) The tunnel lining and surrounding rock are isotropic continuous media. (2) The thermal conductivity, specific heat capacity, and lining and surrounding rock density are fixed constants, which are not affected by other factors.

According to the assumption, since there is no internal heat source in the operating tunnel, the heat transfer differential equation is as follows:

$${\lambda }_c\left[\frac{ {\partial }^{2}T}{\partial x^2}+\frac{ {\partial }^{2}T}{\partial y^2}+\frac{ {\partial }^{2}T}{\partial z^2}\right]={\rho }_c{c}_c\frac{ \partial T}{\partial t^{}}$$

(1)

${\lambda }_c$—Thermal conductivity of lining and surrounding rock, W/(m·°C); $c_c$—Specific heat of lining and surrounding rock, J/(kg·°C); ${\rho }_c$—Lining and surrounding rock density, kg/m³; T—Temperature of any point inside the lining and surrounding rock, °C; t—Time, seconds.

Initial condition. At the initial instant, the tunnel temperature field is a known function (x,y,z) of coordinate T₀(x,y,z), that is, when,

T(x,y,z,0) = T₀(x,y,z)

(2)

In most cases, the initial instantaneous temperature distribution can be considered as a constant, that is, when t = 0,

T(x,y,z,0) = T₀

(3)

Boundary condition. For the first boundary condition, the surface temperature of the lining and surrounding rock is a known function of time:

$$T(t)=f(t)$$

(4)

For the second boundary condition, the heat flux density q on the surface of the lining and surrounding rock is a known function of time:

$$q=-{\lambda }_c\frac{\partial T}{\partial n}=f(t)$$

(5)

For the third boundary condition, the air temperature and convective heat transfer coefficient in contact with the lining are known, and the heat flow on the lining surface is proportional to the difference between the lining surface temperature and the air temperature:

$$-{\lambda }_c\frac{ \partial T}{\partial n}=\beta \left(T-T_a\right)$$

(6)

$T_a$—Air temperature, °C; $\beta$—Coefficient of convective heat transfer, W/(m²·°C).

The finite-element solution needs to transform Formula (1) from the partial differential form to the integral form, and then the discrete format is obtained according to the principle of the Galerkin method.

2.2. Control Equation of Air Thermal Convection in Tunnel in Cold Region

The following assumptions are made in the finite-element model of air thermal convection in a cold tunnel: (1) Air is an incompressible Newtonian fluid. (2) Air thermal conductivity, specific heat capacity, density, dynamic viscosity, and other physical parameters are constants. (3) Air pressure does not change with temperature. (4) Neglect the influence of gravity. (5) Due to the low air velocity, the air viscous dissipation (friction loss) is ignored [31].

The basic governing equations of fluid thermal convection can be derived by using the law of conservation of mass, the law of conservation of momentum, and the law of conservation of energy. Combined with the above assumptions, the final control equation can be obtained as follows:

(1)

Continuity equation

$$\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}=0$$

(7)

$u_x,u_y,u_z$—Velocity of fluid in X, Y, and Z directions, respectively, m/s.

In tensor form:

$$\frac{\partial u_i}{\partial x_i}=0$$

(8)

(2)

The Navier–Stokes (N-S) equations of motion

The motion equation tensor form of viscous fluid:

$${\rho }_a\frac{\partial u_i}{\partial t}+{\rho }_a{u}_j\frac{\partial u_i}{\partial x_j}=-\frac{\partial p}{\partial x_i}+{\mu }_a\frac{\partial }{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$

(9)

${\mu }_a$—Aerodynamic viscosity, $\mathrm{Pa}\cdot \mathrm{s}$; i, j—Three-dimensional coordinates, take 1,2,3 respectively; ${\rho }_a$—Air density, $\mathrm{kg}/{\mathrm{m}}^3$.

(3)

Energy-conservation equation

$${\rho }_a{\mathrm{c}}_a\left(\frac{\partial T_a}{\partial t}+u_x\frac{\partial T_a}{\partial x}+u_y\frac{\partial T_a}{\partial y}+u_z\frac{\partial T_a}{\partial z}\right)={\lambda }_a\left[\frac{{\partial }^2{T}_a}{\partial x^2}+\frac{{\partial }^2{T}_a}{\partial y^2}+\frac{{\partial }^2{T}_a}{\partial z^2}\right]$$

(10)

In tensor form:

$$\frac{\partial T_a}{\partial t}+u_i\frac{\partial T_a}{\partial x_i}={\alpha }_a\frac{{\partial }^2{T}_a}{\partial x_i{}^2}$$

(11)

$T$—Flow temperature, $°C$; $\lambda$—Fluid thermal conductivity, $\mathrm{W}/\left(\mathrm{m}\cdot °C\right)$; $c$—Specific heat of fluid, $\mathrm{J}/\left(\mathrm{kg}\cdot °C\right)$; ${\alpha }_a$—Air thermal conductivity, also known as thermal diffusivity or thermal diffusivity, ${\mathrm{m}}^2/\mathrm{s}$, ${\alpha }_a=\frac{{\lambda }_a}{{\rho }_a{c}_a}$.

Initial condition: In the initial instantaneous, the fluid temperature distribution can generally be set to a constant.

For unsteady flow, the initial state of the flow field should be given, that is, the velocity field, pressure field, and temperature field when t = 0. For steady flow, the flow state is the same at any time, the local derivative term in the basic equation is equal to 0, and the initial conditions are not given.

Boundary condition:

①

Boundary conditions on the fluid–structure interface

Non-slip boundary conditions:

$$\overrightarrow{u}=\overrightarrow{u_b}$$

(12)

$\overrightarrow{u}$—Fluid velocity, m/s; $\overrightarrow{u_b}$—Wall velocity, m/s.

In the steady flow around a fixed object, the wall velocity is zero, so the non-slip condition is

$$\overrightarrow{u}=0$$

(13)

When the viscosity of the fluid is ignored, only the normal velocity on the solid wall is continuous, and there can be relative slip between the fluid and the solid wall. This condition is called the slip boundary condition.

②

Heat balance condition of solid wall temperature

The fluid particle temperature is the same as the wall temperature, that is,

$$T_{\mathrm{w}}=T_b$$

(14)

$T_{\mathrm{w}}$—Fluid particle temperature, $°C$; $T_b$—Wall temperature, $°C$.

③

Adiabatic solid wall conditions

The heat flux density of the adiabatic solid wall is zero, which is the second boundary condition of heat transfer, that is,

$${\left(\frac{\partial T}{\partial n}\right)}_{\mathrm{w}}=0$$

(15)

④

Kinematic, kinetic, and thermodynamic conditions at the fluid interface

If the surface tension is not taken into account, the velocity, stress, and temperature of the fluid at any point on the interface (represented by subscripts “+” and “−“, respectively) should be equal on the non-invasive interface between the two fluids, that is,

$$\overrightarrow{u_{+}}=\overrightarrow{u_{-}},\hspace{1em}{\sigma }_{+}={\sigma }_{-},\hspace{1em}{T}_{+}=T_{-}$$

(16)

3. Three-Dimensional Finite-Element Solution of Air Thermal Convection Model of Tunnel in Cold Region

3.1. Characteristic Line Operator-Splitting Finite-Element Method for N-S Equation

The characteristic-based operator-splitting finite-element method (CBOS method for short) is proposed by Wang et al. [32] to improve the Taylor–Galerkin method by Taylor formula expansion and operator splitting (CBS). In this method, the operator-splitting method is used at each time layer to separately solve the convection term and the diffusion term of the N-S equation. The backward difference scheme is used for the diffusion term’s time discretization, while the standard Galerkin method is used for the diffusion term’s implicitly solved space discretization. The convection term is discretized by the characteristic line–Galerkin method in time, while the space is discretized by the standard Galerkin method, and the explicit solution is obtained [33,34].

According to the above research, the CBOS method is extended to solve the three-dimensional unsteady incompressible fluid.

3.1.1. Operator Splitting of N-S Equation

Based on the idea of CBOS method, the N-S equation of the convective heat transfer model in a cold-region tunnel is split into the diffusion term

$$\rho \frac{u_i^{n+{\theta }_1}-u_i^n}{\Delta t}-\mu \frac{\partial }{\partial x_j}\left(\frac{\partial u_i^n}{\partial x_j}+\frac{\partial u_j^n}{\partial x_i}\right)=\rho f_i$$

(17)

convection item

$$\frac{\partial u_i^{n+{\theta }_2}}{\partial t}+u_j^{n+{\theta }_2}\frac{\partial u_i^{n+{\theta }_2}}{\partial x_j}=0$$

(18)

and pressure term

$$\rho \frac{u_i^{n+1}-u_i^{n+{\theta }_2}}{\Delta t}=-\frac{\partial p^{n+1}}{\partial x_i}$$

(19)

$$\frac{\partial u_i^{n+1}}{\partial x_i}=0$$

(20)

Obtain the pressure Poisson equation:

$$\frac{{\partial }^2{p}^{n+1}}{\partial x_i^2}=\frac{\rho }{\Delta t}\frac{\partial u_i^{n+{\theta }_2}}{\partial x_i}$$

(21)

$u_i^n$—Value of velocity field at time n, $\mathrm{m}/\mathrm{s}$; $p^n$—Value of pressure field at time n, $\mathrm{Pa}$; $u_i^{n+{\theta }_1}$—Solution of the velocity field at the time of the diffusion term n + 1 and the initial value of the convection term n + 1, $\mathrm{m}/\mathrm{s}$; $u_i^{n+{\theta }_2}$—Solution of the velocity field at the moment of convection term n + 1 and the initial value of the pressure correction term n + 1, $\mathrm{m}/\mathrm{s}$; $u_i^{n+1}$—Value of velocity field at time n + 1, $\mathrm{m}/\mathrm{s}$; $p^{n+1}$—Value of pressure field at time n + 1, $\mathrm{Pa}$; $\Delta t$—time step.

3.1.2. The Display Time Discretization and Display Format of the Convection Term along the Characteristic Line

(1)

The display time discretization of the convection term along the characteristic line

Through differential simplification and Taylor formula expansion, a fully explicit time discretization scheme is obtained:

$$\begin{array}{l}{u}_i^{n+{\theta }_2}-u_i^{n+{\theta }_1}=-\Delta t{u}_j^{n+{\theta }_1}\frac{\partial u_i^{n+{\theta }_1}}{\partial x_j}-\frac{\Delta t^2}{2}\frac{\partial u_j^{n+{\theta }_1}}{\partial t}\frac{\partial u_i^{n+{\theta }_1}}{\partial x_j}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{\Delta t^2}{2}{u}_k^{n+{\theta }_1}\frac{\partial }{\partial x_k}\left(u_j^{n+{\theta }_1}\frac{\partial u_i^{n+{\theta }_1}}{\partial x_j}\right)\end{array}$$

(22)

(2)

Multi-step display format of convection term

The multi-step format is used to deal with the convection term [35], and the solution formula of the convection term in the lth sub-time step is obtained:

$$\begin{array}{l}{u}_i^{n+\frac{l}{h}}-u_i^{n+\frac{l-1}{h}}=-\delta t{u}_j^{n+\frac{l-1}{h}}\frac{\partial u_i^{n+\frac{l-1}{h}}}{\partial x_j}+\frac{\delta t^2}{2}{u}_k^{n+\frac{l-1}{h}}\frac{\partial u_j^{n+\frac{l-1}{h}}}{\partial x_k}\frac{\partial u_i^{n+\frac{l-1}{h}}}{\partial x_j}\\ \hspace{1em}+\frac{\delta t^2}{2}{u}_k^{n+\frac{l-1}{h}}\frac{\partial }{\partial x_k}\left(u_j^{n+\frac{l-1}{h}}\frac{\partial u_i^{n+\frac{l-1}{h}}}{\partial x_j}\right)\end{array}$$

(23)

3.1.3. Three-Dimensional Finite-Element Solution

(1)

Three-dimensional finite-element discretization of diffusion term

The diffusion term is discretized by the CBOS method, and the typical element is a hexahedral element with eight nodes by the finite-element solution method.

$$u_i^{(e)}=u_{i\alpha }^{(e)}{N}_{\alpha }^{(e)}; \delta u_{j\alpha }^{(e)}=N_{\alpha }^{(e)}$$

(24)

$\alpha =1,2,\dots ,m$, $m=8$ is the number of velocity nodes in the interpolation unit.

Obtain the unit finite-element equation:

$$A_{\alpha \beta }^{(e)}{u}_{i\beta }^{n+{\theta }_1(e)}=O_{\alpha \beta }^{(e)}{u}_{i\beta }^{n(e)}+F_{i\alpha }^{(e)}+B_{i\alpha }^{(e)}$$

(25)

$\beta =1,2,\dots ,m$;

$$A_{\alpha \beta }^{(e)}=\rho {\displaystyle {\iiint }_{V^{(e)}}{N}_{\beta }^{(e)}{N}_{\alpha }^{(e)}d{V}^{(e)}}+\mu \Delta t{\displaystyle {\iiint }_{V^{(e)}}(\frac{\partial N_{\beta }^{(e)}}{\partial x_j}+\frac{\partial N_{\beta }^{(e)}}{\partial x_i}{\delta }_{ij})\frac{\partial N_{\alpha }^{(e)}}{\partial x_j}}d{V}^{(e)}$$

(26)

$$O_{\alpha \beta }^{(e)}=\rho {\displaystyle {\iiint }_{V^{(e)}}{N}_{\beta }^{(e)}{N}_{\alpha }^{(e)}d{V}^{(e)}}$$

(27)

$$F_{i\alpha }^{(e)}=\rho \Delta t{\displaystyle {\iiint }_{V^{(e)}}{f}_i^{(e)}{N}_{\alpha }^{(e)}d{V}^{(e)}}$$

(28)

$$B_{i\alpha }^{(e)}=\mu \Delta t{\displaystyle {\iiint }_{S^{(e)}}\frac{\partial u_i^{n+{\theta }_1(e)}}{\partial n}{N}_{\alpha }^{(e)}d{S}^{(e)}}$$

(29)

In Equation (29), ${\delta }_{ij}$ is the permutation operator, ${\delta }_{ij}{u}_{i\beta }=u_{j\beta }$.

By synthesizing the whole computational domain, the total finite-element equation is obtained.

$$A_{sq}{u}_{{}_{iq}}^{n+{\theta }_1}=O_{sq}{u}_{iq}^n+F_{is}+B_{is}$$

(30)

s, q—Overall speed node number.

The solution (Equation (30)) can obtain the velocity $u_i^{n+{\theta }_1}$.

(2)

Three-dimensional finite-element discretization of convection term

The convection term is discretized by using the CBOS method, and the finite-element equation is obtained by employing finite-element solution method:

$$O_{\alpha \beta }^{(e)}{u}_{i\beta }^{g(l+1)(e)}=C_{i\alpha }^{(e)}$$

(31)

$$\begin{array}{l}{C}_{i\alpha }^{(e)}=u_{i\beta }^{g(l)(e)}{\displaystyle {\iiint }_{V^{(e)}}{N}_{\alpha }^{(e)}}{N}_{\beta }^{(e)}dV-u_{j\gamma }^{g(l)(e)}{u}_{i\beta }^{g(l)(e)}\delta t{\displaystyle {\iiint }_{V^{(e)}}{N}_{\gamma }^{(e)}}\frac{\partial N_{\beta }^{(e)}}{\partial x_j}{N}_{\alpha }^{(e)}dV\\ +u_{j\eta }^{g(l)(e)}{u}_{i\gamma }^{g(l)(e)}{u}_{k\beta }^{g(l)(e)}\frac{\delta t^2}{2}{\displaystyle {\iiint }_{V^{(e)}}{N}_{\beta }^{(e)}}\frac{\partial N_{\eta }^{(e)}}{\partial x_j}\frac{\partial N_{\gamma }^{(e)}}{\partial x_k}{N}_{\alpha }^{(e)}dV\\ -u_{j\eta }^{g(l)(e)}{u}_{i\gamma }^{g(l)(e)}{u}_{k\beta }^{g(l)(e)}\frac{\delta t^2}{2}{\displaystyle {\iiint }_{V^{(e)}}{N}_{\eta }^{(e)}}\frac{\partial N_{\gamma }^{(e)}}{\partial x_j}{N}_{\beta }^{(e)}\frac{\partial N_{\alpha }^{(e)}}{\partial x_k}dV\end{array}$$

(32)

Synthesizing the whole computational domain, the total finite-element equation can be obtained:

$$O_{sq}{u}_{iq}^{g(l+1)}=C_{is}$$

(33)

The solution (Equation (33)) can obtain the velocity $u_i^{g(l+1)}$.

(3)

Three-dimensional finite-element discretization of pressure Poisson equation

The standard Galerkin method is used to discretize the pressure Poisson equation (Equation (21)), and the typical element is taken as an eight-node hexahedral element. Therefore,

$$p_i^{(e)}=p_{\alpha }^{(e)}{N}_{\alpha }^{(e)}$$

(34)

$$\delta p_i^{(e)}=N_{\alpha }^{(e)}$$

(35)

Obtain the element finite-element equation:

$$H_{\alpha \beta }^{(e)}{p}_{\beta }^{n+1(e)}=-Q_{\alpha \beta }^{(e)}{u}_{i\beta }^{n+{\theta }_2(e)}+D_{\beta }^{(e)}$$

(36)

$$H_{\alpha \beta }^{(e)}={\displaystyle {\iiint }_{V^{(e)}}\frac{N_{\beta }^{(e)}}{\partial x_i}\frac{N_{\alpha }^{(e)}}{\partial x_i}d{V}^{(e)}}$$

(37)

$$Q_{\alpha \beta }^{(e)}=\frac{\rho }{\Delta t}{\displaystyle {\iiint }_{V^{(e)}}\frac{\partial N_{\beta }^{(e)}}{\partial x_i}{N}_{\alpha }^{(e)}d{V}^{(e)}}$$

(38)

$$D_{\beta }^{(e)}={\displaystyle {\iiint }_{S^{(e)}}\frac{\partial p^{n+1}}{\partial n}{N}_{\beta }^{(e)}d{S}^{(e)}}$$

(39)

The total finite-element equation can be obtained by synthesizing the whole computational domain:

$$H_{sq}{p}_q^{n+1}=-Q_{sq}{u}_{iq}^{n+{\theta }_2}+D_q$$

(40)

The solution (Equation (40)) can obtain the pressure $p^{n+1}$.

(4)

Three-dimensional finite-element discretization of velocity correction term

Organize the Formula (19):

$${\displaystyle {\iiint }_{V^{(e)}}{u}_i^{n+1}}\delta u_{i}d{V}^{(e)}={\displaystyle {\iiint }_{V^{(e)}}{u}_i^{n+{\theta }_2}}\delta u_{i}d{V}^{(e)}-\frac{\rho }{\Delta t}{\displaystyle {\iiint }_{V^{(e)}}\frac{\partial p^{n+1}}{\partial x_i}\delta u_i}d{V}^{(e)}$$

(41)

The total finite-element equation can be obtained by synthesizing the whole computational domain:

$$O_{\alpha \beta }^{(e)}{u}_{i\beta }^{n+1(e)}=O_{\alpha \beta }^{(e)}{u}_{i\beta }^{n+{\theta }_2(e)}-Q_{\alpha \beta }^{(e)}{p}^{n+1(e)}$$

(42)

The solution (Equation (42)) can obtain the velocity $u_i^{n+1}$.

3.2. Explicit Characteristic Line–Galerkin Method for Energy-Conservation Equation

3.2.1. Explicit Characteristic Line–Galerkin Method

The fluid issue is solved by converting it into a set of ordinary differential equations on the characteristic line using the explicit characteristic line–Galerkin approach [36]. This method has such advantages as good stability, relatively simple programming, and high computational efficiency. The main idea is the local Taylor expansion.

The explicit characteristic line–Galerkin multidimensional tensor form of the energy equation:

$$T^{n+1}-T^n=-\Delta t{u}_i^n\frac{\partial T^n}{\partial x_i}+\Delta t\alpha \frac{{\partial }^2{T}^n}{\partial x_i^2}+\frac{\Delta t^2}{2}{u}_j^n\frac{\partial }{\partial x_j}\left(u_i^n\frac{\partial T^n}{\partial x_i}\right)$$

(43)

3.2.2. Three-Dimensional Finite-Element Solution

The Galerkin method is used to establish the integral form of Equation (43), and the finite-element equation of the element can be obtained according to the finite-element solution method:

$$A_{\alpha \beta }^{(e)}{T}_{\alpha }^{n+1(e)}=B_{\alpha \beta }^{(e)}$$

(44)

$$A_{\alpha \beta }^{(e)}={\displaystyle {\iiint }_{V^{(e)}}{N}_{\alpha }^{(e)}{N}_{\beta }^{(e)}d{V}^{(e)}}$$

(45)

$$\begin{array}{l}{B}_{\alpha \beta }^{(e)}={\displaystyle {\iiint }_{V^{(e)}}{T}_{\alpha }^{n(e)}{N}_{\alpha }^{(e)}{N}_{\beta }^{(e)}d{V}^{(e)}}-\Delta t{\displaystyle {\iiint }_{V^{(e)}}{u}_i^n\frac{\partial N_{\alpha }^{(e)}}{\partial x_i}{T}_{\alpha }^{n(e)}{N}_{\beta }^{(e)}d{V}^{(e)}}\\ -\Delta t\alpha {\displaystyle {\iiint }_{V^{(e)}}\frac{\partial N_{\alpha }^{(e)}}{\partial x_i}{T}_{\alpha }^{n(e)}\frac{\partial N_{\beta }^{(e)}}{\partial x_i}d{V}^{(e)}}-\frac{\Delta t^2}{2}{\displaystyle {\iiint }_{V^{(e)}}{u}_i^n\frac{\partial N_{\alpha }^{(e)}}{\partial x_i}{T}_{\alpha }^{n(e)}{u}_j^n\frac{\partial N_{\beta }^{(e)}}{\partial x_j}d{V}^{(e)}}\end{array}$$

(46)

The total finite-element equation can be obtained by synthesis:

$$A_q{T}_q^{n+1}=B_q$$

(47)

$q$ is the overall node number of the temperature.

The temperature value $T^{n+1}$ of the energy-conservation equation (Equation (10)) at time $n+1$ can be solved by Equation (47).

3.3. Algorithm Procedure

Assuming that the value of the air in the tunnel at a certain time (denoted by n) is known, the calculation process for solving the thermal convection control equation of tunnels in cold regions is shown as Figure 1 and Figure 2 when the value of the next time (denoted by n + 1) is obtained.

4. Numerical Solution of Heat Conduction-Heat Convection Fluid-Structure Interaction Model of Tunnel in Cold Region

4.1. Model Introduction

In the case of natural ventilation, the temperature field distribution of the tunnel is mainly affected by the internal temperature of the surrounding rock, the temperature in the tunnel, and the change in the external environment temperature. Therefore, the fluid–structure interaction model of tunnels in cold regions is the mutual transfer model of the surrounding rock–lining–air heat. The fluid–structure interaction model of tunnels in cold regions mainly includes:

(1)

Tunnel surrounding rock and lining heat transfer

This mainly calculates the influence of the external environment on the temperature of the tunnel surrounding rock and the temperature transfer between the surrounding rock and the lining.

(2)

Thermal convection of air in tunnel

Firstly, the air velocity field distribution in the tunnel under the inlet wind speed is calculated, and then the air temperature distribution in the tunnel is calculated by regarding the external environment temperature as the boundary condition.

(3)

Heat transfer between tunnel air and lining

Convective heat transfer will occur between the inner surface of the tunnel lining and the air. The finite model of thermal convection and its solution method are proposed above. Through this model, the temperature distribution of air in the tunnel under different wind speeds and wind directions can be calculated, which in turn affects the temperature propagation of the tunnel lining.

In the process of heat transfer between the tunnel air and the lining, the air temperature in contact with the lining at the same time will affect the surface temperature of the lining, and the surface temperature of the lining will also affect the temperature of the air in contact with it. Therefore, the following method is adopted. The temperature of the lining surface at this time is taken as the ambient temperature of the air boundary in the tunnel, and the air temperature in contact with the lining surface at the next time is calculated by using the third boundary condition. Then, the air temperature is used as the ambient temperature of the lining boundary, and the temperature of the lining surface at the next moment is calculated with the help of the third boundary condition.

4.2. Algorithm Procedure

The fluid–structure interaction model of tunnels in cold regions includes heat transfer between the tunnel surrounding rock and lining, air heat convection in the tunnel, heat transfer the between tunnel air and lining, etc. The solution process is shown in Figure 3:

5. Application Example

5.1. Project Profile

The Hekashan tunnel is located in Xinghai County, Hainan Tibetan Autonomous Prefecture, Qinghai Province. It is a straight tunnel with a total length of 2300 m. The starting and ending pile numbers are K232 + 090~K234 + 390, with an average altitude of 3795 m and a maximum buried depth of 210 m. It adopts a composite lining structure.

The annual freezing period of the tunnel site lasted from October to April of the following year, and the entrance and exit of the tunnel were permafrost sections. January is the coldest month of the year, with an average temperature of −9.9 °C and a minimum temperature of −29 °C. The highest temperature is in July, with a maximum temperature of 16.0 °C. The annual average wind speed is about 3.5 m/s. Due to the large temperature difference between day and night in the tunnel area, the freezing and thawing frost heave phenomenon of the tunnel is extremely serious.

The tunnel site is located in Heka Mountain, Xinghai County, which is mainly composed of Lower Triassic feldspar sandstone, sandstone, calcareous slate, silty clay slate, and Quaternary Upper Pleistocene alluvial and diluvial deposits. The grade of the surrounding rock is grade III~V, in which the rock mass at the inlet and outlet is broken into gravel, and the integrity is poor, all of which are grade V surrounding rock. The average temperature of the surrounding rock is 0.7 °C.

5.2. Monitoring Scheme and Result Analysis

The first thing to be determined is the range of monitoring length, using the Kurokawa Xifan formula. It is concluded that the length range of the tunnel monitoring section is 618 m. According to the recommended length of the insulation layer of the Zhangzhidao tunnel (

Table 1. Suggested length table for laying insulation section of long tunnel in cold region.

Elevation above Sea Level (m)	Average Temperature of Coldest Month (°C)	Length of Insulation Section (m)
3300	−10	680
3600	−10.5	690
3800	−11	710
4000	−12	750
4200	−13	830
4400	−14	860
4600	−15	900
4800	−16	930

), combined with the actual situation of the site, the length range of the monitoring and insulation section of the Hekashan tunnel is determined to be 680 m [37].

The monitoring contents include air temperature at the tunnel entrance, wind speed perpendicular to the tunnel entrance, and surface temperature of the inner wall of the tunnel lining.

Two fixed monitoring points are set up on the left and right sides of the entrance of the Hekashan tunnel. The temperature and wind speed of each monitoring point were measured at 07:30 and 14:30 every day, and the arithmetic mean of the test results of the two monitoring points was taken as the average temperature and wind speed of the day.

Five observation sections of 0 m, 200 m, 400 m, 600 m, and 680 m are set up in the depth direction of the Hekashan tunnel from 0 m to 680 m. Five monitoring points of the left and right arch foot, the left and right arch waist, and the vault of the secondary lining are determined for each observation section. The longitudinal and section measurement layout of the tunnel is shown in Figure 4.

The temperature of the monitoring points was measured at 07:30 and 14:30 daily, and the average value of the two test results was taken as the average temperature of the day.

The measurement date of air temperature, wind speed, and lining wall temperature of the Kashan tunnel of the Tuohe river starts from 1 January 2016 to 31 December 2016. The monitoring results are as follows:

(1)

Temperature change

By fitting the daily average temperature of 365 days, the fitting formula of the temperature change with time at the tunnel entrance is as follows:

$$T_a=-6.5+22.5\sin (\frac{2\pi }{365}t+\frac{3\pi }{2})$$

(48)

T_a—Daily average temperature, °C; t—time, day.

Figure 5 shows the temperature change curve at the entrance of the Hekashan tunnel in 2016, which is a sine curve.

(2)

Wind speed variation

Through the test and calculation of the wind speed of the Hekashan tunnel, the average wind speed at the entrance of the tunnel in winter is about 3.8 m/s, the average wind speed at the entrance of the tunnel in summer is 3.1 m/s smaller than that in winter, and the wind speed at other times is between 3.2 and 3.5 m/s. The one-year wind speed is 3.47 m/s by using the arithmetic average method.

(3)

Lining wall temperature change

The arithmetic mean value of the temperature measured at five monitoring points in each section of the tunnel every day is taken as the daily average temperature of the secondary lining surface. The five sections of the tunnel are arranged and plotted, as shown in Figure 6.

5.3. Calculation Model

The numerical calculation model is established by taking the length of the right line of the Hekashan tunnel from the entrance to the 680 m from the entrance. The right line of the tunnel is one-way two-lane, and the surrounding rock grade of the entrance section is V grade. The initial lining is made of C25 shotcrete with a thickness of 26 cm. The secondary lining is C30 reinforced concrete with a thickness of 50 cm. The net height of the tunnel lining is 6.60 m, the net width is 11.10 m, and the cross-section adopts a four-center circle.

The influence of steel bars on the lining structure on heat transfer is ignored, and it is assumed the surrounding rock, lining, and air in the tunnel are homogeneous and isotropic continuous media. As shown in Figure 7, the model takes the tunnel lining contour to the upper, lower, left, and right surrounding rock boundaries of 40 m each. A three-dimensional calculation model is established, and the eight-node hexahedral element is used to divide the surrounding rock, primary lining, secondary lining, and air.

There are many factors affecting the distribution of the air velocity field in tunnels in cold regions, including tunnel plane alignment, cross-section shape, initial wind speed, lining wall roughness, air temperature in the tunnel, and boundary conditions. For the inlet and outlet of short tunnels or long tunnels in cold regions, due to its short length and smooth wall surface, this paper ignores the air resistance loss, dynamic pressure loss, and local resistance loss in order to simplify the calculation.

The secondary lining is divided into three layers along the radial direction of the tunnel, and the primary lining is one layer. As shown in Figure 8, the temperature variation of Layers A, B, C, D, and E in the secondary lining and the primary lining of the tunnel is compared and analyzed. The radial distances from each layer to the inner surface of the secondary lining are 0 cm, 17 cm, 33 cm, 50 cm, and 76 cm.

Along the depth direction of the tunnel, the temperature variation of the five sections of 0 m, 200 m, 400 m, 600 m, and 680 m is mainly studied and analyzed, as shown in Figure 9.

Through on-site sampling and investigation, the thermal parameters of the surrounding rock, lining concrete, and air in the tunnel located in the frozen soil section are summarized in

Table 2. Calculation parameter table of fluid–structure interaction model of Hekashan tunnel.

Parameter	Surrounding Rock	Primary Lining	Lining Concrete	Air in the Hole
Thermal conductivity λ (W/(m·°C))	3.50	1.70	1.85	2.30 × 10−2
Density ρ (kg/m3)	2120	2300	2500	1.40
Specific Heat Capacity c (J/(kg·°C))	877	950	970	1000
Dynamic Viscosity μ (Pa/·s)	-	-	-	1.82 × 10−5

at the section of K232 + 120 pile at the right entrance of the Hekashan tunnel.

Initial conditions: Assuming that the initial velocity, pressure, and initial temperature of air in the tunnel are 0, the initial temperature of the surrounding rock and lining is measured to be −0.7 °C. For the boundary conditions, the wind speed in the direction perpendicular to the tunnel entrance is 3.5 m/s, and the speed in the other two directions of the tunnel entrance is assumed to be 0. The contact part between the inner wall surface of the tunnel lining and the air is the third boundary condition.

5.4. Verification and Analysis of Calculation Results

According to the field monitoring results, the ambient temperature of the Hekashan tunnel reached its highest value around 1 July 2016 (Day 183) and reached its lowest value on 31 December 2016 (Day 365). The overall temperature field and profile distribution results are shown in 4.7–4.8. It can be seen from the figure that the temperature at the bottom of the surrounding rock is basically about 3 °C in the range of 10 m upward, indicating that the temperature at the bottom of 3 °C has little influence on the surrounding rock. The temperature field distribution of the whole and the section is shown in Figure 10 and Figure 11.

The temperature-calculated values of the inner surface of the secondary lining of different sections (Layer A in Figure 8) are compared with the field-measured values, as shown in Figure 12. It can be seen from the diagram that the calculated value is basically consistent with the trend of the measured value over time.

Through comparative analysis, there is a certain difference between the calculated value and the field-measured value. The reason is that in the calculation of the air velocity field, the coupling model ignores the influence of air resistance loss, dynamic pressure loss, and local resistance loss on the flow field. In this paper, the temperature propagation law at the 0 m position of the tunnel entrance is studied and analyzed. The temperature changes of each layer of lining A, B, C, D, and E are shown in Figure 13.

It can be seen from the figure that (1) the lining temperature of each layer of the tunnel changes in a sine curve with time, and its change law is similar to the change law of the external environment temperature. (2) The average temperature difference between the lining layers is about 2 °C. (3) It can be seen from the temperature change curve that the temperature of each layer has a gradual lag.

The temperature extreme value of each layer of lining and its occurrence time are shown in

Table 3. The temperature extreme value and its occurrence time of each layer of tunnel portal lining.

Research Location	Maximum Temperature (°C)	Occurring Time (d)	Minimum Temperature (°C)	Occurring Time (d)
Secondary lining A layer	14.90	3 July	−27.63	31 Dec.
Secondary lining B layer	13.52	6 July	−25.85	31 Dec.
Secondary lining C layer	12.00	11 July	−23.77	31 Dec.
Secondary lining D layer	10.85	14 July	−22.13	31 Dec.
Secondary lining E layer	9.09	21 July	−19.48	31 Dec.

. It can be seen from the table that (1) the temperature variation amplitude of the temperature extreme value of each layer gradually decreases with the increase in the radial distance of the lining. (2) With the increase in the radial distance of the lining, the occurrence time of the temperature extreme value of each layer gradually lags behind, and the lag time is 3–7 days.

Figure 14 shows the variation law of the temperature of each lining layer with the depth of the Hekashan tunnel at the highest temperature of 16 °C and the lowest temperature of −29 °C. It can be seen from the diagram that (1) when the inlet temperature is 16 °C, the temperature of each layer of the tunnel lining remains unchanged in the first 500 m, and then the temperature gradually decreases with the increase in the depth. (2) When the inlet temperature is −29 °C, the temperature of each layer of the tunnel lining remains unchanged at the first 500 m, and then the temperature gradually increases with the increase in the depth.

6. Conclusions

This paper is based on the heat transfer mechanism of the tunnel surrounding rock–lining–air. The coupling calculation model of air heat convection in cold tunnels is combined with the tunnel heat transfer control model, and the fluid–structure interaction finite-element model is then established. Taking the Hekashan tunnel as an example, the temperature field of the tunnel portal section in the cold region is studied. The main conclusions are as follows:

(1)

In terms of the idea of the characteristic-based operator-splitting (CBOS) finite-element method and the explicit characteristic–Galerkin method, a finite-element method for solving three-dimensional N-S equations and thermal convection governing equations is developed., and the velocity and temperature are coupled. A three-dimensional finite-element calculation model of air thermal convection in tunnels in cold regions is established. After the velocity and pressure are decoupled in the model, the same order interpolation function can be used, which greatly improves the calculation efficiency. Consequently, the calculation results can be of high accuracy.

(2)

Considering the dynamic influence of air velocity and temperature on the lining surface in the tunnel, the fluid–structure interaction finite-element calculation model of tunnels in cold regions is established by coupling the heat transfer model of the tunnel lining and surrounding rock with the air heat convection model. Compared with other calculation models, this model can fit the actual situation and more accurately reflect the distribution law of the lining temperature field and air temperature field in each layer of the tunnel.

(3)

The fluid–structure interaction calculation model of tunnels in cold regions is used to simulate the temperature field of 680 m at the portal section of the Hekashan tunnel and then compared with the measured data. It is found that the calculated value is basically consistent with the measured value over time, which indicates that the fluid–structure interaction model established in this paper has certain reliability.

(4)

Temperature monitoring points were set up in the sections of 0 m, 200 m, 400 m, 600 m, and 680 m of the secondary lining of the Hekashan tunnel. The measured temperature curves of each section of the secondary lining of the tunnel were obtained by averaging the five monitoring data of the left and right arch foot, the left and right arch waist, and the vault. It is found that the temperature of each layer of the tunnel lining changes in a sine curve with time, and the temperature of each layer gradually lags behind. The temperature variation amplitude of the extreme value of the layer temperature gradually decreases with the increase in the radial distance of the lining. Near the entrance end of the tunnel, the temperature of each layer of lining remains unchanged, and the temperature gradually decreases or increases with the increase in the depth.

(5)

In this paper, the distribution law of the temperature field in tunnels in cold areas is studied, on the basis of which the stress–strain relationship between the surrounding rock and the lining of the tunnel can be further studied.