有限差分基础白

有限差分基础白第一页，共七十四页，2022年，8月28日差分法的基础是用差商代替微商。若y=f(x)是连续函数，则它的导数为△f/△

x—差商，df/dx—微商。在△

x到达零以前，△

f/△

x只是df/dx的近似，两者的差值|△

f/△

x-df/dx|表示差商代替微商的偏差。用差商代替微商，则微分方程就变成了差分方程。1.差商与微商1有限差分法基础第二页，共七十四页，2022年，8月28日2.差分公式偏微分方程数值解法的基本原理是用几个相邻点的函数值和相邻点的间距来表示某点的导数。邻点间的距离可以相等，也可以各不相等。考虑函数f(x)，将自变量x等间距离散化，取步长为△x，令xi=i△x，fi=f(xi)(i=0,1,…)则依据所取函数值的不同，可得到不同形式的差分公式。现只讨论等间距即均匀网格中函数的导数。△x=xi+1-xi(i=0,1,…)△x△xxixi+1xi-1xfi-1fifi+1f(x)f(x)O第三页，共七十四页，2022年，8月28日向前差公式（导数在点xi计算，而差商取fi及向前一点fi+1）向后差公式（导数在点xi计算，而差商取fi及向后一点fi-1）中心差公式（两侧差分平均值）函数f(x)在x=xi处的二阶导数函数f(x)的一阶导数(xi-1,xi)和(xi,xi+1)两区间的一阶导数差除以Δx得到△x△xxixi+1xi-1xfi-1fifi+1f(x)f(x)O第四页，共七十四页，2022年，8月28日一般地说，当差分公式的截断误差E=O(△xp)时，则称其具有p阶精度。向前差公式在x=xi展开得，E=O(△x)；向后差公式在x=xi展开得，E=O(△x)；中心差公式在x=xi展开得，E=O(△x2)；二阶导数公式在x=xi展开得，E=O(△x2)。对差分公式按泰勒级数展开，可得各自的截断误差E。可见，后两个公式比前两个公式精度高一阶。截断误差第五页，共七十四页，2022年，8月28日2有限差分的基本原理存在初值的一维热传导问题，可以用下式表示在给定条件下，上述偏微分方程有唯一确定的解。(1)定解区域的离散化用网格线将定解区域离散化为节点集，是将微分方程定解问题离散化为差分方程的基础。（2）一维热传导或扩散方程：（1）其中：称为导温系数或扩散系数第六页，共七十四页，2022年，8月28日图1定解区域网格线节点:网格线的交点空间步长：平行于t轴的网格线间距时间步长：平行于x轴的网格线间距网格线:其中节点()常简记为()初值问题的解u是依赖连续变化的变量x和t的函数。采用有限差分法求解u在节点上的近似值。也就是说，把依赖连续变化x和t的问题归结为依赖离散变化i和j的问题。第七页，共七十四页，2022年，8月28日(2)差分方程的建立（差分格式的构造）对于节点(i,j)，u的偏微商与差商之间有以下关系将上面两式代入式（1），并去掉O(Δx2+Δt)项，可得（5）该式称为方程（1）的有限差分方程。（3）（4）第八页，共七十四页，2022年，8月28日改写成便于计算的形式:称为网格比其中式（1）中的初始条件：其中差分格式：通常把定解问题中的微分方程的差分方程和定解条件的离散形式统称为定解问题的一个差分格式

显式差分格式下一时刻节点的函数值可由当前时刻直接计算得到隐式差分格式差分格式在t=(j+1)△t时间层上包含多于一个节点的未知数第九页，共七十四页，2022年，8月28日传热分析用到的物理参数及其单位:温度传导系数（导热系数）密度比热容导温系数（热扩散系数）时间热流密度：单位时间通过单位面积的热流量内热源强度表面放热系数第十页，共七十四页，2022年，8月28日3热传导问题1．热传导基本方程T—t时刻点(x,y,z)处的温度；λ—为导热系数，α=λ/ρc—导温系数或热扩散率；ρ—密度，c—比热容，H—内热源强度（单位体积的产热量）。热物性参数不随温度变化，且各向同性。

