?2?u2?x?u3?x=?3?2?u3?x表示介质II与介质III的两个耦合条件； ，u2?=u3? ，
?222?3=?4?3?u4?x表示介质III与介质IV的两个耦合条件； ，u3?=u4? ，
?333ui(x,0)=37 ， i=1,2,3,4 ，表示热传导方程的初值条件； u1(0,t)=75 ，表示热传导方程的左边界Dirichlet边值条件；
?u4=h(uh?u5) ，表示热传导方程的右边界Robin边值条件。 ?x5.1.5 基于热传导方程的温度分布模型的求解 （1）求解方法分析
在5.1.4中建立的温度分布模型（20）属于抛物型方程，因边值条件条件复杂难以求得解析解。本文采用有限差分法：将连续的定解区域用有限个离散点构成的网格来代替；把定解区域上的连续变量的函数用网格上定义的离散变量的函数来近似；把原方程和定解条件中的微商用差商来近似。最终，把原微分方程和定解条件用代数方程组来代替，即有限差分方程组。解此方程组就可以得到原问题在离散点上的近似值。
常用的差分格式有：向前差分格式、向后差分格式和C-N差分格式（即Crank-Nicolson差分格式）。
在向前差分格式中，要求时间步长与空间步长的比r?0.5 ，即时间步长要比空间步长小得多。若不满足此条件，极易发生解的爆破。在本文温度分布模型中，时间步长与空间步长达到50，远远大于0.5，因此不适合用向前差分格式。此外，向前差分格式为显式格式，难以进行多层介质的耦合。
在C-N差分格式中，稳定性好精度高，计算在交叉点处的函数值，具体如下图7所示：
??4 图7 C-N差分格式图
在向后差分格式中，具有无条件稳定的特点，选取第j+1时间层相邻的三个结点进行u对x二阶偏导的近似，选取第j+1时间层与第j时间层的相邻两个结点进行对u对t一阶偏导的近似，具体如下图8所示：
13
图8 向后差分格式图
（2）温度分布模型的具体求解过程：
经过上述三种差分格式的比较，本文采用向后差分格式求解，并借助迭代法求解。求解流程如下图9所示：
图9 温度分布模型求解流程图
Step1：对求解区域进行网格剖分。
求解区域记为? ?=?(x,t)0?x?L,0?t?T?，其中L=?Li且Li为第i层介
i=14质的厚度。对此求解区域进行网格剖分，具体如下：
① 在空间维度上，将区间[0,L1]做M1等分，将区间[L1,L2]做M2等分，将区
间[L2,L3]做M3等分，将区间[L3,L4]做M4等分。
i=1,2,3,4 。 记mi=?Ml ，l=1i记?xi=Li ，?xi 表示为第i层介质的空间步长。 Mi14
② 在时间维度上，将区间[0,T]作n等分，并记?t=T，?t为时间步长。 n??x=xi,1?i?m③ 用两组平行直线簇? 将?分割成矩形网格，见下图10：
t=t,1?j?n?j? 图10 网格剖分图
Step2：建立隐式向后差分格式。
定义?上的网格函数，u=uij1?i?m4+1,1?j?n+1，其中uij=u(xi,tj)，且有1?i?m4+1,1?j?n+1 。
考虑温度分布模型（20）中的热传导方程中的右端项，略去小量项，用二阶
中心差商代替u对x的偏导数，得到：
??
?2u1(x,t)=[u(xi?1,tj)?2u(xi,tj)+u(xi+1,tj)] ij?x2(?x)2(23)
考虑温度分布模型（20）中的热传导方程中的左端项，用一阶向前差商代替
u对t的一阶偏导数，得到：
?u1(xi,tj)=[u(xi,tj+1)?u(xi,tj)] ?t?t(24)
建立差分格式的具体步骤如下：
① 在结点处考虑不同介质下的热传导方程，得到介质内部的差分格式。
2?uk2?uk=ak , k=1,2,3,4 ，则相应的差分格第i层介质的热传导方程为：?t?x2式为：
i+1iukj?ukji+1i+1i+1uk,j?1?2ukj+uk,j+1
?t?a?2k(?x)15
2=0 (25)
i其中，ukj 表示第k层介质在(xi,tj)结点处的温度，k=1,2,3,4 ; j=1,2, ,m。
② 依据热流量密度的耦合条件，得到介质边界处的差分格式。 在临界面?1、?2、?3 处，关于热流量密度的耦合条件为：
?k?uk?x=?k+1?k?uk+1?x,k=1,2,3
?k则得到相应的差分格式为：
i+1iuk,mk+1?uk,mki+1i+1uk+1,mk+2?uk+1,mk+1
?k?xk=?k+1??xk+1 , k=1,2,3 (26)
③ 依据第IV层介质右侧Robin边值条件，得到边界?4 上的差分格式。 第IV层介质边界?4处，满足温度分布模型（20）的条件??4可得到相应的差分格式为：
i+1i+1um?um44?1?u4=h(uh?u5) ， ?x
??4?x4i+1=h(um?37) 4(27)
综合①②③，得到温度分布模型（20）的隐式向后有限差分近似如下：
i+1ii+1i+1i+1?ukj?ukju?2u+ukjk,j+1?ak2?k,j?1=0 k=1,2,3,4 , j=1,2,,m?2?t(?x)?i+1i+1?ui+1?uiu?u??kk,mk+1k,mk=?k+1?k+1,mk+2k+1,mk+1 k=1,2,3?xk?xk+1??i+1i+1? (28) ???um4?um4?1=h(ui+1?37)m4?4?x4??u1i+1=75 i=1,2,,n?1?ukj=37 k=1,2,3,4 , j=1,2,,m+1?jju=u,n+1 k,m+1k+1,mk+1 k=1,2,3 , j=1,2,??k其中，（26）式中第一式、二式及三式在步骤①、步骤②及步骤③中已做详细说明；（26）式第四式，为基于热传导方程的温度分布模型（20）中左边界Dirichlet边值条件；（26）式第五式，为温度分布模型（20）中初值条件；（26）式第六式，为温度分布模型（20）中基于温度的耦合条件条件。
Step3：将差分格式整理为代数方程组。 ①差分格式（26）中，第一式可整理为：
16
jj?1?rkukj,i?1+(1+2r)uki?rukj,i+1=uki
(29)
其中，k=1,2,3,4 , 2?i?m , 2?j?m 。
rk=ak?tk ，表示为第k层介质剖分的步长比。 2(?x)②差分格式（26）式中，第二式可整理为：
??k?xki+1uk,mk+(?k?xk+?k+1?xk+1i+1)uk,mk+1??k+1?xk+1i+1uk+1,mk+2=0
(30)
③差分格式（26）式中，第三式可整理为：
??4?x4i+1u4,m4?1+(?4?x4i+1+h)u4,m4=37h
(31)
由①②③可将差分格式（26）写成如下矩阵形式：
?1+2r1?r1??r1+2r?r11?1???r11+2r1?r1???1?1?2?1?+???x1?x1?x2?x1????3?3?4?3??+???x3?x3?x4?x3??r41+2r4?r4?????r41+2r4????4??x4?i+1ii+1?u12??u12?+ru111??i+1???i+1??u13??u13??????????????????i+1??i+1???u1m1??u1m1???ui+1??0???1m1+1????????=?????????i+1?i+1??u3,u??3,m3+1?3+1??i+m??i+1?1u??u4,m4,m+13??3+2????????????r4???????i+1??i+1??4uu?4,m?14,m?14+h?4???i+1?x4?i+1??u??u4,m?4??4,m4??
Step4：利用追赶法解三对角线性方程组
在step3中，得到的代数矩阵为三对角线性方程组AU=b，一般可用matlab求逆命令求解U，但因本题中系数矩阵规模较大，且除主对角线和两个次对角线之外，其余元素均为0，因而求逆解方程组会导致浪费内存、程序运行缓慢。因此本文利用追赶法求解三对角线性方程组。
17
2024年全国大学生数学建模大赛A题高温作业专用服装设计及三篇优秀论文精选 - 图文
