8.9 二维非瞬态热传导问题的差分格式
在直角坐标系内二维热传导问题可以用如下方程描述：
其中k=k（x，y）为热传导系数，c=c（x，y）为物体的比热容，ρ=ρ（x，y）为物体的密度，f（x，y，t）为热源。
如果物体均匀且各向同性，即k、c和ρ都为常数时方程（8.47）变为
其中
当f（x，y，t）=0时，方程（8.48）中的
用对t偏导数的向前差分即〔式（8.34）〕来表示
用对x二阶偏导数的中心差分式〔式（8.42）〕来表示
用对y二阶偏导数的中心差分〔式（8.46）〕来表示，代入式（8.48）并整理，得
其中
。
例8.4 在方程（8.48）中，当f（x，y，t）=0时，初始条件和边界条件分别为T（x，y，0）=sin（πx）+sin（πy）和T（x，y，0）
=0（Ω=[0，1]×[0，1]），则式（8.48）的解析解为
第一步，调试迭代计算。在“文件”菜单下选择“选项”，在“选项”对话框中选择“公式”副对话框。选择启用迭代计算，最多迭代计算次数设为1，最大误差设为0.001并单击“确定”。
第二步，输入已知数据。在单元格A2、B2、C2和D2中分别输入常数α，区间的长度，空间x、y方向和时间t方向的网格步长参数（在这里分别输入“1”“1”“0.1”和“0.001”）；在单元格E2、F2和G2中分别输入常数
、计算停止时间的公式和计算总计算次数的公式（这里分别输入“=D2/（C2^2）”“0.5”和“=F2/D2”）；在单元格H2中输入实现循环的表达式“=IF（H2＜G2/2，H2+1，0）”；在单元格I2中输入计算当前计算时刻的表达式“=2*H2*D2”，如图8.16所示。
图8.16 已知数据
第三步，输入初始条件。在单元格区域A5：A15中输入自变量y的值；在单元格区域B16：L16中输入自变量x的值；在单元格B5中输入初始条件T（x，y，0）=sin（πx）+sin（πy）的计算公式“=SIN（PI（）*$A5）+SIN（PI（）*B$16）”。注意，单元格A5对应自变量y，其应用列绝对引用，行相对引用；单元格B16对应自变量x，其应用行绝对引用，列相对引用。选中单元格B5，把鼠标指针移动到单元格B5的右下角，待其变为填充柄（“
”符号）时，按住鼠标左键向右拖动到单元格L5，完成单元格区域B5：L5的自动填充。选中单元格区域B5：L5，把鼠标指针移动到单元格L5的右下角，待其变为填充柄（“
”符号）时，按下鼠标左键向下拖动到单元格L45，实现单元格区域B5：L15的自动填充，完成初始条件的输入，如图8.17所示。
图8.17 初始时刻的值(https://www.daowen.com)
第四步，从初始时刻计算第一时刻。和第三步一样，在单元格区域A19：A29中输入自变量y的值；在单元格区域B30：L30中输入自变量x的值；在单元格B19中输入计算边界条件表达式“=（SIN（PI（）*$A19）+SIN（PI（）*B$30））*EXP（-$A$2*PI（）*PI（）*（$I$2-$D$2））”；用自动填充功能拖动填充单元格区域C19：L19、L20：L29、B20：B29和C29：K29；用式（8.40）在单元格C20中输入“=C6+$A$2*$E$2*（C5+C7+B6+D6-4*C6）”；用自动填充功能拖动填充单元格区域C20：K28，这样就得到第一时刻的数值解，如图8.18所示。
图8.18 第一时刻的数值解
第五步，从第一时刻计算第二时刻。和第三步一样，在单元格区域A33：A43中输入自变量y的值；在单元格区域B44：L44中输入自变量x的值；在单元格B19中输入计算边界条件表达式“=（SIN（PI（）*$A47）+SIN（PI（）*B$58））*EXP（-$A$2*PI（）*PI（）*$I$2）”；用自动填充功能拖动填充单元格区域C33：L33、L34：L43、B34：B43和C43：K43；用式（8.49）在单元格C20中输入“=C20+$A$2*$E$2*（C19+C21+B20+D20-4*C20）”；用自动填充功能拖动填充单元格区域C34：K42，这样就得到第二时刻的数值解，如图8.19所示。
图8.19 第二时刻的数值解
第六步，显示正在计算时刻，将单元格区域A18：L18合并为一个单元格并输入“=＂n=＂&2*H2-1&＂，t=＂&I2-D2&＂时刻的数值解＂”；将单元格区域A32：L32合并为一个单元格并输入“=＂n=＂&2*H2&＂，t=＂&I2&＂时刻的数值解＂”。
第七步，实现循环计算。为了实现循环计算，将单元格C20中的计算公式改为“=IF（$H$2=1，C6+$A$2*$E$2*（C5+C7+B6+D6-4*C6），C34+$A$2*$E$2*（C33+C35+B34+D34-4*C34））”，用自动填充功能拖动填充单元格区域C20：K28。这时，如果我们按功能键F9，则在单元格H2的值每一次增加1，循环在1～250之间，当H2的值不等于1时，单元格区域C20：K28用单元格区域C34：K42的值进行计算。而且单元格区域C34：K42的值随着单元格区域C20：K28值的更新而更新，这样在单元格区域C34：K42内得到单元格I2的值在对应时刻的数值解，图8.20是n=100，t=0.1时刻的数值解。
图8.20 n=100，t=0.1时刻的数值解
第八步，计算解析解。将单元格区域A46：L46合并为一个单元格并输入“=＂n=＂&2*H2&＂，t=＂&I2&＂时刻的解析解＂”。同第三步一样，在单元格区域A47：A57中输入自变量y的值；在单元格区域B58：L58中输入自变量x的值；在单元格B47中输入“=（SIN（PI（）*$A47）+SIN（PI（）*B$58））*EXP（-$A$2*PI（）*PI（）*$I$2）”；用自动填充功能拖动填充单元格区域B47：L57，如图8.21所示。
图8.21 解析解
第九步，绘制图像。选中单元格区域B33：L43，依次单击“插入”、图表下拉按钮、“更多散点图（
）...”，选择“插入图表→曲面图→曲面图”，单击“确定”按钮，就得到数值解的云图，如图8.22（a）所示；选中单元格区域B47：L57，依次单击“插入”、图表下拉按钮、“更多散点图
）...”，选择“曲面图→曲面图”，单击“确定”按钮，就得到解析解的云图，如图8.22（b）所示，为了更好地观察，在一个直角坐标系内绘制直线y=0.5上的数值解和解析解的比较图，如图8.23所示。
图8.22 n=100，t=0.1时刻的云图
图8.23 y=0.5上的数值解和解析解的比较图