稳态时，，有若温度场内无内热源，即H=0,该式即为拉普拉斯方程(Laplace)。第十一页，共七十四页，2022年，8月28日初始条件指某一时刻导热物体的温度分布。对于稳定导热问题，温度场不随时间变化，时间条件自然消失。温度随时间变化时，给出某一瞬时物体内部各点温度。t=0时物体内部的温度分布规律通常为T|t=0=T0(x,y,z)2．导热问题的定解条件边界条件即物体边界上的换热条件。常分为三类：第一类：已知物体边界的温度，即Ts=T0(x,y,z,t)第十二页，共七十四页，2022年，8月28日第二类：已知物体边界上各点的热流密度，即第三类：已知物体边界与周围介质的热交换，即式中，n为边界外法线方向，为外法向导数，h为表面放热系数，Ta为周围介质的温度。当时，即表示与外界无热交换，即绝热条件.实际问题往往是上述三类边界条件的组合。第十三页，共七十四页，2022年，8月28日4稳态传热问题的有限差分方程对于多变量函数T=T(x,y)，涉及到求一阶和二阶偏导数的近似值。若把y看作常数，则函数T对于x的偏导数就是T对x的普通导数。同样，若把x看作常数，则函数T对于y的偏导数就是T对y的普通导数。因此，可以直接应用前面介绍的所有导数的概念和公式。当然，在应用前面的公式对x求偏导时，必须保持y=y0。第十四页，共七十四页，2022年，8月28日首先只考察内部节点。图4所示是直角坐标系下一个三维导热区域中的网格点P及其六个相邻点，它们分别记为N，S，E，W，I，O。令网格间距Δx=Δy=Δz=Δ。1.内节点差分方程图4在均匀网格的三维直角坐标中典型点P及其六个相邻点式中，H—内热源，为单位体积内热量产生的速率。稳态基本方程为第十五页，共七十四页，2022年，8月28日PIOEWNS第十六页，共七十四页，2022年，8月28日上式可简化为(三维)利用式可得近似式一维导热的公式为二维导热区域的公式为第十七页，共七十四页，2022年，8月28日二维稳态问题差分方程为若以Ti,j表示(i,j)点温度，则同样，一维稳态问题的差分方程为第十八页，共七十四页，2022年，8月28日(边界条件的差分形式)采用差商代替微商的办法把定解问题中的各种边界条件表示成差分形式。给定温度边界换热边界条件用T对x的向前差商代替T对x的一阶微商，则或写成Ti,j=Tsλ(Ti+1,j－Ti,j)/Δx=h(Ti,j–Ta)(Bi+1)Ti,j

－Ti+1,j

=Bi

TaBi=hΔx/λ—毕欧数；h—表面放热系数，λ—导热系数，Ta是环境温度；Ti,j—边界节点温度。内部导热；边界换热、对流或定温2.边界节点差分方程第十九页，共七十四页，2022年，8月28日热流边界条件（沿y方向流入体内时）用T对y的向前差商代替微商，则或写成绝热边界可以写成－λ(Ti,j+1－Ti,j)/Δy=qTi,j－Ti,j+1=qΔy/λTi,j－Ti-1,j=0或Ti,j－Ti,j-1=0或第二十页，共七十四页，2022年，8月28日每一个边界节点只应属于一种边界条件。在两种边界条件交接的节点，可人为规定属于哪一种的边界条件。对应边界的差分方程均采用一阶向前差商代替一阶微商得到，其截断误差为O(Δx)或O(Δy)量级，比内节点差分方程的截断误差低一个数量级。为了提高整个差分格式的计算精度，可对上述边界条件作进一步处理，如用中心差商代替微商等。注意：第二十一页，共七十四页，2022年，8月28日5非稳态的有限差分方程非稳态或瞬变传热问题的特征是热流和温度场随时间而变，因此离散化包含两个方面：空间域离散几何区域离散化，确定内节点、边界节点时间域离散热过程经历的时间区域离散化。在构造非稳态传热的差分方程时，必须特别注意它的稳定性，因为用不稳定的差分方程进行求解是没有意义的。此外，在边界条件差分形式的处理上，也有新的特点需要考虑。第二十二页，共七十四页，2022年，8月28日几何区域离散化。假定区域离散化后，距离步长Δx=xi+1-xi,Δy=yj+1-yj，且Δx=Δy。显然，xi=iΔx；yj=jΔy，i,j=0,1,2,…。时间域离散化。用n（n=0,1,2,…）将时间区域t≥0离散化，两个时刻的间隔（时间步长）Δt=tn+1-tn，tn=nΔt。1.二维非稳态热传导方程(1)离散化规定：(i,j)—(xi,yj)，n—tn

