有限差分基础白.ppt

有限差分法主要步骤： 利用网格将区域离散处理； 构造差分格式 用差分、差商来代替微分、微商，将微分方程离散化为差分方程，并将定解条件离散化； 解线性代数方程组 建立差分格式后，微分方程的求解就可归结为求解一个线性代数方程组，通过解线性代数方程组，得到的是数值解。,本章主要以传热问题的有限差分计算为基础，介绍有限差分法求解的基本原理和过程。,第七章 有限差分基础,有限差分法：偏微分方程的一种数值解法。,差分法的基础是用差商代替微商。 若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, ),向前差公式（导数在点xi计算，而差商取fi及向前一点fi+1）,向后差公式（导数在点xi计算，而差商取fi及向后一点fi-1）,中心差公式（两侧差分平均值）,函数f(x)在x=xi处的二阶导数,函数f(x)的一阶导数,(xi-1,xi)和(xi,xi+1)两区间的一阶导数差除以x得到,一般地说，当差分公式的截断误差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内热源，为单位体积内热量产生的速率。,稳态基本方程为,P,I,O,E,W,N,S,上式可简化为(三维),利用式,可得近似式,一维导热的公式为,二维导热区域的公式为,二维稳态问题,差分方程为,若以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 Ta,Bi= hx/ 毕欧数 ； h表面放热系数，导热系数，Ta是环境温度；Ti,j边界节点温度。,内部导热；边界换热、对流或定温,2. 边界节点差分方程,热流边界条件,（沿y方向流入体内时）,用T对y的向前差商代替微商，则,或写成,绝热边界,可以写成,(Ti,j+1Ti,j)/y = q,Ti,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）,tn+1时刻(i, j)节点的温度Ti,jn+1 ，可以根据自身及其相邻节点在tn时刻的温度来计算，而tn时刻的温度是已知的。因此，结合初始条件和边界条件，就可以计算区域内各节点随时间t增长的温度值Ti,jn。,显式格式的优点 每个节点方程均可独立求解，整个计算过程十分方便。 缺点 若F0值取的不当，计算得到的解可能不稳定。因此，对时间步长的选取及网格的划分等要求比较严格。,若不考虑换热边界条件的影响，为保证稳定，必须要求,或写成,对于一维热流公式,等号右端用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 / F0,T2n = 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 / F0,T0n+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）误差 有限差分方程的近似解和偏微分方程精确解之间的差值。,指计算误差随步数的增加是否会积累到超出所允许的范围；或者说，最后计算结果对初始条件和边界条件的数据误差及计算中的舍入误差是否敏感的问题。,（2）差分格式的稳定性,保证差分格式的稳定性很重要。这是因为初始条件和边界条件不可避免地包含着误差，在数值计算中也有舍入误差。如果这些误差在计算过程中不断被放大，导致求解不稳定，那么计算结果就失去了意义。,差分格式的稳定性,显式有限差分方程稳定性条件,在计算起始阶段，容易产生振荡，出现不稳定现象。 缩短时间步长可以提高精度。但步长如果太小，精度虽然得到保证，但计算工作量加大。 比较好的方法是采用变步长计算。在开始阶段选小步长，经短时间后，逐步加大步长。既能保证精度又可节约计算时间。,（3）初始精度和步长选择,7.6 有限差分方程的应用,（1）二维稳态问题求解,考察图示平板问题。假设有一块长Lx，宽Ly，导热系数为的平板，且Lx=Ly。假设平板内热源强度为H=0。二维稳态热流问题的基本偏微分方程为：,换热边界,热流边界,绝热边界,固定温度边界,y=Ly, 0xLx,x=0, 0yLy,y=0, 0xLx,x=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=0,Ti,j = (Ti+1,j+BiTa)/(1+Bi),Bi =x /,对流换热边界,给定热流边界条件,绝热边界条件,给定温度边界条件,差分格式的建立,为了熟悉差分方程的建立和表示方法，本例将内节点（i, j）标为5。采用图示网格，并假定x=y=1。 针对每个节点列出差分方程，形成线性方程组：,Ti,jTi,j+1 = q/,Ti,j = Ti-1,j,Ti,j = Ta,图7-21 求解区域网格,Bi = /,图7-21 求解区域网格,把线性代数方程组写成矩阵形式,采用高斯消去法或迭代法求解该线性方程组。 上例中，网格划分很粗，方程数仅有9个。实际上差分计算时，网格划分得很细，节点数很多。因而在有限个节点上求得的温度值和连续的温度分布就相当接近。,（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）点，其有限差分方程的显式格式为,边界条件,根据问题的边界条件，节点的温度已知，即,根据显式表达式的稳定性判据，即 0F01/4,取F0=0.25时，根据傅立叶数公式，时间步长,T1,j=1000 j=1,2,3,4,T5,j=860 j=1,2,3,4,Ti,4=1000 i=2,3,4,Ti,1=300 i=2,3,4,求解,t =c2F0/=22s,故时间步长t不能超过22s。,7.7 有限差分方程的计算机解法,含有n个未知数和n个方程的方程组，其一般形式为,将方程简化为,n为节点数，aij和bi都是常数（i=1,2, , n，j=1,2, n）,（7-91）,A矩阵，通常是方阵；T节点温度列向量；B右端列向量,差分格式通常用矩阵形式表示,AT=B,矩阵A是一个稀疏的带状矩阵，含有大量的0元素，而非0元素位于主对角线两侧。,方程组有唯一解的充分必要条件是系数矩阵A是非奇异的，即它的行列式|A|0。,迭代法 一般比直接法好，系数矩阵是稀疏矩阵时，迭代解法更快些，并且没有舍入误差的积累。,直接解法 作有限次运算就可得到解，但是有舍入误差的积累，这种误差积累可能淹没计算的解。,用计算机求解方程组的方法：直接法和迭代法。,1. 高斯消去法,高斯消去法是线性代数方程组直接接法之一，其解题步骤分两步：第一步是消元过程，第二步是回代过程。,假定把要求的n阶线性代数方程组写成如下形式,（1）消元过程,首先消去第二个方程以后的n-1个方程中的未知量T1。,方法是：第二个方程减去第一个方程乘以a21(1)/a11(1)，第三个方程减去第一个方程乘以a31(1)/a11(1) ，直到最后一个方程。这样方程组变为等价方程组,方程中的系数aij(2)和bi(2)的计算公式为,其次，从等价方程的第3个方程开始，消去后面n-2个方程中的未知量T2。即第三个方程减去第二个方程乘以a32(2)/a22(2) ，第四个方程减去第二个方程乘以a42(2)/a22(2) ，如此下去，即可得到,方程中的系数计算公式为,上述步骤重复n-1次后，可得到如下等价三角形方程组,表示成矩阵形式为,A(n)T=B(n),式中，A(n)为上三角形式的矩阵，B(n)为n1阶向量,（1）,可以把消元过程的计算归结为对于k=1,2, ,n-1的递推公式,（2）,该式就是消元过程的基本公式。,（2）回代过程,经过消元以后的方程组为式（1）。按照由下往上的顺序，可依次解出Tn, Tn-1, , T1。,由式（1）最后一个方程可得,（3）,将解得的Tn代回倒数第2式，可解得Tn-1。依次类推，由第i式,可解得,（4）,式（3）和式（4）便是回代过程的基本公式。,2. 迭代法,迭代法的基本思想是，构造一个由T1，T2，Tn列向量序列，使其收敛于某个极限向量T1*, T2*, , Tn*，而且T1*, T2*, , Tn*就是方程组的精确解。 根据构造列向量的方法不同，有简单迭代法、高斯赛德尔迭代法和超松弛迭代法。,迭代法不需要存储系数矩阵中的零元素，所以占用的存储单元少，程序也比较简单。其缺点是需要进行多次迭代才能达到收敛指标要求，耗费机时较多。,（1）简单迭代法,迭代的最终目的是求解方程组(7-91)中的T1,T2,Tn。 若系数矩阵A对角线元素不为0，即aii0 (i=1,2,n)，则可将式(7-91)改写为：,其中任一方程均可写成,任意给定Ti (0)作为解的第零次近似，把它们代入上式右端，得,（7-99）,（7-98）,作为解的第一次近似。把第一次近似得到的解再代入（7-99）的右端，得到解的第二次近似 。,一般地，将已得到的解