不可压库埃特流的数值解

第9章不可压库埃特流的数值解主讲人：彭赛，韩昭辉时间：2016.5.8•9.1引言•局限性：•1）之前第7,8章的数值方法是双曲型偏微分方程的显示解法，这限制了推进的步长；•2）之前讲的都是无粘的流动。•本章的改进或不同之处：•1）求解控制方程的方法是隐式的差分解法；•2）问题的控制方程为抛物线方程；•3）考虑的问题是粘性的流动。•本章具体考虑的是不可压库埃特流动。•1）有解析解；•2）粘性流动，与边界层的流动具有很多相似的物理性质；•3）提出压力修正的方法处理二维不可压的N-S方程。•9.2物理问题及其解析解•平行平板的间距为,•上平板以运动，•下平板静止，，•粘性流动。•不可压流体的质量守恒方程为•(9-1)•如果我们认为方向无限长，则有•(9-2)•则可以将(9-1)式简化为•(9-3)•而上下边界处有•(9-4)•故整个流场中都有•(9-5)•eu0u0uvxyx0x0vy000vvD0vD•考虑方向的动量方程•(9-6)•不考虑体积力，则可以简化为•(9-7)•而牛顿流体切应力与速度的关系有•(9-8)•因此(9-7)式可以简化为•(9-9)•所以库埃特流动，在和方向都没有压力梯度。•考虑方向的动量方程•(9-10)•由于(9-11)•(9-12)•将(9-11)和(9-12)带入(9-9)有•(9-13)yxyyyyzyDvpfDtyxyz0yypyy20yyuvvxyy0pyxyxxxuupuvxyxxyxy20xxuvuxyxyxvuuxyy0uyy•流动不可压、恒温的，则常数，有•(9-14)•对积分两次有•(9-15)•其中和是积分常数，考虑边界条件，则有•(9-16)220uyy12ucyc1c2ceuyuD•9.3数值方法：隐式克兰克—尼克尔森方法•初始速度剖面设定•(9-17)•速度剖面随时间演化如图所示。•经过足够长的时间，速度剖面的形状趋近于9-2d所示的定常流。•图9-2中，速度剖面随时间变化的流动称为非定常的库埃特流动。这里任然假设和，沿方向的非定常库埃特流动，控制方程为(9-18)这是一个抛物线型的偏微分方程，可以在时间上推进求解。0for0fory=DeyDuu/0x0vx22uuty•9.3.1数值方法•把(9-18)表示成无量纲的形式。定义如下的无量纲变量••把方程(9-18)无量纲化为••或(9-19)•在(9-19)中，我们发现••是按照两平板之间距离计算的雷诺数，于是方程(9-19)变为•(9-20)•省略“ˊ”号，(9-20)可改写为•(9-21)•根据克兰克-尼克尔森格式，方程(9-21)的有限差分表达式为••euuuyyD/ettDu2222/////eeeeeuuuuuuDDtDuyD22euutuDy1ReeDuDReDD221ReDuuty221ReDuuty111111112111221222RennnnnnnnjjjjjjjjDuuuuuuuuty•(9-22)•方程(9-22)可以表示成如下形式：•(9-23)•其中(9-24)•(9-25)•(9-26)•求解方程(9-24)的网格如图9-3所示。•(9-27)•和可以通过边界条件得到•(9-28a)•(9-28b)•由式(9-23)可以得到一个方程组•1111111222221112ReRe2ReRe2RennnnnnjjjjjjDDDDDtttttuuuuuuyyyyy11111nnnjjjjAuBuAuK22ReDtAy12BA1112nnnjjjjKAuAuuDyN1u1Nu10u10Nu•第一个方程组是•(9-29)•因为，因此方程(9-29)变为•(9-30)•式(9-23)所代表的最后一个方程为•(9-31)•因为。因此方程(9-31)变为•(9-32)•于是，由(9-23)所代表的方程组可以用矩阵的形式表示为：1111232nnnAuBuAuK10u11232nnBuAuK11111nnnNNNNAuBuAuK11Nu111nnNNNAuBuKA•(9-33)•这是一个三对角的方程组，可以用托马斯算法。•把托马斯的算法应用到式(9-33)所表示的方程组上，将得到，•，…，。它们是时刻的速度值。然后，反复进行这个过程，直到速度剖面收敛到一个稳定的状态，如图9-2所示。•9.3.2问题的提法•这里取了21个网格点，因此••本问题的特点，稳定性好，但是我们如果要研究开始流动的瞬时变化，那么应取得小一些，以降低时间上的截断误差。12213314415511110000000000000000000000000...............0000000000000nnnnnNNnNNBAKABAuKABAuKABAuKuKuKAABAuAB12nu13nu1nNu1n120yt•相应的显示算法的稳定性条件是•(9-34)•也就是(9-35)•在这里我们取(9-36)•来通过数值的实验考察在这个范围内得到的结果。•(9-37)•将(9-37)式带入(9-24~26),我们得到•(9-38)•(9-39)•(9-40)•我们通过解析解可以发现最终稳态的时候速度剖面与是没有关系的，但趋近定常状态的瞬时过程却是依赖于的。这里定义了参数,排除了雷诺数的影响，得到更本质的特性。211Re2Dtx21Re2Dty2Re14000DtEyEE2ReDtEy2EA1BE1112nnnjjjjEKEuuuReDReDE•9.3.3中间结果•取和,考察1个时间步之后的速度剖面•9.3.4最终的结果1ERe5000D•表9-2后面几个时间步的速度剖面•这里我们都是将,我们下面考虑取更大的时间步长会有什么影响？•1)从稳定性的角度来看，克兰克-尼克尔森方法是无条件稳定的；•2)时间步长的增加，会使得我们在离散化的时候产生更大的截断误差，降低瞬时结果的精度。•我们下面看看运算过程是不是这样。1E•表9-3中给出了,无量纲化的时刻的速度剖面，稳态时刻。•可以看出时间步长越长，我们计算的瞬态速度精度越差。1,5,10E31.510t34.510t•我们下面考虑一个极端的情况。••圆点是解析解，定常的解，我们可以看到，在这种极端情况下，瞬态速度已经严重偏离精确解，已经没有物理上的意义了。•所以我们希望考察一个物理变化瞬时结果,我们要将时间步长设定为比较小。4000E•9.4另一种数值方法：压力修正法•之前我们考虑板是无限延伸的，但是计算的区域有限的，计算区域如下图所示：•边界条件设定：在入口处，和固定，可以浮动，出口处，只有压强固定。•压力修正法是一种pvup