Tni,j

—n时刻节点（i,j）处的温度T(xi,yj,tn)。第二十三页，共七十四页，2022年，8月28日显示差分格式将导热微分方程应用于时刻n的节点（i,j），可写成（n>0）（2）式(2)等号两侧的偏微分用差商来近似(2)差分格式采用不同的差分公式，可建立不同形式的差分方程。温度对时间一阶向前差商来近似二阶偏微分用中心差商来近似第二十四页，共七十四页，2022年，8月28日将三式代入式（2），得相应的差分方程为该式即为微分方程的差分方程，截断误差为O(Δt+Δx2+Δy2)。令Δx=Δy=Δ，代入式（6）并整理，得F0—傅立叶数，（6）αΔt(Δx)2αΔt(Δy)2αΔtΔ2F0===第二十五页，共七十四页，2022年，8月28日tn+1时刻(i,j)节点的温度Ti,jn+1，可以根据自身及其相邻节点在tn时刻的温度来计算，而tn时刻的温度是已知的。因此，结合初始条件和边界条件，就可以计算区域内各节点随时间t增长的温度值Ti,jn。显式格式的优点每个节点方程均可独立求解，整个计算过程十分方便。缺点若F0值取的不当，计算得到的解可能不稳定。因此，对时间步长的选取及网格的划分等要求比较严格。若不考虑换热边界条件的影响，为保证稳定，必须要求αΔt(Δx)2F0x=≤1/4αΔt(Δy)2F0y=≤1/4第二十六页，共七十四页，2022年，8月28日或写成对于一维热流公式第二十七页，共七十四页，2022年，8月28日等号右端用n+1时刻的一阶向后差商来近似，而等号左端温度对距离的二阶偏微商则对应tn+1时刻，故相应的差分方程为:差分方程的截断误差也是O(Δt+Δx2+Δy2)。完全隐式格式将导热微分方程应用于时刻n+1的节点（i,j），可写成第二十八页，共七十四页，2022年，8月28日上式包含邻点tn+1时刻的温度值。因此，从tn时刻的值不能简单地计算出（i,j）点tn+1时刻的温度，必须在每一个时间步长内求解一组联立方程才能求得Ti,jn+1（这组方程的数目等于待求温度的节点总数）。故称这种差分格式为隐式差分格式。隐式差分格式多种多样，式（12）的差分形式称为完全隐式差分格式。其优点是它不受边界条件、步长的影响，是无条件稳定的格式。Δx=Δy=Δ时，该式可简写为（12）第二十九页，共七十四页，2022年，8月28日对于一维热流公式为第三十页，共七十四页，2022年，8月28日将对应节点（i,j）的微分方程写成如下形式式中，θ为加权系数，取值范围为0≤θ≤1。加权差分格式加权差分格式为第三十一页，共七十四页，2022年，8月28日简化为当θ=0时，显式格式。当θ=1时，完全隐式格式。当θ=1/2时，C-N格式。当θ=2/3时，加辽金格式。θ从0到1变化时，可得到不同的差分格式。对于给定的Δt和Δ，随着θ的增长，计算精度下降，稳定性却越能得到保证。隐式格式第三十二页，共七十四页，2022年，8月28日对于二、三类边界条件，在边界外设立虚节点，使边界节点变换为内节点：用中心差商近似一阶微商；边界节点取内节点差分方程。2.非稳态问题的边界条件在开始进行计算的一瞬间（t=0），边界温度突然由T0变为Tw不大合理。因此，实际计算时，应作适当处理。（1）给定温度Tw第一步计算(t=0时)边界节点温度为Tw/2；完成第一步计算之后固定温度边界节点保持Tw温度第三十三页，共七十四页，2022年，8月28日为提高整个差分格式的计算精度，常用中心差商来代替边界上的一阶微商。为此，在边界外设虚假节点。（2）给定换热边界条件在具有边长Δ为正方形网格的二维矩形区域中，有两类节点：一类是边上（如节点1），另一类在角上（如节点0）。对于边节点1在边界外与节点2对称的位置设以虚假节点2’。这样，换热边界条件中的偏微商可用中心差商近似。第三十四页，共七十四页，2022年，8月28日节点1变为内节点，其显式差分方程为将T2’n表达式代入上式，消去T2’n后，可得上式即为边界节点1的差分方程，其截断误差为O(Δ2)，与内节点差分方程的相一致。或写成用中心差商代替边界条件中的偏微商，得λ(T2n－T2’n)/(2Δ)=h(T1n－Ta)T2n+T2’n

