流体力学论文

C++程序实现非定常COUETTE流速度模拟及古典显式有限差分格式收敛性验证院系名称能源与动力工程学院专业名称飞行器动力工程作者姓名洪启臻(13041222)2015年5月13日C++程序实现非定常COUETTE流速度模拟及古典显式有限差分格式收敛性验证摘要非定常Couette流（库艾特流）是指两平行平板间的非定常流，是流体力学中简单而又基础的一种流动方式。有限差分法是CFD(计算流体力学）中常用的方法，其中本文考虑的是扩散式抛物线方程，即热传导方程，运用古典显式有限差分格式，也就是时间导数采用向前差分，空间二阶导数采用中心差分来研究Couette流。进而通过C++编程来进行Y方向上的速度分布模拟，通过EXCEL做成表格。分别考察了△t等于20,200,2000的情况，得出相应结果。最后，就古典显式有限差分格式的收敛性给出了讨论，当系数u大于0.5时，求解将不再收敛，进而验证了数值分析中VonNeumann得出的结论。关键字：Couette流（库艾特流）C++古典显式有限差分收敛性一．引言1.1基础知识[1]下面，我们先引入流体力学的一些基本假设及控制方程。我们认为，流体质点是几何尺寸同流动空间相比是极小量，又含有大量分子的微元体。连续介质则是质点连续地充满所占空间的流体或固体。进而提出连续介质模型，即把流体视为没有间隙地充满它所占据的整个空间的一种连续介质，且其所有的物理量都是空间坐标和时间的连续函数的一种假设模型：u=u(t,x,y,z)。流体的压缩性是用体积压缩率k来量度d/d/ddVVkpp(1-1)式中：p为外部压强。在研究流体流动过程中，若考虑到流体的压缩性，则称为可压缩流动，相应地称流体为可压缩流体，例如高速流动的气体。若不考虑流体的压缩性，则称为不可压缩流动，相应地称流体为不可压缩流体，如水、油等。本文考虑的流体均为不可压缩流体。流体的粘性指在运动的状态下，流体所产生的抵抗剪切变形的性质。粘性大小由粘度来量度。流体的粘度是由流动流体的内聚力和分子的动量交换所引起的。粘度有动力粘度和运动粘度之分。动力粘度由牛顿内摩擦定律导出：dduy(1-2)式中：为切应力，Pa；为动力粘度，Pas；d/duy为流体的剪切变形速率。根据流体是否满足牛顿内摩擦定律，将流体分为牛顿流体和非牛顿流体。牛顿流体严格满足牛顿内摩擦定律且保持为常数。非牛顿流体的切应力与速度梯度不成正比，一般又分为塑性流体、假塑性流体、胀塑性流体3种。本文考虑的流体均为牛顿流体。描述流体物理量有两种方法，一种是拉格朗日描述；一种是欧拉描述。由于本文需要只介绍欧拉描述。欧拉描述，也称空间描述，它着眼于空间点，认为流体的物理量随空间点及时间而变化，即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为(q1,q2,q3)，用欧拉坐标表示的各空间点上的流体物理量如速度、压强等，在任一时刻t的值，可写为q1、q2、q3及t的函数。从数学分析知道，当某时刻一个物理量在空间的分布一旦确定，该物理量在此空间形成一个场。因此，欧拉描述实际上描述了一个个物理量的场。若以f表示流体的一个物理量，其欧拉描述的数学表达式是(设空间坐标取用直角坐标)(,,,)(,)fFxyztFrt(1-3)进一步的，根据流体流动过程以及流动过程中的流体的物理参数是否与时间相关，可将流动分为定常流动与非定常流动。定常流动是流体流动过程中各物理量均与时间无关，这种流动称为定常流动。非定常流动是流体流动过程中某个或某些物理量与时间有关，则这种流动称为非定常流动。本文考虑的是非定常Couette流，即其流动过程中物理量与时间有关。流体流动所遵循的物理定律，是建立流体运动基本方程组的依据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考虑不同的流态，如层流与湍流。在流体力学中，系统是指某一确定流体质点集合的总体。系统与外界无质量交换，但可以有力的相互作用，及能量(热和功)交换。控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动形状任意的空间体积。包围这个空间体积的边界面，称为控制面。控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能量交换，也可与控制体外的环境有力的相互作用。下面，我们引入计算流体力学中的控制方程，我们仅介绍与本文有关的方程。质量守恒方程(连续性方程)：质量的守恒写作：(1-4)其中ρ是流体的密度。在不可压缩流体的情况ρ不是时间或空间的函数。方程简化为：(1-5)也就是0uvwxyz(1-6)动量守恒方程(运动方程)：动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或N-S方程，其微分形式表达如下：ddddddyxxxzxbxxyyyzybyyzxzzzbzpppuFtxyzpppvFtxyzpppwFtxyz(1-7)式中：bxF、byF、bzF分别是单位质量流体上的质量力在三个方向上的分量；yxp是流体内应力张量的分量。若不采用简化书写的完整形式非常繁琐，分别为：上组就是著名的N-S方程，以克劳德-路易·纳维和乔治·加布里埃尔·斯托克斯命名，这些方程建立了流体的粒子动量的改变率(加速度)和作用在液体内部的压力的变化和耗散粘滞力(类似于摩擦力)以及重力之间的关系。1.2问题背景本文考虑的是这样的Couette流动：假设在两相距1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0，上平板突然以恒定的速度U=1m/s平动，求在任意时刻t两板间的速度分布。该问题的控制方程最后可总结为热传导方程：22uutx(1-8)上式是抛物型方程，因此可用有限差分法求解。并且收敛性由VonNeumann收敛性定理保证。Couette流动是流体力学中基础而又十分重要的一种流动，讨论及模拟它的方法有很多。在数值分析方面，Crank-Nicolson(克兰克-尼科尔森)格式(Crank-Nicolson差分在求解抛物线方程中广泛应用），显式古典有限差分格式，MacCormack格式[2]甚至是格子Boltzmann方法[3]等都可以很好的解决Couette流动。而在模型模拟方面，商业软件如Fluent,CFX均能很好的模拟并求解该问题。本文采用的显式古典有限差分格式，并且应用C++程序实现求解，是一种既简单有十分有效的方法，因此是有实际意义的。另一方面，对于模型方程，采用一维热传导方程，即方程(1-8)：22uutx利用VonNeumann稳定性方法可得：212txs(1-9)。这一结果是以Fourier方法，VonNeumann稳定性判据得出的。当然，也可以直接运用矩阵方法得出，得到(1-9)的方法有很多，本文只介绍以Fourier方法得到的VonNeumann稳定性。二．正文让我们回顾一下本文研究的问题：假设在两相距1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0，上平板突然以恒定的速度U=1m/s平动，求在任意时刻t两板间的速度分布。如下图所示。图一假设整个流场内部没有压力梯度，流动中运动平板对流体产生的剪力维持流体流动由于两板设为无限长，可忽略端部效应，这样可以知道每一个等x的截面速度分布相同，因此只需求x=0截面的速度分布即可。首先，我们来确定控制方程和定解条件。运用非定常二维不可压黏性流的基本方程，即：)()(022222222yuxuyuvxuutuyuxuyuvxuutuyvxu(2-1)根据此前提假设，令v=0，则由连续性方程得：0xu(2-2)。这样y方向的方程可以自动满足.而x方向的方程则可以简化为22yutu(2-3)上式其实就是(1-8)。此时，一个方程解一个未知数，方程封闭可解。近似的，取为十的负六次方。边界条件为y=0,u=0;y=1,u=U=1。而初始条件（t=0)为y1,u=0;y=1,u=1。接着，开始进行网格的划分。这是CFD中的关键步骤。简单的，采用均匀网格，网格点数设为k,即边界点分别为1和k，因此，网格空间尺度上步长为△x=1/(k-1)。为了更好的进行之后的有限差分离散方程组，先引入以下几种有限差分表达式及其模型[2]：1．X方向一阶导数的一阶向前差分1,,,ijijijuuuxx2.X方向一阶导数的一阶向后差分,1,,ijijijuuuxxx()()1,ij,ijx()(),ij1,ij3.X方向一阶导数的二阶中心差分1,1,,2ijijijuuuxx4.X方向二阶导数的二阶中心差分21,,1,22,2ijijijijuuuuxx由上面几种形式可知，对于方程22yutu(2-3)有显式古典差分格式，也就是时间导数采用向前差分，空间二阶导数采用中心差分的提法)2()()(211212111ninininininininininiuuuxtuuxuuutuu(2-4)设2)(xt。下面介绍当取何值时解是稳定的。利用有限差分格式进行计算时是按时间层逐层推进的。如果考虑x()1,ijx()1,ij()(2),ij1,ijxx()1,ijx1ioxti1it1nn1n二层差分格式，那么计算第1n层上的值1nju时，要用到第n层上计算出来的结果值nlju，nlju1，…，nlju．而计算nlju，nlju1，…，nlju时的舍入误差（包括0n的情况，不过此时是由于初始数据不精确而引起的）必然会影响到1nju的值。从而就要分析这种误差传播的情况.希望误差的影响不是越来越大，以致掩盖差分格式的解的面貌，这就是所谓稳定性的问题。对于稳定性的分析，我们一般运用Fourier方法[4]，有时候我们会用更简便的方法。在前面我们将方程(2-3)化成了便于计算的差分形式(2-4))2()(1121nininininiuuuxtuu实际上取ikjhnnjevu，代入相应的差分方程。现在代入上式中得)2()1()1(1hjiknikjhnhjiknikjhnikjhnevevevaevevikjhikhikhnikjhneeeavev)]2(1[1消去公因子ikjhe有)]2(1[1ikhikhnneeavv这时出现一个引子)2(1ikhikheea，我们把它称为是增长因子，记作),(tG．如果在实际的问题中出现矩阵，称为是增长矩阵。下面我们给出判别准则vonNeumann条件[2]：定理差分格式稳定的必要条件是当0，Tn,对所有的RK有MkGj1)),((，1j，2，…，,p(2-5)其中)),((kGj表示),(kG的特征值，M为常数．对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设ijnjnure,则错误!未找到引用源。式可写为11112ijijijijnnnnresresresre,即112iinnrrsesse放大因子1212122cos121cos14sin2niinsessesrsrssG所以141sG,为保证1G,应有212txs只要满足错误!未找到引用源。式，差分格