偏微分一维热传导问题

偏微分大作业偏微分大作业 一维热传导方程问题 运用隐式格式求解数值解 目录目录 问题描述 3 1 解析解 分离变量法 3 2 数值解 隐式格式 5 3 证明隐式格式的相容性与稳定性 5 4 数值解 分析与 Matlab 实现 6 5 数值解与解析解的比较 9 6 随时间变化的细杆上的温度分布情况 11 7 稳定后细杆上的温度分布情况 12 参考文献 13 附录 14 有限长杆的一维热传导问题 问题描述问题描述 一根单位长度的细杆放入 100 的沸水中 当细杆的温度达到 100 时取出 假设细杆四周绝热 在时间 t 0 时 细杆两端浸入 0 的冰水中 一维热传导 方程 现在令 从而可知本题 现 2 0 txx ua u 2 1a 0 txx uu 在要求细杆温度分布 u x t 1 解析解解析解 分离变量法分离变量法 热传导偏微分方程 0 txx uu 1 0 t 1 t 0uu 0 u xx 其中 001xx 或 x 100 0 1 x 首先令 2 u x tX x T t 将 2 式带入 1 式得 T 0X xtT t X x 于是可得 T tX x T tX x 可以得到两个微分方程 T 0tT t 0X xX x 先求解空间项 当时 0 xx X xAeBe 4 由于 0 t 1 t 0 uut 可知 由于解的收敛性 0B 0 1 00XXAAeA 则此时是平庸解 当时 0 X xABx 0 1 00 0XXAABAB 则此时是平庸解 当时 其中 0 cossinX xAkxBkx k 0 00XAA 1 sin0 1 2 3 XBkknn 所以 sin n X xBn x 1 2 3 n 因为 22 n 所以 22 nt n T tC e 1 2 3 n 则 22 1 sin nt n n u x tD en x 初始条件 0 u xx 1 0 sin n n u xDn xx 1 0 2 sin n Dxn x dx 1 2100sin 1 200 cos 1 cos n x dx nn n 200 0lim 1 cos n Dn n 当时 最终 5 22 1 200 1 1 sin nnt n u x ten x n 1 2 3 n 2 数值解数值解 隐式格式隐式格式 目前 研究热传导问题特别是非稳态热传导问题十分重要 这里使用隐式 格式 1 利用 关于 t 进行向前差商 关于 x 进行二阶中心 u x t 1kk jj UU t A 差 111 11 2 2 kkk jjj UUU x A 代入偏微分方程可以得到隐式差分格式 1111 11 2 2 kkkkk jjjjj UUUUU tx AA 1 3 证明隐式格式的相容性与稳定性证明隐式格式的相容性与稳定性 1 相容性 相容性 2 122 2 22 12233 1 22 22 12233 1 22 Taylor 1 2 11 22 11 22 kk jj kk jj kk jj UU UUttt tt UUUU UUttxxtx ttxx UUUU UUttxxtx ttxx AAA AAAAAA AAAAAA 根据展开 代入隐式格式得 2 22 22 22 1 2 UUU tttx ttx AAAA 将 2 与原微分方程相减 得到截断误差 6 2 2 1 0 2 U tx t AA 所以此隐式格式与原微分方程相容 2 稳定性 稳定性 令网格比为 则可以将 1 式改写得到 2 r t x A A 111 11 12 kkkk jjjj rUr UrUU 3 首先令 11 1 1 11 11 1 1 kkIj j kkIj j kkIj j kkIj j UUe UUe UU e UUe 4 将 4 代入 3 式 根据欧拉公式化简得 1 1 22 cos kk UrrU 5 故得放大因子是 1 1 1 1 2 1 cos k k U G Ur 所以根据 Fourier 方法 隐式格式恒稳定 4 数值解数值解 分析与分析与 Matlab 实现实现 1 边值与初值离散化边值与初值离散化 将边值与初值离散化 与式 3 联立得差分线性方程组 111 11 12 kkkk jjjj rUr UrUU 0 1 2 M 1 j 0 1 2 N 1 k 7 0 jj Ux 0 1 2 M j 0 0 k U 0 1 2 N k 0 k L U 0 1 2 N k 再将方程组改写成的形式 AUB 11 110 1 22 1 33 1 22 11 11 1 M 1 1 2 12 12 12 12 kkk kk kk kk MM kkk MMMM rrUUrU rrrUU rrUU rrrUU rrUUrU 本题的边界条件均为零 所以可以将上式改写 1 11 1 22 1 33 1 22 1 11 1 M 1 1 2 12 12 12 12 kk kk kk kk MM kk MMM rrUU rrrUU rrUU rrrUU rrUU 2 Matlab 的实现的实现 杆长 1 米 时间 2 秒 设计空间步长 h 0 1 和时间步长 t 0 01 网格比是 2 t r h 从而得到划分的空间网格点数是 M1 1 时间网格点数是 M2 1 先设初 始的温度矩阵 U M2 1 M1 1 再将边界条件和初始条件编写到表示温度 分布的矩阵中 具体代码可见最后附录 编写矩阵 A 核心代码 对角线 A i i 1 2r 对角线的右方和下方 A i i 1 r A i 1 i r 8 下面就要运用进行迭代 1 A U kjU k j 当 k 1 时 A U 2 j U 1 j 当 k 2 时 A U 3 j U 2 j 当 k 3 时 A U 4 j U 3 j 以此迭代下去直到 k M2 就可以得到整个温度随时间和空间的分布矩阵 U 数值解画图 如图 1 a 和图 1 b 所示 图 1 a 数值解的温度分布图 现在将着色平稳过渡 9 图 1 b 着色平稳过渡的数值解的温度分布图 5 数值解与解析解的比较数值解与解析解的比较 首先 我们需要将解析解离散化 解析解中有一项 当 n 越来越 22 nt e 大时 会快速趋于 0 故我们可以取 n 8000 现在来证明可行性 在 matlab 里的工作空间运算 将解析解的温度分布画出来 数值解画图 如图 2 所示 2 10 图 2 解析解的温度分布图 将数值解与解析解相减 得到误差图 如图 3 a 和图 3 b 我们从图 3 a 上可以看出空间上的误差空间上的误差 在边界处误差比较大 11 图 3 a 数值解与解析解空间误差 我们从图 3 a 上可以看出时间的误差时间的误差 在时间的最开始 处误差最大 然 后又有一个小的波动 最后就误差渐渐变小 最后趋于 0 图 3 b 数值解与解析解时间误差 6 随时间变化的细杆上的温度分布情况随时间变化的细杆上的温度分布情况 从数值解的温度分布三维图 如图 4 a 和图 4 b 可以看出随着时间的增加 细杆温度下降最后趋于 0 从物理角度来说 细杆的温度会不断地向两端扩散 热量会慢慢散失 最 终随着时间的增加 细杆的温度会趋于 0 12 图 4 a 细杆温度随时间的变化图 现取细杆中心处一点 观看它随时间的温度变化情况 图 4 b 细杆中央 x 0 5 温度随时间的变化图 13 7 稳定后细杆上的温度分布情况稳定后细杆上的温度分布情况 从图像上可以看出 最后稳定的情况下 细杆的温度是 0 参考文献参考文献 1 冯立伟 热传导方程几种差分格式的 MATLAB 数值解法的比较 J 沈阳化 工大学 辽宁沈阳 2011 6 2 一维热传导方程数值解法及 Matlab 实现 EB OL 2014 11 20 14 附录附录 代码 此程序用于解决一维热传导方程 ut a 2uxx 0 边界条件 u 0 t u L t 0 初始条件 u x 0 100 x 0 和 L u 0 0 0 u L 0 0 其中 a 2 1 L 1 clc clear all 区域及划分网格 L 1 单位长度的细杆 T 2 时间 h 0 1 空间的划分 t 0 01 时间的划分 r t h h 网格比 15 设计步长 M1 L h M2 T t 构造边界条件 构造的矩阵 U 时间 空间 U zeros M2 1 M1 1 编程包含边值 如 U k 1 u 0 t for k 1 M2 1 时间划分了 M2 份 有 M2 1 个节点 U k 1 0 两个边界处温度恒为零 U k M1 1 0 end 构造初始条件 for j 2 M1 位置划分了 M1 份 有 M1 1 个节点 U 1 j 100 end U 1 1 0 U 1 M1 1 0 差分格式的矩阵形式 A U k 1 j U k j 构造矩阵 A A zeros M1 1 for i 1 M1 1 A i i 1 2 r end for i 1 M1 2 A i i 1 r A i 1 i r end 构造 AU B 中的 B 本题边值的特殊 矩阵 B 大大简化了 B zeros M1 1 1 for k 1 M2 j 2 M1 B j 1 1 U k j x A B for j 2 M1 U k 1 j x j 1 k 1 时刻的不同位置的温度分布 end end 作图 x 0 h 1 y 0 t 2 xx yy meshgrid x y 16 figure 1 surf xx yy U shading flat title 一维热传导方程 数值解 温度分布图 xlabel 位置 x ylabel 时间 t zlabel 温度 T figure 2 s 0 for i 1 8000 s s 200 1 1 i i pi sin i pi xx exp i 2 pi 2 yy end surf xx yy s title 一维热传导方程 解析解 温度分布图 xlabel 位置 x ylabel 时间 t zlabel 温度 T figure 3 x 0 h 1 y 0 t 2 xx yy meshgrid x y dd U s surf xx