+T4n+T0n－(4－1/F0)T1n=

T1n+1/F0T2’n

=T2n－2Bi(T1n－Ta)T1n+1

=F0[T0n+2T2n

+T4n+2BiTa+(1/F0－4－2Bi)T1n

]第三十五页，共七十四页，2022年，8月28日同样可得换热边界条件的隐式差分格式：T1n

=F0[-T0n+1-2T2n+1-T4n+1-2BiTa+(1/F0+4+2Bi)T1n+1

]对于角节点0在1和3的对称位置设虚节点1’和3’，在x，y方向分别应用中心差商格式λ(T3n－T3’n)/(2Δ)=h(T0n－Ta)T3’n

=T3n－2Bi(T0n－Ta)λ(T1n－T1’n)/(2Δ)=h(T0n－Ta)T1’n

=T1n－2Bi(T0n－Ta)或第三十六页，共七十四页，2022年，8月28日内节点0显式差分格式为同样可以隐式形式表示T3n+T3’n

+T1n+T1’n－(4－1/F0)T0n=

T0n+1/F0T0n+1

=

2F0[T1n+T3n

+

2BiTa+(1/(2F0)－2－2Bi)T0n

]T0n

=

2F0[-T1n+1-T3n+1-2BiTa+(1/(2F0)+2+2Bi)T0n+1

]第三十七页，共七十四页，2022年，8月28日设在x=0处为给定热流q的边界条件，且保持不变，则边界条件可写成用中心差商代替微商（方法下同）代入内节点显式差分格式（3）给定热流边界得热流边界在边界外设立虚节点，使边界节点变换为内节点第三十八页，共七十四页，2022年，8月28日隐式格式代入下式得显式格式隐式格式（4）绝热边界第三十九页，共七十四页，2022年，8月28日3.求解的精确性和稳定性（1）误差有限差分方程的近似解和偏微分方程精确解之间的差值。误差类型说明对精度的影响截断误差用有限差分代替导数引起。取决于初始给定的温度分布、边界条件、差分格式和傅立叶数Fo影响显著数值误差即舍入误差。计算中对有效数字的限制引起的很少或根本不影响第四十页，共七十四页，2022年，8月28日指计算误差随步数的增加是否会积累到超出所允许的范围；或者说，最后计算结果对初始条件和边界条件的数据误差及计算中的舍入误差是否敏感的问题。（2）差分格式的稳定性保证差分格式的稳定性很重要。这是因为初始条件和边界条件不可避免地包含着误差，在数值计算中也有舍入误差。如果这些误差在计算过程中不断被放大，导致求解不稳定，那么计算结果就失去了意义。第四十一页，共七十四页，2022年，8月28日显式有限差分方程稳定性条件直角坐标内部节点边界节点(表面放热)备注一维F0≤1/2F0≤1/[2(1+Bi)]傅立叶数F0=Δt/(αΔ2)毕欧数Bi=hΔ/λα-导温系数h-换热系数λ-导热系数Δ-距离步长二维F0≤1/4F0≤1/[2(2+Bi)]F0≤1/[4(1+Bi)]

