数学建模——永冻土层上关于路基热传导问题

问题描述

在永冻土地上铺设道路、飞机跑道和某些结构的地基。分析这类地基的结构，有沥青层、混凝土层、干砂层、石子层、绝缘材料层，再下面是湿的沙土层和冻土层（参见图1）。已知，外界的温度是关于时间的函数，冻土层的温度是不变的零下温度。请回答以下问题：

• 问题1：分析空气温度传入路基规律，各材料层的温度分布；
• 问题2：由于一些设备的支架不能固定在解冻土层上，必须固定在永冻土层中，因此确定地下土层的解冻位置非常重要，请给出解冻砂土与冻结砂土的分界线；
• 问题3：请给出各层材料的最佳厚度（结合温度分布、成本和耐用性 ）；
• 问题4：请结合我国青藏铁路永冻土层地基进行仿真，并给施工单位给出合理建议。

模型求解

问题一的分析

对于问题一，考虑到冻土面上各点处曲率足够大，将地基层视为互相平行的平面层，建立基于傅里叶定律的一维抛物型偏微分方程模型，由于外界的温度是关于时间的函数，因此可以通过搜集拉萨四季的日温度数据，进行正弦函数拟合，从而得到外界温度随时间变化的拟合函数。将其作为热传导方程的边界条件，根据查阅冻土土层温度为0℃及以下，本组选取 $T_0=0$ ℃ 作为热传导方程右边界，在给定一个初值条件下，通过有限差分法，给出精度较高的Crank-Nicolson差分格式，用微商代替微分得到离散的线性方程组，通过追赶法求解线性方程组，得到时间和一维空间维度下的温度分布情况，并对结果进行稳定性分析以及误差检验。

一维耦合抛物型偏微分方程模

热传导方程

根据傅里叶热传导定律得到,导体获取热量微元等式：
$$dQ=-\lambda \nabla u\cdot d\overrightarrow{S}dt$$
考虑问题研究为一维，通过积分以及高斯公式进行化简后得到$t_1$~$t_2$时间段内流入导体热量为：
$$Q_1={\int }_{t_1}^{t_2}\left[{\int }_{x_1}^{x_2}k\frac{{\partial }^{2}u}{\partial x^2}dx\right]dt$$
又由于导体内温度升高吸收热量为：
$$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$$
根据能量守恒得到能量关系等式：
$$Q_1=Q_2$$
从而得到：
$$\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}$$
其中：$k$为热传导率，$c$为导体比热，$\rho$为导体干密度。
微元法示意图如图2。

耦合抛物型偏微分方程

考虑到5层不同地基结构的热传导率、比热以及干密度的不同，建立五个一维热传导方程：
$$\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}$$
根据Fourier热传导定律，在导热过程中，两相邻地基材料临界面处热流量密度$\vec{q}$相同，得到热流量密度耦合条件:
$$\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}$$
由任意时刻下相邻导体临界面温度相等，从而得到导体传热的临界面温度耦合条件：
$$\begin{array}{r}{u}_i\left(x,t\right){|}_{{\Gamma }_i}=u_{i+1}\left(x,t\right){|}_{{\Gamma }_i},\end{array}$$
通过对搜索的拉萨四季日温度数据进行正弦函数二阶拟合得到:
$$\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}$$
由此通过拟合函数作为外界温度的变化函数作为热传导方程的边界条件:
$$\begin{array}{r}{u}_0\left(0,t\right)=a\left(t\right),\end{array}$$
根据查阅得知冻土层温度为0℃及以下，本组选取$T_0=0$，作为右边界条件:
$$\begin{array}{r}{u}_5\left(L_5,t\right)=T_0,\end{array}$$
最终整理上述公式得到初边值问题的热传导方程：
$$\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}$$

采用CN格式处理

$$\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}$$

问题一结果展示

后续问题

。。。