传热与流体流动的数值计算(4-5章)

传热与流体流动的数值计算[美]S.V.帕坦卡著同济大学机械工程学院朱彤第四章热传导4-1本章的对象着手构建一个求解通用微分方程的数值方法构成一个求解通用微分方程的数值方法，略去对流项。其他一些物理过程也由非常类似于热传导方程的数学方程所控制。本章完成了随后几章所需要的若干预备性的工作，提出代数方程的求解方法。基本方程稳态一维问题的控制微分方程：推导出离散化方程4-2一维稳态热传导-网格间距网格点距离(δx)e与(δx)w没有必要相等。虽然只有在网格相当细时才可能得到精确的解，但是在因变量随x变化相当慢的区域没有必要采用细的网格；在T~x变化较陡的区域则需要细的网格。误区：不均匀网格的准确度比均匀网格差。设计一个合适的非均匀网格：从解的定性预计得到指导。用初步粗网格的解求得T~x变化形式；然后构成合适的非均匀网格。先进行预备性的实验或探索性试验，然后应用得到的数据资料确定在最终的实验中所应安装的探头位置和数目。-界面导热系数ke最直截了当的方法是假设k在P点和E值间呈线性变化：其中插入因子在某些情况下这种简单化会导致相当不准确的结果；而且这样做不可能精确处理组合材料中可能遇到的导热系数的突然变化。一种替代方法：得到一个通过下式描述的界面热流密度qe的良好达式：(4.5)(4.6)(4.7)讨论这样一种情况：围绕着网格点P的控制容积由具有均匀导热系数kP的材料填满，围绕着E点的控制容积由导热系数kE的材料填满，对于P点和E点之间的组合板，根据稳态无内热源一维导热的分析，有：合并得：当界面l位于P和E之间的中点时，有fe=0.5，有：上式说明ke是kP和kE的调和平均值，而非给出的平均值。(4.8)(4.9)(4.10)应用于系数的定义式，得到aE：其效能可由两种极限情况看出：令kE0，则有ke0（4.12），即一个绝热层表面上的热流密度为0。令kP>>kE，那么ke＝kE/fe（4.13）。表明界面的导热系数ke完全与kP无关；ke不等于kE，而是它的1/fe。目的是通过（4.7）得到一个正确的qe，应用(4.13)，得：当kP>>kE时，温度Tp将一直扩展到界面e，而温降Tp-TE将实际上发生在距离(x)e+内。两个极限情况讨论表明这个公式可以适用于导热系数突然变化的情况，而无需在发生突变的邻近区域采用极细的网格。(4.11)(4.14)-非线性即便是在热传导问题中我们也经常遇到非线性的情况。如离散化方程中的系数本身与T有关。我们用迭代的方法来处理。过程包括：一开始在所有各个网格点上，猜测或估计一个T值。由这些估计的T值，计算出离散化方程中的系数的试探值。解名义上的线性化方程组，得到一组新的T值。以这些T值作为较好的估计值，返回到第二步并重复整个过程，直到这种进一步的重复计算（迭代）不再引起T值任何有意义的变化为止。这种最终不变的状态叫做迭代的收敛。与之相反，迭代永远也不会收敛到一个解的状态称为发散。-源项的线性化当源项S与T有关时，用方程（4.4）给出的线性形式表达。当S是T的一个非线性函数时，必须把它线性化，即规定SC和SP的值。有很多方法可以把给定的S表达式分解成SC和SP。如：已知S=4-5T3。某些可能的线性化：1.SC=4-5Tp*3，Sp=0。这种做法不能很好利用已知S~T关系的有利条件。2.SC=4，Sp=-Tp*3。看起来像准确的线性化，但已知的曲线比这一关系所反映的曲线要陡。3.推荐的方法：在点Tp*，所选择的直线与S~T曲线相切。4.SC=4+20Tp*3，Sp=-25Tp*2。这一线性化比已知的S~T曲线陡，使收敛速度降低。四种可能的线性化与实际曲线比较如图：-边界条件讨论图中所示网格点组。在两个边界上各有一个网格点。其余网格点称为内点。围绕每个内点有一个控制容积。对每一个控制容积可以写一个像方程（4.2）那样的离散化方程，如果看作是关于Tp的方程，那么就有了对所有内网格点上未知温度所必要的方程。其中有两个方程包含着边界网格点上的温度。通过处理这些边界温度，就把已知的边界条件引入到数值解法中。热传导问题中有三类典型边界条件，对于每一种有：1.已知边界温度。此时不需要外加任何方程。2.已知边界热流密度。得到：(4.16)Boundaryconditions如果边界上的热流密度qB已知，则要求的对TB的方程变成：3.通过放热系数和周围流体的温度来规定边界的热流密度。如果热流密度qB是放热系数h以及环境流体温度Tt规定，那么，方程TB方程变为：这样就构成了对所有未知温度的足够数量的方程。-线性代数方程的解一维离散化方程的解可以用的高斯消去法得到。当写这些方程的系数矩阵时，所有的非零系数均排列在矩阵的三角对角线上，这种算法称为TDMA（三对角矩阵算法）。设网格点标号为1，2，3，…，N，其中1和N代表边界点。有：边界温度已知时，对边界点的方程只剩下一个无意义的形式。(如c1=0,bN=0)T2可以用T3的一个关系式表示，….，TN可以由TN+1表示，回代就是TDMA的要点。(4.22)TDMAwhathappenswhenaboundarytemperatureisgiven?TDMATDMATDMA算法计算：P1=b1/a1andQ1=d1/a1使用迭代关系式，获得Pi和Qi，i=2,3…N.设TN=QN使用迭代关系式，得到Ti=PiTi+1+Qi，i=N-1,N-2…3,2,1从而依次得到TN-1,TN-2,…T3,T2,T1.通用的离散化方程时间是一个单向坐标，由一已知的初始温度分布开始，沿着时间坐标逐步向前求解：已知t时刻T在网格点上的值，求得t+Δt时刻值。对整个控制容积积分方程：得到：4-3非稳态一维热传导假设在网格点上的T值代表整个控制容积上的值，最后得到：假设用下式归纳一般化有关TP、TE和TW如何随时间由t到t+Δt而变化的关系：其中f是在0和1之间变化的加权因子。于是：改写后得：-显式，克兰克-尼科尔森模式，以及全隐式模式对于某个特定的加权因子f的值，离散化方程可简化为适用于抛物线型微分方程的我们所熟悉的模式之一。f=0导致显式模式；f=0.5导致克兰克-尼科尔森模式；f=1导致全隐式模式。不同f值可以由图中所示关系来说明：显式模式假设老的值代表除了时刻t+Δt以外整个时间间隔上的Tp值；全隐式模式假设在时刻t，Tp的值突然降了，而后整个时间步上保持为降后的值，于是整个时间步期间温度为新值所确定。克兰克-尼科尔森模式假设Tp呈线性变化。如果我们要求方程（4.36）中系数务必永不为负，只有f=1。即全隐式模式能够满足我们简单而物理上又满意的要求。f=1显式格式ExplicitschemeinordertogiverealisticsolutionsCrank-Nicolson格式cangiveunrealisticsolutions隐式格式Implicitschemealwaysgivesrealisticsolutions-全隐式离散化方程线性化源项结果：Δt趋近于无穷大时，这个方程简化为稳态的离散化方程。全隐式模式主要原则：Tp的新值代表整个时间步上的值。因此如果导热系数kp与温度有关，就应当反复由它迭代算得新值。稳态程序的其他环节，如边界条件、源项线性化处理以及TDMA也都完全适用于不稳态问题。4.4Unsteady2-DheatconductionDiscretizedunsteady2-DheatconductionequationUnsteady3-DheatconductionequationDiscretizedunsteady3-Dheatconductionequation-三维问题的离散化方程加入两个z方向的相邻点T和B（项和底）构成三维的网格图形。相邻系数aE、aW、aN…、aB代表P点与相邻点之间的热导；a0PT0P是t时刻控制容积内部所包含的内能（除以Δt）常数b由这一内能项与由Sc所造成的在控制容积内的发热率组成。中心点系数ap是所有相邻点系数之和，并包括一项由线性的源项所作的贡献。-代数方程的解迭代法高斯-赛德尔逐点计算法按一定的顺序逐个访问每一个网格点，以计算那里的变量值。在计算机内值需要存储一组T值。开始，这些值代表最初的估计值或上一次迭代得到的值，在访问每一个网格结点时，在计算机存储中相应的T值交替改变。这种方法不是总可以得到收敛解的。斯卡巴勒准则：高斯-赛德尔法收敛的充分条件是：中，这种方法的主要缺点是收敛速度太慢，特别是网格点数很大时。逐行法把TDMA和高斯-赛德尔法结合起来。选择一条网格行（设在y方向选取这样的网格行），假定沿相邻的行上的T值批最新值构成。用TDMA法求得所选行上的T值。将在同一方向的所有行进行这种计算。如果想做的话，再按相同的方法在其他方向重复上述程序。以二维为例，如图所示的情况需要注意：其它一些迭代方法ADI（方向交替的隐式）的逐行求解法；解多维离散化方程的强隐式法（SIP）依前后二次迭代之间因变量的变化究竟是被加速还是被减慢的过程称为超松弛或欠松弛。超松弛常用于和高斯-赛德尔法相结合，叫做持续超松弛（SOR）；欠松弛在强烈非线性方程组的迭代求解中用来避免发散。取T*p作为前一次迭代所得Tp值。引进松弛因子，得到：可以根据经验以及对所给定的问题所作的试探性计算求得一个合适的值。4-5超松弛和欠松弛(overrelaxationandunderrelaxation)通用惯量进行松弛。用下面公式代替离散化方程：式中i是所谓的惯量。对于正的i值，方程具有欠松弛作用；对于负的i则产生超松弛。控制容积面的位置讨论控制容积面构成的两种不同的替代形式，并讨论它们各自有关的优点。为方便起见，描述针对二维问题。方法A：控制容积面放在两个网格之间的中点。4-6某些几何上的考虑结果是：一个典型的网格P并不落在包围该点得控制容积的几何中心上。方法B：网格点放在控制容积的中心：克服了A的缺点。具有方便性。我们所提出的这种方法不只限于直角坐标系，还可以用于任意一种正交坐标系。以二维极坐标问题为例，与方程对应的r形式是：其中的网格与控制容积如图示：设控制容积在z方向厚度为1，方程两边乘以r，并在整个控制容积范围内对r和进行积分，得到下面的离散化方程：4-6其它坐标系由一个新的坐标系引入的补充特征主要是几何上的特征。5-1任务在通用微分方程中将对流项考虑进去，只要对流项的加入不改变离散化的形式，同样的处理方法仍然适用。本章任务是：在已知的流场（即速度分量和密度）的情况下，求得对φ的解。已知流场必须满足连续性方程：第五章对流与扩散通用微分方程也可以改写为：对于已知的、uj、以及S的分布，任何解φ及其变体（φ加一常量）将同时满足方程，关于系数和的基本原则仍然适用。讨论只有对流与扩散这两项存在的情况下的一维稳态问题。控制微分方程：应用图示三网点群：5-2一维稳态对流与扩散-预备性的推导在整个控制容积内对方程（5.4）积分：由对φ的一个分段线性分布表示项Γdφ/dx。结果是：因子1/2出自界面位于中点的假设；对不同的界面位置要采用其它内插因子。则方程（5.6）写成：定义两个新的符号：两者具有相同因次，F表示对流或流动的强度；D是扩散传导性。（注意，D永远为正，而F不同）离散化方程变为:讨论由于连续性Fe=Fw，得到ap=aE+aW上述离散化方程隐含着φ分段线性分布——中心差分假定De=Dw=1及Fe=Fw＝4，则若φE＝200及φW＝100，φP＝50若φE＝100及φW＝200，φP＝250方程（5.11)表明系数可能出现负值。当|F|小于2D时，系数才可能始终为正，即中心差分只能限于低Reynold数流动。若扩散项为零，则中心差分格式导致ap=0，无法使用逐点法求解-上风亦称为上风差分格式、迎风格式、上游差分格式以及供体（施主）室法等。方案认为预备性公式的弱点在于假设：界面上的对流性质φe是φE和φP的平均值。提出：保留扩散项的公式不变，而对流项则按下列假设计算：界面上φ的值等于界面上风侧网格点上的φ值。于是类似方法可确定φw值。定义代表A，B中大者。则上风方案意味：离散化方程可写为：不会产生负的系数；可以把这个方案说成是建立在“槽与管”的模型基础上，管内的流体不会“知道”将要流入那个槽内的任何情况，但它却携带了它所来自那个槽内的全部信息。这就是上风方案的本质。精确解(Exactsolution)如果Γ取作常数，且边界条件为：则其中为贝克列数（PecletNumber），是对流与扩散强度之比。不同的贝克列数时的φ~x变化如图P为0的极限条件下，问题成了纯扩散（或热传导）问题。除非lPllPl值非常小，曲线均偏离线性很远。当lPl值很大时，符合上风方案假设，但上风方案假设用于所有lPl。当lPl值很大时，x=L/2处d/dx几乎为0，扩散几乎不存在。上风方案总是由一线性的～x分布计算扩散项，从而在大的lPl值条件下过高估计了扩散项。-指数方案由对流流量密度和扩散流量密度所组成的总流量密度由dJ/dx=0，在整个控制容积内积分，得到：可推得Je表达式：同理得到Jw表达式代入，写成标准形式：在应用于一维稳态问题时，该指数方案保证得到精确解。但：费时；对于二维、三维及源项不为零时，不准确。-混合方案由可以看到aE/De准确变化的某些特殊性质：代表这些极限情况的三条直线如图中所示，它们构成准确曲线的一根包络，并代表着这一准确曲线的合理近似。混合方案实际就是有着三条直线组成。用特殊符号代表其中包含的所有量的最大值。于是在贝克列数为-2≤Pe≤2时，混合方案同中心差分格式一致；在该范围之外，混合方案简化为上风方案。可把混合方案的对流-扩散离散化方程写成：-幂函数方案在Pe=±2时，混合方案偏离准确曲线相当大，一个更好的近似由幂函数方案给定。aE的幂函数表达式可以写成：紧凑形式可以写作：幂函数与准确的指数方案之间的差异非常小。-一个通用化的公式讨论图中所示由距离δ分开的网格点i和i+1。有：建议：α、β是与P有关的无因次乘数，这样，就有：其中A、B是无因次系数。它们是贝克列数的函数。如果φi和φi+1相等，扩散流为0，就有B=A+P。如果将坐标轴方向反转，有A(-P)=B(P)或B(-P)=A(P)。A和B随贝克列数P的准确变化如图所示。应用流量关系式（5.37）于界面e和w，并利用方程得到通用的对流扩散公式：于是可以把前面所推得的各种方案看成是选择不同的函数A(lPl)而已。如图、表所示：-各种方案（格式）的结果在结束一维问题的讨论之前，检查一下对于φE和φW由各个方案所计算出来的φP值。令φE=1，φW=0，(δx)e与(δx)w相等，于是φP将是P的函数。如图所示：讨论图中的控制容积。5-3二维问题的离散化方程推导的细节方程（5.2）的二维形式：其中在控制容积内对方程（5.48）积分，得到：类似积分连续性方程得到：后面四项示通过控制容积面的质量流量。以φP乘以方程（5.51）并从方程（5.50）中减去所乘结果，得：在控制容积面上均匀性的假设使我们可以利用一维做法解决二维问题，有：-最终的离散化方程二维的离散化方程写成：相应的传导性定义为：贝克列数定义为：幂函数方案是推荐的，有：5-4三维问题的离散化方程流量与传导性定义为：贝克列数取为F/D，于是Pe=Fe/De，依此类推。幂函数公式为：5-5单向空间坐标使空间坐标成为单向坐标的条件当贝克列数大时，下游相邻点的系数变小，当贝克列数超过10时，幂函数方案将取下游相邻点系数为0。这样，由于在任何点上的φ值将不受x方向下游值的影响，x就成为一个单向的坐标。即便一个空间坐标就整个计算域而言并不是单向的，但在实际处理边界条件时，往往应用其局部的单向特性。出流边界条件在流动出口的边界上，不需要有关边界条件的任何信息。如果贝克列数足够大，系数aE将为0，因而系数乘边界值为0。即在出流边界附近的区域，对大的贝克列数而言，呈现局部的单向状态。5-6假扩散关于假扩散的一般观点中心差分格式具有二阶的精度，而上风方案只具有一阶的精度；上风方案引起严重的假扩散由于在对流-扩散问题中所产生的φ~x变化是指数的，除了极小的Δx值外，任何情况下泰勒级数不再是一个好的表达式。较大Δx值时，泰勒级数的分析给人误解；此时上风方案给出比中心差分格式更为合理的结果。把中心差分格式同上风方案的系数比较可以看到，上风方案似乎把真正的扩散系数扩大了一个虚拟的扩散系数ρuδx/2。但需要注意的是，争辩的麻烦建立在把中心差分格式假设为精确而又标准的基础上，并且用这个参考来观察上风方案。按这种做法，就会发现本身就是精确解的指数方案也有某种假扩散。实际上所谓的假扩散系数是在大的贝克列数条件下的一种理想补充。有关假扩散的正确看法假扩散是一种多维的现象，在稳态的一维问题中绝对不会有相应形式的假扩散。为了具体化，讨论图中所示状态。两股速度相等而温度不同的平行流相遇。如果扩散系数不等于0，就形成一温度逐渐由高温到低温变化的混合层；如果扩散系数不为0，就不会形成混合层，并且在流线方向将维持温度得不连续性。如果对Γ=0情况的数值解产生一个逐渐变化的温度分布，我们可以得出结论：该数值解方案引起假扩散。对两种不同网格方向用上风方案求解5.13b中的问题。在x方向的均匀流动。Φp=φW，结果，每一条水平线上上游的已知值将构成在该线上所有点上的值。没有假扩散现象发生。在与网格线成45。角的方向上的均匀流动。有设左边界温度为100，底边界温度为0来表示不连续。在内点上所得的解写在每个网格点的旁边。所得到的实际值却代表着一个逐渐变化的温度分布。要点流场方向与网格点线成一倾斜角，并在与流动方向相垂直的方向上存在有非零的因变量梯度时，就会有假扩散出现。对二维情况的假扩散系数的一个近似表达式由公式给定：其中，Γ假是假扩散系数，U是合速度，θ是速度向量与x方向间夹角。可以看出：合流动方向与其中一组网格线相重时，不存在假扩散；流动方向与网格线间夹角45。时，假扩散最严重。减少Δx和Δy可以减少假扩散大小，只要可能，应当把网格线布置在接近于流动的方向。应使假扩散与实际值相比足够小。采用中心差分格式不是解决假扩散的灵丹妙药。产生假扩散的基本原因是由于把流过每一控制容积面的流动处理成局部一维流动这样的一种做法。那些可能给出较小假扩散的方案应当考虑流动的多维性质。