（角节点）三维F0≤1/6F0≤1/[2(3+Bi)]第四十二页，共七十四页，2022年，8月28日在计算起始阶段，容易产生振荡，出现不稳定现象。缩短时间步长可以提高精度。但步长如果太小，精度虽然得到保证，但计算工作量加大。比较好的方法是采用变步长计算。在开始阶段选小步长，经短时间后，逐步加大步长。既能保证精度又可节约计算时间。（3）初始精度和步长选择第四十三页，共七十四页，2022年，8月28日6有限差分方程的应用例（1）二维稳态问题求解考察图示平板问题。假设有一块长Lx，宽Ly，导热系数为λ的平板，且Lx=Ly。假设平板内热源强度为H=0。二维稳态热流问题的基本偏微分方程为：换热边界热流边界绝热边界固定温度边界y=Ly,0≤x≤Lxx=0,0<y<Lyy=0,0≤x≤Lxx=Lx,0<y<Ly边界条件环境温度Ta求解平板内温度分布T（x，y）。第四十四页，共七十四页，2022年，8月28日将求解区域进行网格划分。x,y方向的步长分别为Δx,Δy，节点坐标为(i,j)，i,j为整数。节点温度为Ti,j。离散化内节点的差分方程对于内节点（i，j）的差分方程为边界节点的差分形式Ti-1,j+Ti+1,j+Ti,j-1+Ti,j+1－4Ti,j=0Ti,j

=(Ti+1,j+BiTa)/(1+Bi)Bi=αΔx/λ对流换热边界第四十五页，共七十四页，2022年，8月28日给定热流边界条件绝热边界条件给定温度边界条件差分格式的建立为了熟悉差分方程的建立和表示方法，本例将内节点（i,j）标为5。采用图示网格，并假定Δx=Δy=1。针对每个节点列出差分方程，形成线性方程组：Ti,j－Ti,j+1=qΔ/λTi,j

=Ti-1,jTi,j

=Ta图21求解区域网格第四十六页，共七十四页，2022年，8月28日节点编号差分方程1T1-T2=q/λ2(Bi+1)T2-T5=BiTa3T3=Ta4T4-T5=q/λ5-T2-T4+4T5-T6-T8=06T6=Ta7TT8=q/λ8T8-T5=09T9=TaBi=αΔ/λ图21求解区域网格Ti-1,j+Ti+1,j+Ti,j-1+Ti,j+1－4Ti,j=0Ti,j

=(Ti+1,j+BiTa)/(1+Bi)Ti,j－Ti,j+1=qΔ/λTi,j

=Ti-1,jTi,j

=Ta内节点换热热流绝热固定第四十七页，共七十四页，2022年，8月28日把线性代数方程组写成矩阵形式采用高斯消去法或迭代法求解该线性方程组。上例中，网格划分很粗，方程数仅有9个。实际上差分计算时，网格划分得很细，节点数很多。因而在有限个节点上求得的温度值和连续的温度分布就相当接近。=1-1

B+1-111-1-1-14-1-111-1-111T1T2T3T4T5T6T7T8T9q/λBTaTaq/λ0Taq/λ0Ta第四十八页，共七十四页，2022年，8月28日例（2）二维非稳态问题求解考察大小为0.12×0.15m的矩形钢板。假定任何边界的温度已知，钢板初始温度为20℃，没有内热源和边界热流的作用。钢的导热系数λ=41J/(m·s·℃)，比热c=504J/kg·℃，密度ρ=8×103kg/m3。求5min后钢板的温度分布。图23矩形钢板Ta=1000℃(x=0,0≤y≤0.12)Tb=1000℃(y=0.12,0<x<0.15)Tc=300℃(y=0，0<x<0.15)Td=860℃(x=0.15,0≤y<0.12)已知边界温度第四十九页，共七十四页，2022年，8月28日将的求解区域划分成正方形网格。Δx=Δy=0.03m，节点(i,j)的温度用Ti,j表示。离散化内部节点的差分方程内部节点的显式差分方程傅立叶数例如对于（2，2）点，其有限差分方程的显式格式为第五十页，共七十四页，2022年，8月28日边界条件根据问题的边界条件，节点的温度已知，即

