有限差分法主要步骤：利用网格将区域离散处理；构造差分格式 用差分、差商来代替微分、微商，将微分方程离散化为差分方程，并将定解条件离散化；解线性代数方程组 建立差分格式后，微分方程的求解就可归结为求解一个线性代数方程组，通过解线性代数方程组，得到的是数值解。本章主要以传热问题的有限差分计算为基础，介绍有限差分法求解的基本原理和过程。第七章 有限差分基础有限差分法：偏微分方程的一种数值解法。差分法的基础是用差商代替微商。若y=f(x)是连续函数，则它的导数为f/ x差商，df/dx微商。在 x到达零以前， f/ x 只是df/dx的近似，两者的差值| f/ x - df/dx |表示差商代替微商的偏差。用差商代替微商，则微分方程就变成了差分方程。1. 差商与微商7.1 有限差分法基础2. 差分公式偏微分方程数值解法的基本原理是用几个相邻点的函数值和相邻点的间距来表示某点的导数。邻点间的距离可以相等，也可以各不相等 。考虑函数f(x)，将自变量x等间距离散化，取步长为x，令xi=ix ，fi=f (xi) (i=0,1, )则依据所取函数值的不同，可得到不同形式的差分公式。 现只讨论等间距即均匀网格中函数的导数。 x=xi+1-xi (i=0,1, )xxxixi+1xi-1xfi-1fifi+1f (x)f (x)O向前差公式（导数在点xi计算，而差商取fi及向前一点fi+1） 向后差公式（导数在点xi计算，而差商取fi及向后一点fi-1） 中心差公式（两侧差分平均值） 函数f(x)在x=xi处的二阶导数函数f(x)的一阶导数(xi-1,xi)和(xi,xi+1)两区间的一阶导数差除以x得到xxxixi+1xi-1xfi-1fifi+1f (x)f (x)O一般地说，当差分公式的截断误差E=O(x p)时，则称其具有p阶精度。 向前差公式 在x=xi展开得， E=O(x)；向后差公式 在x=xi展开得，E=O(x)；中心差公式 在x=xi展开得， E=O(x2)；二阶导数公式 在x=xi展开得， E=O(x2)。对差分公式按泰勒级数展开，可得各自的截断误差E。可见，后两个公式比前两个公式精度高一阶。截断误差注意：中心差公式虽然具有二阶截断误差，但这并不意味着它一定能得出比具有一阶截断误差的向前、向后差分公式更为精确的结果。7.2 有限差分的基本原理存在初值的一维热传导问题，可以用下式表示 在给定条件下，上述偏微分方程有唯一确定的解。(1) 定解区域的离散化用网格线将定解区域离散化为节点集，是将微分方程定解问题离散化为差分方程的基础。（7-2）一维热传导或扩散方程：（7-1）其中：称为导温系数或扩散系数图 7-1定解区域网格线节点: 网格线的交点空间步长：平行于t轴的网格线间距 时间步长：平行于x轴的网格线间距 网格线: 其中 节点()常简记为( ) 初值问题的解u是依赖连续变化的变量x和t的函数。采用有限差分法求解u在节点上的近似值。也就是说，把依赖连续变化x和t的问题归结为依赖离散变化i和j的问题。(2) 差分方程的建立（差分格式的构造）对于节点(i, j)，u的偏微商与差商之间有以下关系将上面两式代入式（7-1），并去掉O(x2+t)项，可得（7-6） 该式称为方程（7-2）的有限差分方程。（7-3）（7-4）改写成便于计算的形式: 称为网格比 其中 式（7-1）中的初始条件： 其中 差分格式：通常把定解问题中的微分方程的差分方程和定解条件的离散形式统称为定解问题的一个差分格式 显式差分格式 下一时刻节点的函数值可由当前时刻直接计算得到隐式差分格式 差分格式在t=(j+1) t时间层上包含多于一个节点的未知数传热分析用到的物理参数及其单位: 温度传导系数（导热系数）密度比热容导温系数（热扩散系数）时间热流密度：单位时间通过单位面积的热流量 内热源强度表面放热系数7.3 热传导问题 1热传导基本方程Tt时刻点(x,y,z)处的温度；为导热系数， =/c导温系数或热扩散率； 密度，c比热容，H内热源强度（单位体积的产热量）。热物性参数不随温度变化，且各向同性。 稳态时， ，有 若温度场内无内热源，即H=0, 该式即为拉普拉斯方程(Laplace)。初始条件 指某一时刻导热物体的温度分布。对于稳定导热问题，温度场不随时间变化，时间条件自然消失。温度随时间变化时，给出某一瞬时物体内部各点温度。t=0时物体内部的温度分布规律通常为T|t=0=T0(x,y,z)2导热问题的定解条件边界条件 即物体边界上的换热条件。常分为三类： 第一类：已知物体边界的温度，即Ts=T0(x,y,z,t) 第二类：已知物体边界上各点的热流密度，即 第三类：已知物体边界与周围介质的热交换，即 式中，n为边界外法线方向， 为外法向导数，h为表面放热系数，Ta为周围介质的温度。 当 时，即表示与外界无热交换，即绝热条件.实际问题往往是上述三类边界条件的组合。 7.4 稳态传热问题的有限差分方程对于多变量函数T=T( x, y)，涉及到求一阶和二阶偏导数的近似值。若把y看作常数，则函数T对于x的偏导数就是T对x的普通导数。同样，若把x看作常数，则函数T对于y的偏导数就是T对y的普通导数。因此，可以直接应用前面介绍的所有导数的概念和公式。当然，在应用前面的公式对x求偏导时，必须保持y = y0。首先只考察内部节点。图7-4所示是直角坐标系下一个三维导热区域中的网格点P及其六个相邻点，它们分别记为N，S，E，W，I，O。令网格间距x= y=z =。1. 内节点差分方程图7-4 在均匀网格的三维直角坐标中典型点P及其六个相邻点式中，H内热源，为单位体积内热量产生的速率。稳态基本方程为PIOEWNS上式可简化为(三维) 利用式可得近似式一维导热的公式为 二维导热区域的公式为二维稳态问题差分方程为若以Ti,j表示 (i, j)点温度，则同样，一维稳态问题的差分方程为(边界条件的差分形式)采用差商代替微商的办法把定解问题中的各种边界条件表示成差分形式 。给定温度边界 换热边界条件 用T对x的向前差商代替T对x的一阶微商，则 或写成Ti,j = Ts(Ti+1,jTi,j)/x=h(Ti,j Ta)(Bi+1)Ti,j Ti+1,j =Bi TaBi= hx/ 毕欧数 ； h表面放热系数，导热系数，Ta是环境温度；Ti,j边界节点温度。 内部导热；边界换热、对流或定温2. 边界节点差分方程热流边界条件 （沿y方向流入体内时） 用T对y的向前差商代替微商，则或写成绝热边界 可以写成(Ti,j+1Ti,j)/y = qTi,j Ti,j+1= qy/Ti,j Ti-1,j=0 或 Ti,j Ti, j-1=0 或每一个边界节点只应属于一种边界条件。在两种边界条件交接的节点，可人为规定属于哪一种的边界条件。对应边界的差分方程均采用一阶向前差商代替一阶微商得到，其截断误差为O(x)或O(y)量级，比内节点差分方程的截断误差低一个数量级。为了提高整个差分格式的计算精度，可对上述边界条件作进一步处理，如用中心差商代替微商等。注意：7.5 非稳态的有限差分方程非稳态的有限差分方程非稳态或瞬变传热问题的特征是热流和温度场随时间而变，因此离散化包含两个方面：空间域离散 几何区域离散化，确定内节点、边界节点时间域离散 热过程经历的时间区域离散化。在构造非稳态传热的差分方程时，必须特别注意它的稳定性，因为用不稳定的差分方程进行求解是没有意义的。此外，在边界条件差分形式的处理上，也有新的特点需要考虑。几何区域离散化。假定区域离散化后，距离步长x=xi+1-xi, y=yj+1-yj，且x=y。显然，xi=ix；yj=jy，i,j=0,1,2,。时间域离散化。用n（n=0,1,2,）将时间区域t0离散化，两个时刻的间隔（时间步长）t=tn+1-tn，tn=nt。1. 二维非稳态热传导方程(1) 离散化规定： (i, j )(xi, yj) ，ntn Tni,j n时刻节点（i,j）处的温度T(xi, yj, tn)。显示差分格式 将导热微分方程应用于时刻n的节点（i, j），可写成（n0） （2） 式(2)等号两侧的偏微分用差商来近似(2) 差分格式采用不同的差分公式，可建立不同形式的差分方程。温度对时间一阶向前差商来近似二阶偏微分用中心差商来近似将三式代入式（2），得相应的差分方程为 该式即为微分方程的差分方程，截断误差为O(t+x2+y2)。令x=y=， 代入式（6）并整理，得F0傅立叶数， （6） t(x )2t(y )2t2F0 = = =tn+1时刻(i, j)节点的温度Ti,jn+1 ，可以根据自身及其相邻节点在tn时刻的温度来计算，而tn时刻的温度是已知的。因此，结合初始条件和边界条件，就可以计算区域内各节点随时间t增长的温度值Ti,jn。显式格式的优点每个节点方程均可独立求解，整个计算过程十分方便。缺点 若F0值取的不当，计算得到的解可能不稳定。因此，对时间步长的选取及网格的划分等要求比较严格。 若不考虑换热边界条件的影响，为保证稳定，必须要求t(x )2F0x = 1/4t(y )2F0y = 1/4或写成对于一维热流公式 等号右端用n+1时刻的一阶向后差商来近似，而等号左端温度对距离的二阶偏微商则对应tn+1时刻，故相应的差分方程为: 差分方程的截断误差也是O(t +x2 +y2) 。完全隐式格式将导热微分方程应用于时刻n+1的节点（i, j），可写成上式包含邻点tn+1时刻的温度值。因此，从tn时刻的值不能简单地计算出（i, j）点tn+1时刻的温度，必须在每一个时间步长内求解一组联立方程才能求得Ti,jn+1（这组方程的数目等于待求温度的节点总数）。故称这种差分格式为隐式差分格式。隐式差分格式多种多样，式（12）的差分形式称为完全隐式差分格式。其优点是它不受边界条件、步长的影响，是无条件稳定的格式。x=y=时，该式可简写为（12） 对于一维热流公式为将对应节点（i, j）的微分方程写成如下形式 式中，为加权系数，取值范围为01。加权差分格式加权差分格式为简化为当=0时，显式格式。当=1时，完全隐式格式。当=1/2时，C-N格式。当=2/3时，加辽金格式。从0到1变化时，可得到不同的差分格式。对于给定的t和，随着的增长，计算精度下降，稳定性却越能得到保证。隐式格式对于二、三类边界条件，在边界外设立虚节点，使边界节点变换为内节点：用中心差商近似一阶微商；边界节点取内节点差分方程。2. 非稳态问题的边界条件在开始进行计算的一瞬间（t=0），边界温度突然由T0变为Tw不大合理。因此，实际计算时，应作适当处理。（1）给定温度Tw第一步计算(t=0时) 边界节点温度为Tw/2；完成第一步计算之后 固定温度边界节点保持Tw温度为提高整个差分格式的计算精度，常用中心差商来代替边界上的一阶微商。为此，在边界外设虚假节点。（2）给定换热边界条件在具有边长为正方形网格的二维矩形区域中，有两类节点：一类是边上（如节点1），另一类在角上（如节点0）。对于边节点1 在边界外与节点2对称的位置设以虚假节点2。这样，换热边界条件中的偏微商可用中心差商近似。 节点1变为内节点，其显式差分方程为将T2n表达式代入上式，消去T2n后，可得上式即为边界节点1的差分方程，其截断误差为O(2)，与内节点差分方程的相一致。 或写成用中心差商代替边界条件中的偏微商，得(T2nT2n) /(2)=h(T1nTa)T2n + T2n + T4n + T0n(41/F0)T1n= T1n+1 / F0T2n = T2n2Bi(T1nTa)T1n+1 = F0T0n + 2T2n + T4n + 2BiTa + (1/F042Bi)T1n 同样可得换热边界条件的隐式差分格式：T1n = F0-T0n+1 - 2T2n+1 - T4n+1 - 2BiTa + (1/F0+4+2Bi)T1n+1 对于角节点0在1和3的对称位置设虚节点1和3，在x，y方向分别应用中心差商格式(T3nT3n) /(2)=h(T0nTa)T3n = T3n2Bi(T0nTa)(T1nT1n) /(2)=h(T0nTa)T1n = T1n2Bi(T0nTa)或内节点0显式差分格式为 同样可以隐式形式表示T3n + T3n + T1n + T1n(41/F0)T0n= T0n+1 / F0T0n+1 = 2F0 T1n + T3n + 2BiTa + (1/(2F0)22Bi)T0n T0n = 2F0 -T1n+1 - T3n+1 - 2BiTa + (1/(2F0)+2+2Bi)T0n+1 设在x=0处为给定热流q的边界条件，且保持不变，则边界条件可写成用中心差商代替微商（方法下同）代入内节点显式差分格式 （3）给定热流边界得热流边界在边界外设立虚节点，使边界节点变换为内节点隐式格式代入下式 得显式格式隐式格式（4）绝热边界3. 求解的精确性和稳定性 （1）误差 有限差分方程的近似解和偏微分方程精确解之间的差值。 误差类型说 明对精度的影响截断误差用有限差分代替导数引起。取决于初始给定的温度分布、边界条件、差分格式和傅立叶数Fo影响显著数值误差即舍入误差。计算中对有效数字的限制引起的很少或根本不影响指计算误差随步数的增加是否会积累到超出所允许的范围；或者说，最后计算结果对初始条件和边界条件的数据误差及计算中的舍入误差是否敏感的问题。 （2）差分格式的稳定性保证差分格式的稳定性很重要。这是因为初始条件和边界条件不可避免地包含着误差，在数值计算中也有舍入误差。如果这些误差在计算过程中不断被放大，导致求解不稳定，那么计算结果就失去了意义。 差分格式稳定性和精度应 用显式差分格式有条件稳定。t2/(k)（k=2,4,6; 1D,2D,3D）完全隐式格式无条件稳定。精度较低，但其稳定性最好C-N格式无条件稳定。精度比加辽金格式高，但初始精度差F0较小或边界换热较缓慢的场合加辽金格式无条件稳定。初始精度比C-N好F0较大或边界换热较剧烈的场合加权格式1/21时，无条件稳定差分格式的稳定性显式有限差分方程稳定性条件 直角坐标内部节点边界节点 (表面放热)备 注一维F01/2F0 1/2(1+Bi)傅立叶数F0=t/(2)毕欧数 Bi=h/-导温系数h-换热系数-导热系数-距离步长二维F01/4F0 1/2(2+Bi)F0 1/4(1+Bi) （角节点）三维F01/6F0 1/2(3+Bi)在计算起始阶段，容易产生振荡，出现不稳定现象。缩短时间步长可以提高精度。但步长如果太小，精度虽然得到保证，但计算工作量加大。比较好的方法是采用变步长计算。在开始阶段选小步长，经短时间后，逐步加大步长。既能保证精度又可节约计算时间。（3）初始精度和步长选择7.6 有限差分方程的应用（1）二维稳态问题求解考察图示平板问题。假设有一块长Lx，宽Ly，导热系数为的平板，且Lx=Ly。假设平板内热源强度为H=0。二维稳态热流问题的基本偏微分方程为： 换热边界 热流边界 绝热边界 固定温度边界 y=Ly, 0xLxx=0, 0yLyy=0, 0xLxx=Lx, 0yLy边界条件环境温度Ta求解平板内温度分布T（x，y）。 将求解区域进行网格划分。x, y方向的步长分别为x, y，节点坐标为(i, j)，i, j为整数。节点温度为Ti,j。 离散化内节点的差分方程 对于内节点（i，j）的差分方程为 边界节点的差分形式Ti-1,j+Ti+1,j+Ti,j-1+Ti,j+14Ti,j=0Ti,j = (Ti+1,j+BiTa)/(1+Bi)Bi =x / 对流换热边界 给定热流边界条件 绝热边界条件 给定温度边界条件 差分格式的建立 为了熟悉差分方程的建立和表示方法，本例将内节点（i, j）标为5。采用图示网格，并假定x=y=1。针对每个节点列出差分方程，形成线性方程组： Ti,jTi,j+1 = q/Ti,j = Ti-1,jTi,j = Ta图7-21 求解区域网格节点编号差分方程1T1-T2=q/2(Bi+1)T2-T5=BiTa3T3=Ta4T4-T5=q/5-T2-T4+4T5-T6-T8=06T6=Ta7T7-T8=q/8T8-T5=09T9=TaBi = /图7-21 求解区域网格Ti-1,j+Ti+1,j+Ti,j-1+Ti,j+14Ti,j=0Ti,j = (Ti+1,j+BiTa)/(1+Bi)Ti,jTi,j+1 = q/Ti,j = Ti-1,jTi,j = Ta内节点换热热流绝热固定把线性代数方程组写成矩阵形式采用高斯消去法或迭代法求解该线性方程组。上例中，网格划分很粗，方程数仅有9个。实际上差分计算时，网格划分得很细，节点数很多。因而在有限个节点上求得的温度值和连续的温度分布就相当接近。=1 -1 B+1 -1 1 1 -1 -1 -1 4 -1 -1 1 1 -1 -1 1 1 T1T2T3T4T5T6T7T8T9q/BTaTaq/0Taq/0Ta（2）二维非稳态问题求解 考察大小为0.120.15m的矩形钢板。假定任何边界的温度已知，钢板初始温度为20，没有内热源和边界热流的作用。钢的导热系数= 41J/(ms)，比热c =504J/kg，密度=8103kg/m3。求5min后钢板的温度分布。图7-23 矩形钢板Ta=1000(x=0,0y0.12)Tb=1000(y=0.12,0x0.15)Tc=300 (y=0，0x0.15) Td=860 (x=0.15,0y0.12)已知边界温度将的求解区域划分成正方形网格。x=y =0.03m ，节点(i, j)的温度用Ti,j表示。 离散化内部节点的差分方程内部节点的显式差分方程 傅立叶数 例如对于（2，2）点，其有限差分方程的显式格式为边界条件 根据问题的边界条件，节点的温度已知，即 根据显式表达式的稳定性判据，即 00)，即满足： （2）高斯赛德尔迭代法简单迭代法计算程序简单，但计算每一个Ti(k+1)都要用全部的Ti(k)值。因此在计算机上必须有两套工作单元来存放全部未知的旧值和新值。为了节省工作单元并加快收敛速度，对上述计算作如下处理： 方程组（7-98）的第一式仍改写成方程组（7-98）中第二式中换成第一次得到的，即按此规律，可得Ti(k+1)的一般关系式 （7-100） 可见，用高斯赛德尔迭代法求解时，只需用一套工作单元存放近似值Ti(k)或Ti(k+1)，这样节省了工作单元，同时在迭代过程中，有一半用迭代的新值，加快收敛速度。因此，该法是一种常用的方法 （3）超松弛迭代法当未知数很多时，高斯赛德尔迭代法的收敛速度仍然较慢，作为改进，人们又提出了超松弛迭代法。 对于方程组（7-91），先假定第零次近似向量序列为T1(0),T2(0),Tn(0)。现在分两部来求第一次近似。第一步用高斯赛德尔迭代法 求出第一个近似解Zi(1)；第二步按下式来改善Zi(1)，并得第一次的近似解Ti(1) 是一个适当选择的参数，用它来改善 Zi(1)。由上述两式消去Ti(1) ，可以得出由Ti(0) 直接求出Ti(1)的公式。一般来说，由Ti(k)算得Ti(k+1)的公式为： 这就是超松弛迭代法的计算式，称为松弛因子。（7-103）当 =1时，式（7-103）简化为高斯赛德尔迭代法的计算公式。因此，超松弛法是更为一般的迭代公式。取不同值时，超松弛法收敛的快慢也不同。可以找出一个最佳的松弛因子*，使超松弛法收敛得比取别的值更快。 *值范围为1*2。对于01的情况称为超松弛。 能量平衡法非均匀网格非直角坐标系（柱状、球坐标）三角形网格界面、复合体下列内容自学（一般了解）