瞬态热传导方程-时域有限元
简述
热学中，稳态方程仅表示热传导的最终状态，而实际问题中，分析热传导的过程也是非常重要的，这就涉及瞬态热传导方程，也就是控制方程中加入时间的变量。
本文以二维矩形平板为例，实现时域有限元方法的瞬态热传导。
1.控制方程
考虑一个尺寸为 Lx×Ly的矩形平板，初始温度 T0。左侧边界施加恒定温度 T h，右侧边界对流换热（环境温度 Tinf，换热系数 h），上、下边界绝热。材料参数：导热系数k、密度ρ、比热容c。
二维瞬态热传导方程：
边界条件：
2.有限元公式推导
推导过程与传统有限元大同小异，首先乘以试探函数，并在求解域积分，然后通过分部积分，整理得到：
矩形模型左侧为第一类边界条件，上下边界为绝热边界不处理，右侧边界项代入对流条件：
使用矩形单元进行网格划分，形函数表示为N(x,y)，对温度离散，得到，可以发现温度T离散后保留的时间变量，仅对空间进行离散：
带入到有限元弱形式中，得到：
进一步化简：
空间上的离散过程与传统有限元一直，然后才是对时间上离散：
取theta=0.5，得到：
时间离散满足稳定性条件，即时间步长dt迭代过程中，时间步长导致的系统误差是收敛的，稳定性条件为：
以第一次迭代为说明，其初始条件的加载：
已知T0状态为初始状态，即整个区域为恒定温度T0,只在左端边界突然给定Th的值。此时，由此可以求解得到右端项，进而求解得到T1时刻的温度场。然后以T1为右端项输入值，计算得到T2，以此类推...直到整个区域温度再次达到稳态状态。
3.基函数推导
对于矩形单元而言，这里使用标准单元，将实际单元通过雅可比矩阵映射到标准单元的方法。
已知标准单元的基函数：
标准参考单元与实际单元之间存在等参变换关系：
对上述x,y在标准坐标系下进行偏导数处理，可以得到雅可比矩阵：
继续推导，可以得到形函数的导数在实际坐标系与参考坐标系的关系：
因此，有限元积分公式可以进一步写成：
如此，将实际单元的积分转化到参考单元的积分，进而可以使用高斯积分，求解每个单元的单元系数矩阵。
系数矩阵组装与经典有限元组装流程已知，这里不再过多介绍，具体可参考之前文章。
4.结果展示
物理模型参数如下：
% 参数定义Lx = 1.0; Ly = 1.0; % 区域尺寸nx = 20; ny = 20; % 单元数量k = 1000; % 导热系数rho = 7800; % 材料密度c = 450; % 比热容h = 100; % 对流系数T_inf = 25; % 常温工况Th = 100; % 左侧加热温度T0 = 20; % 初始温度
时间迭代步长2秒，求解总时长510秒，绘制期间时间迭代过程中部分时间的温度变换规律如下：
下面是时间迭代过程中，左端某个点与右端某个点的温度随时间的曲线如图：
最后
相比于传统的没有时间项的有限元求解，时域的有限元求解需要在每次时间迭代过程中多次求解方程组，其右端项是通过上一次迭代解更新得到的新的右端项。
本质上，在方程组装中，空间域还是基于泊松方程，时间域对时间进行差分离散，从而得到与时间迭代相关的系数矩阵方程。