根据显式表达式的稳定性判据，即0<F0≤1/4取F0=0.25时，根据傅立叶数公式，时间步长T1,j=1000j=1,2,3,4T5,j=860j=1,2,3,4Ti,4=1000i=2,3,4Ti,1=300i=2,3,4求解Δt=ρcΔ2F0/λ=22s故时间步长Δt不能超过22s。第五十一页，共七十四页，2022年，8月28日7有限差分方程的计算机解法含有n个未知数和n个方程的方程组，其一般形式为将方程简化为n为节点数，aij和bi都是常数（i=1,2,…,n，j=1,2,…,n）（91）第五十二页，共七十四页，2022年，8月28日[A]—矩阵，通常是方阵；{T}—节点温度列向量；{B}—右端列向量差分格式通常用矩阵形式表示[A]{T}={B}矩阵[A]是一个稀疏的带状矩阵，含有大量的0元素，而非0元素位于主对角线两侧。方程组有唯一解的充分必要条件是系数矩阵[A]是非奇异的，即它的行列式|A|≠0。第五十三页，共七十四页，2022年，8月28日迭代法一般比直接法好，系数矩阵是稀疏矩阵时，迭代解法更快些，并且没有舍入误差的积累。直接解法作有限次运算就可得到解，但是有舍入误差的积累，这种误差积累可能淹没计算的解。用计算机求解方程组的方法：直接法和迭代法。第五十四页，共七十四页，2022年，8月28日1.高斯消去法高斯消去法是线性代数方程组直接接法之一，其解题步骤分两步：第一步是消元过程，第二步是回代过程。假定把要求的n阶线性代数方程组写成如下形式（1）消元过程首先消去第二个方程以后的n-1个方程中的未知量T1。第五十五页，共七十四页，2022年，8月28日方法是：第二个方程减去第一个方程乘以a21(1)/a11(1)，第三个方程减去第一个方程乘以a31(1)/a11(1)，直到最后一个方程。这样方程组变为等价方程组方程中的系数aij(2)和bi(2)的计算公式为第五十六页，共七十四页，2022年，8月28日其次，从等价方程的第3个方程开始，消去后面n-2个方程中的未知量T2。即第三个方程减去第二个方程乘以a32(2)/a22(2)，第四个方程减去第二个方程乘以a42(2)/a22(2)，如此下去，即可得到方程中的系数计算公式为第五十七页，共七十四页，2022年，8月28日上述步骤重复n-1次后，可得到如下等价三角形方程组表示成矩阵形式为[A(n)]{T}={B(n)}式中，[A(n)]为上三角形式的矩阵，{B(n)}为n×1阶向量（1）第五十八页，共七十四页，2022年，8月28日可以把消元过程的计算归结为对于k=1,2,…,n-1的递推公式（2）该式就是消元过程的基本公式。（2）回代过程经过消元以后的方程组为式（1）。按照由下往上的顺序，可依次解出Tn,Tn-1,…,T1。第五十九页，共七十四页，2022年，8月28日由式（1）最后一个方程可得（3）将解得的Tn代回倒数第2式，可解得Tn-1。依次类推，由第i式可解得（4）式（3）和式（4）便是回代过程的基本公式。第六十页，共七十四页，2022年，8月28日2.迭代法迭代法的基本思想是，构造一个由{T1，T2，…，Tn}列向量序列，使其收敛于某个极限向量{T1*,T2*,…,Tn*}，而且{T1*,T2*,…,Tn*}就是方程组的精确解。根据构造列向量的方法不同，有简单迭代法、高斯—赛德尔迭代法和超松弛迭代法。迭代法不需要存储系数矩阵中的零元素，所以占用的存储单元少，程序也比较简单。其缺点是需要进行多次迭代才能达到收敛指标要求，耗费机时较多。第六十一页，共七十四页，2022年，8月28日（1）简单迭代法迭代的最终目的是求解方程组(91)中的T1,T2,…,Tn。若系数矩阵[A]对角线元素不为0，即aii≠0(i=1,2,…,n)，则可将式(91)改写为：其中任一方程均可写成任意给定Ti

(0)作为解的第零次近似，把