Mathematical modeling -- heat conduction of subgrade on Permafrost

Problem description

Pave roads on permafrost 、 Foundation of aircraft runway and some structures . Analyze the structure of this kind of foundation , There is an asphalt layer 、 Concrete layer 、 Dry sand layer 、 Gravel layer 、 Insulating material layer , Below is the wet sandy soil layer and frozen soil layer （ See chart 1）. It is known that , The outside temperature is a function of time , The temperature of the frozen soil layer is constant below zero . Please answer the following questions ：

• problem 1： Analyze the law of air temperature entering the subgrade , Temperature distribution of each material layer ;
• problem 2： Because the supports of some equipment cannot be fixed on the thawed soil , It must be fixed in permafrost , Therefore, it is very important to determine the thawing location of the underground soil , Please give the dividing line between thawed sand and frozen sand ;
• problem 3： Please give the best thickness of each layer （ Combined with temperature distribution 、 Cost and durability ）;
• problem 4： Please simulate the permafrost foundation of Qinghai Tibet Railway in China , And give reasonable suggestions to the construction unit .

Model solving

Analysis of problem 1

For question 1 , Considering that the curvature at each point on the frozen soil surface is large enough , Treat the foundation layer as a plane layer parallel to each other , A one-dimensional parabolic partial differential equation model based on Fourier law is established , Because the outside temperature is a function of time , Therefore, we can collect the daily temperature data of the four seasons in Lhasa , Sine function fitting , Thus, the fitting function of the external temperature changing with time is obtained . Take it as the boundary condition of the heat conduction equation , According to the reference, the temperature of frozen soil layer is 0℃ And the following , Selected by this group $T_0=0$ ℃ As the right boundary of the heat conduction equation , Given an initial value , By finite difference method , Give a high precision Crank-Nicolson The difference scheme , The discrete system of linear equations is obtained by replacing differential with derivative , Solving linear equations by catch-up method , Get the temperature distribution in time and one-dimensional space dimension , And the stability analysis and error test of the results .

One dimensional coupled parabolic partial differential equation module

Heat conduction equation

According to Fourier heat conduction law , The differential equation of the conductor acquiring heat ：
$$dQ=-\lambda \nabla u\cdot d\overrightarrow{S}dt$$
Consider that the problem is studied in one dimension , After simplification by integral and Gauss formula $t_1$~$t_2$ The heat flowing into the conductor during the time period is ：
$$Q_1={\int }_{t_1}^{t_2}\left[{\int }_{x_1}^{x_2}k\frac{{\partial }^{2}u}{\partial x^2}dx\right]dt$$
Due to the increase of temperature in the conductor, the absorbed heat is ：
$$Q_2={\int }_{x_1}^{x_2}c\rho \left[{\int }_{t_1}^{t_2}\frac{\partial u}{\partial t}dt\right]dx={\int }_{t_1}^{t_2}\left[{\int }_{t_1}^{t_2}c\rho \frac{\partial u}{\partial x}dx\right]dt$$
According to the conservation of energy, the energy relationship equation ：
$$Q_1=Q_2$$
To get ：
$$\frac{\partial u\left(x,t\right)}{\partial t}=\frac{k}{c\rho }\frac{{\partial }^2{u}_i\left(x,t\right)}{\partial x^2}$$
among ：$k$ Is the thermal conductivity ,$c$ Is the specific heat of conductor ,$\rho$ Is the dry density of the conductor .
The schematic diagram of microelement method is shown in Figure 2.

Coupled Parabolic Partial Differential Equations

in consideration of 5 Thermal conductivity of different foundation structures 、 Differences in specific heat and dry density , Establish five one-dimensional heat conduction equations ：
$$\begin{array}{r}\frac{\partial u_i\left(x,t\right)}{\partial t}=\frac{k_i}{c_i{\rho }_i}\frac{{\partial }^2{u}_i\left(x,t\right)}{\partial x^2},i=1\cdots 5,\end{array}$$
according to Fourier Law of heat conduction , In the process of heat conduction , Heat flux density at the critical interface of two adjacent foundation materials $\vec{q}$ identical , The coupling condition of heat flux density is obtained :
$$\begin{array}{r}{k}_i\frac{\partial u_i\left(x,t\right)}{\partial n_i}{|}_{{\Gamma }_i}=k_{i+1}\frac{\partial u_{i+1}\left(x,t\right)}{\partial n_{i+1}}{|}_{{\Gamma }_i},i=1\cdots 4,\end{array}$$
The temperature at the critical interface of adjacent conductors is equal at any time , Thus, the critical interface temperature coupling condition of conductor heat transfer is obtained ：
$$\begin{array}{r}{u}_i\left(x,t\right){|}_{{\Gamma }_i}=u_{i+1}\left(x,t\right){|}_{{\Gamma }_i},\end{array}$$
Through the sine function second-order fitting of the searched temperature data of the four seasons in Lhasa :
$$\begin{array}{r}a(t)={\alpha }_i\sin \left(w_{1}t+{\phi }_i\right)+{\beta }_i\sin \left(w_{2}t+{\theta }_i\right),i=1\cdots 4,\end{array}$$
Therefore, the fitting function is used as the change function of the external temperature as the boundary condition of the heat conduction equation :
$$\begin{array}{r}{u}_0\left(0,t\right)=a\left(t\right),\end{array}$$
According to the review, the temperature of frozen soil layer is 0℃ And the following , Selected by this group $T_0=0$, As the right boundary condition :
$$\begin{array}{r}{u}_5\left(L_5,t\right)=T_0,\end{array}$$
Finally, the heat conduction equation of the initial boundary value problem is obtained by sorting out the above formula ：
$$\begin{array}{r}\{\begin{array}{l}\frac{\partial u_i\left(x,t\right)}{\partial t}=\frac{k_i}{c_i{\rho }_i}\frac{{\partial }^2{u}_i\left(x,t\right)}{\partial x^2},x\in \left(L_{i-1},L_i\right),\\ k_i\frac{\partial u_i\left(x,t\right)}{\partial n_i}{|}_{{\Gamma }_i}=k_{i+1}\frac{\partial u_{i+1}\left(x,t\right)}{\partial n_{i+1}}{|}_{{\Gamma }_i},\\ u_i\left(x,t\right){|}_{{\Gamma }_i}=u_{i+1}\left(x,t\right){|}_{{\Gamma }_i},\\ u_0\left(0,t\right)=a\left(t\right),\\ u_5\left(L_5,t\right)=T_0,\end{array}\end{array}$$

use CN Format processing

$$\begin{array}{r}\{\begin{array}{l}\frac{u_j^{n+1}-u_j^n}{\tau }=\frac{k_i}{c_i{\rho }_i}\cdot \frac{\left[\left(u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}\right)+\left(u_{j+1}^n-2u_j^n+u_{j-1}^n\right)\right]}{2h^2}\\ k_i\frac{u_J^n-u_{J-1}^n}{h}=k_{i+1}\frac{u_{J+1}^n-u_J^n}{h}\\ u_i\left(J,t\right)=u_{i+1}\left(J,t\right)\\ u_1\left(0,t\right)=a\left(t\right)\\ u_5\left(N,t\right)=T_0\\ J=\frac{L_i}{L}N\end{array}\end{array}$$

Question 1 result display

Follow-up questions

...