第八章-边界层方程的数值解

1第八章边界层方程的数值解对于边界层，第五章中给出了其微分方程和积分方程。这两类方程均可用数值解法求解。在高速数字计算机尚未普遍地得到应用之前，常采用积分方程的数值解。为了更准确地预测边界层内的计算结果，以及在电子计算机的速度和储存量都大大增加的有利情况下，对于边界层的微分方程也完全有能力进行数值计算。边界层的微分方程比积分方程给出了更为详细的物理描述，但需要把流体所受的剪切应力和流体的传热量分别与贯穿整个边界层的速度分布和温度分布联系起来，而不像积分方程那样，只牵涉到壁面的情况。对于层流，情况较为简单；可对于紊流，必需有一些假设来提供必须的关系式。如常用的假设为有效粘度模型，即yueff∂∂=μτ（8－1）其中，effμ称为有效动力粘度，teffρεμμ+=（8－2）层流时，0=tε；紊流时，必须由半经验关系式给出effμ，不同的模型给出的半经验关系式也不同，因而也有不同的数值解法。布腊德晓（Bradshaw）等提出了将有效粘度与紊流动能K相联系的模型：21Klkeffρμ=（8－3）其中，kl是距壁面距离的函数，K是流体脉动的平均动能)212121(222wvuK′+′+′=ρρρ。由于式（8－3）中包含有紊流动能K，因此必须联立求解紊流的能量方程和动量方程，因此，布腊德晓（Bradshaw）等提出了一个将紊流的能量方程和动量方程联立求解的数值解程序。另一种称为斯帕丁－帕坦克(Spalding－Patankar)的方法，是采用普朗特混合长度的模型，即假定yulmeff∂∂+=2ρμμ（8－4）其中，ml是普朗特混合长度。有了这种紊流模型后，便可将边界层的微分方程化成差分方程，然后进行数值计算。边界层微分方程的另一种方法是横流向积分，即在一个给定地点（给定的x），将偏微分方程用有限差分近似转换成常微分方程，再用迭代法进行数值积分。还有一种是参数积分法，将偏微分方程乘以自变量或他变量的权函数后，横越边界层进行积分（对y积分），得出一组常微分方程。这些方程可以用矩阵表示，矩阵的微分系数为未知量，经过矩阵逆变后用已知量表示，然后再解此方程组。这里为了便于大家对边界层偏微分方程数值解法的了解，主要介绍最常用的数值计算方法――斯帕丁－帕坦克(Spalding－Patankar)的方法。许多数值计算的方法现在都已编译成2了标准的计算软件，如Ansys软件、Fluent软件等。8—1不可压缩流体层流流动动量方程的有限差分解不可压缩流体在层流流过平板时的边界层的连续方程和动量方程（忽略体积力）为：0=∂∂+∂∂yvxu（8－5）221yudxdpyuvxuu∂∂+−=∂∂+∂∂νρ（8－6）边界条件为：0=y，0=u，0=vδ=y，∞=uux和y分别是顺流向和横流向的坐标，u和v分别是顺流向和横流向的分速度。如用矩形网格，则层流边界层的有限差分网点如图8－1所示。其中，xΔ和yΔ分别是顺流向和横流向的步长，均取为常数。对节点),(ji处的速度jiu,和jiv,分别用泰勒级数展开，可得：3,222,,,1)()(!2)()(xOxuxxuxuujijijijiΔ−∂∂Δ+∂∂Δ−=−3,222,,,2)()(!2)2()(2xOxuxxuxuujijijijiΔ−∂∂Δ+∂∂Δ−=−由此可得，2,2,1,)(2xOuuujijijiΔ+−=−−类似地，对于jiv,，可得2,2,1,)(2xOvvvjijijiΔ+−=−−同样，对1,+jiu、1,−jiu用泰勒级数展开，可得：3,222,,1,)()(!2)()(yOyuyyuyuujijijijiΔ+∂∂Δ+∂∂Δ+=+3,222,,1,)()(!2)()(yOyuyyuyuujijijijiΔ−∂∂Δ+∂∂Δ−=−对y的导数用中心差分公式，对x的导数用向后差分公式，则得：yx1,−jiji,1,+ji1,1−−jiji,1−1,1+−ji1,2−−jiji,2−1,2+−ji图8－1层流边界层的有限差分网点32,2,1,,)(243)(xOxuuuxujijijijiΔ−Δ+−=∂∂−−21,1,,)(2)(yOyuuyujijijiΔ+Δ−=∂∂−+及221,,1,,22)()(2)(yOyuuuyujijijijiΔ+Δ+−=∂∂−+将上述各有限差分代入式（8－6），原动量微分方程即成如下线性体系：jjijjijjijFuCuBuA=+++−1,,1,（8－7）其中，2,2,1)()2(2yxvvyxAjijijΔΔ−−ΔΔ−=−−ν2,2,1)(2)2(23yxuuBjijijΔΔ+−=−−ν2,2,1)()2(2yxvvyxCjijijΔΔ−−ΔΔ−=−−νijijijijijdxdpxuuuuF)()4)(2(21,2,1,2,1ρΔ−−−=−−−−对于横贯边界层的N个网络点（y方向除y=0和y=δ以外）而言，就有)2(−N个)1,,3,2(,−=NjujiL的线性代数方程组。因系数jA、jB、jC、jF中只包含已知量（即预先算出的或初始的u和v值，以及已规定的dxdp值），jiu,的各未知值即可由形成的矩阵式计算得出，边界条件为01,=iu，∞=uuNi,。得出jiu,后，可由连续方程求jiv,，dyxuvyji∫∂∂−=0,（8－8）这个积分可用梯形积分计算。被积函数由一群直线段来近似，每个线段下的面积即为梯形。故⎥⎦⎤⎢⎣⎡∂∂+∂∂Δ−=−−jijijijixuxuyvv,1,1,,)()(2（8－9）其中，01,=iv。在给定的i处，如1,iu及下游站的速度已确定后，可用一个线性插值函数来计算壁面处的切应力。434,3,2,)()2918(6)(yOuuuyyuiiissΔ++−Δ=∂∂=μμτ（8－10）此法带有一个重要的缺点，即坐标系（笛卡尔直角坐标系或极坐标系）一般不能和边界层的外形相符（y=δ，而δ随x变化）。当边界层增长时，必须添加更多的网络点，或从一开始就带有额外的网点，这样就使得在下游方向每走一步，都要核对横流向的每一步的jiu,值的变化，来试出边界层的外缘，因此在计算时间和储存方面都是不经济的。更好的办法是用一个随着边界层增长而增长的坐标系。8—2流线坐标系对于流体沿着平板的流动，如喷嘴或扩压器内流体的流动、沿对称物体的流动、圆管和扁管内的流动以及各种自由射流均可建立轴对称极坐标系，使问题简化。对于稳态、可压缩流体的流动，用轴对称极坐标系时，其连续方程和动量方程（忽略体积力）为：0)()(=∂∂+∂∂yvrxurρρ（8－11）yrrdxdpyuvxuu∂∂+−=∂∂+∂∂)(1)(τρ（8－12）其中，x和y仍分别是顺流向和横流向的距离，u和v仍分别是顺流向和横流向的分速度。r为离对称轴的半径，1=r时，上式即简化为笛卡尔坐标系的方程。采用轴对称坐标，仍避免不了矩形网格的缺点，为此需采用轴对称的流线坐标系。引进流函数ψ作为横流向的自变量，顺流向的坐标用ζ，x≡ζ。这样，以),(yx表示的轴对称坐标系就变成了以),(ζψ表示的轴对称流线坐标系了。这种坐标变换称为冯·米塞兹变换。由于引入了流函数ψ，则yur∂∂=ψρ，xvr∂∂−=ψρ。显然，流函数能自动满足连续方程，同时可在边界层动量方程中除去横流向的速度v，使之简化。由于坐标的变换，则式（8－12）中的xu∂∂、yu∂∂等项成为：ζψρζζψψ∂∂+∂∂−=∂∂∂∂+∂∂∂∂=∂∂uuvrxuxuxuψρζζψψ∂∂=∂∂∂∂+∂∂∂∂=∂∂uuryuyuyu5ζddpdxdp=（忽略横流向压力变化，即忽略ψddp）ψτρτ∂∂=∂∂)()(1ruyrr经过转换，式（8－12）就成为：ζρψτζddpuru1)(−∂∂=∂∂当把yueff∂∂=μτ代入后，即得：ζρψμρψζddpuuurueff1)(2−∂∂∂∂=∂∂（8－13）这就是以轴对称流线坐标系表示的动量方程。同样，可将能量方程222)()(yuytytvxtucp∂∂+∂∂=∂∂+∂∂μλρ也改写成轴对称流线坐标系中的形式：⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎥⎥⎦⎤⎢⎢⎣⎡∂∂−+∂∂∂∂=∂∂)21()Pr(~Pr~22uhurheffeffeffeffeffψμμψμρψζ（8－14）其中，221~uTChp+≡称为全焓，effPr为有效普朗特数。经过这样的变换后，边界层的动量方程和能量方程具有相同的形式，即同为抛物型的偏微分方程，因此能用相同的数值计算程序来求解。事实上，若对质量扩散方程2*2**yCDyCvxCujjj∂∂=∂∂+∂∂（C为质量浓度，D为质扩散系数）。进行变换，也是变成抛物型的偏微分方程，也能用该数值计算程序来求解，如下式：)(*2*uRCSeurCjjeffeffjρψμρψζ+∂∂∂∂=∂∂（8－15）其中，jR为由于化学反应等原因组份j的产生率，effSe为有效斯密特数。在对边界层微分方程进行数值计算时，为了使边界层的有限差分网格会自动增长或收敛，与定义的边界层外缘相一致，这里特定义一个无量纲流函数η，IEIψψψψη−−≡r1rφηψηη,,yζ,x对称轴线等ηI边界，0=ηE边界，1=η图8－2无量纲的流线坐标系6其中，Iψ和Eψ分别为所考虑的内边界I和外边界E处的流函数，通常Iψ和Eψ是ζ的函数。这样，顺流向的变量为ζ，而横流向的变量则为η，等η线与等ζ正交。边界层位于0=η和1=η之间。如图8－2所示，其中Ir为内边界离对称轴的距离，r为边界内的某一点离对称轴的距离，φ为母线与对称轴的夹角。由于引入了无量纲流函数η，式（8－13）、（8－14）和（8－15）均可用),(ηζ坐标表示成：dcba+∂Φ∂∂∂=∂Φ∂++∂Φ∂)()(ηηηηζ（8－16）此处，Φ为一充代变量，可以是u、h~或*jC，系数分别定义为：IEIImraψψ−≡&IEIIEEmrmrbψψ−−≡&&22)(IEeffeffurcψψσμρ−≡表8－1给出了和Φ相应的effσ和称为源项的d值。表8－1和Φ相应的effσ和称为源项的d值他变量Φeffσdu1ζρddpu11−h~effPr⎥⎥⎦⎤⎢⎢⎣⎡∂∂−−∂∂)21()Pr()(222uureffeffeffIEημμψψρη*jCeffSeuRjρ8—3边界层方程的有限差分形式7将基本的偏微分方程用式（8－16）的形式表示后，需将转换成有限差分方程，以便求出数值解。一、网格的划分用数值法计算边界层内的动量、热量和质量交换情况时，必须把计算区域划分成网格，并沿横流向解有限差分方程。网点之间的距离，亦即网格的尺寸应根据边界层内的速度或焓、浓度等的分布规律来确定。当横流向的速度或焓、浓度等的梯度不大时，网格可以划分得较稀一些，相反，如果横流向的梯度很大，则需要密集的网格。仔细分析紧靠壁面处流体的流动与传热和传质情况，可以发现，不论是速度分布，还是焓的分布，在紧靠壁面处，它们的梯度较大，因而需要将该处进行特殊的处理。一般有两种处理方法可供选择。一种是“越过壁面函数”（BypasstheWallFunction）法。在靠近壁面处与壁面之间，用逐渐变细的网格，通常由于紧靠壁面处的梯度很大，以致在紧靠壁面20％的边界层厚度内的网点数需与其余80％厚度内的网点数相同，从而大大增加了需用的计算机储存量和计算的时间。这种方法用于压力梯度较大的场合。另一种方法是“壁面函数”（theWallFunction）法。在整个高速度梯度区利用库特流公式进行数值计算。这种方法的主要优点是大大减少近壁面区域的有限差分网格数，特别在计算高雷诺数的流动时，用“壁面函数”法特别方便。这里用“壁面函数”法来划分网格。在任何顺流向位置（=ζ常数）处，网格线将0=η和1=η之间的区域划分成N个条带。由于在边界处的速度梯度和温度梯度比内部区域大得多，因此，在0=η处设置1η、2η两个点，使021==ηη，1η和2η间的条带称为I边界。同样，令112==++NNηη，1+Nη和2+Nη间的条带称为E边界，而2η和1+Nη间的)2(−N个条带则称为内部区域，如图8－3所示。二、内部区域的有限差分方程对于内部区域，这里将通过