流体力学论文

1、 C+程序实现非定常COUETTE流速度模拟及古典显式有限差分格式收敛性验证院系名称能源与动力工程学院专业名称飞行器动力工程作者姓名洪启臻(13041222)2015 年 5月 13日C+程序实现非定常COUETTE流速度模拟及古典显式有限差分格式收敛性验证摘要 非定常Couette流（库艾特流）是指两平行平板间的非定常流，是流体力学中简单而又基础的一种流动方式。有限差分法是CFD(计算流体力学）中常用的方法，其中本文考虑的是扩散式抛物线方程，即热传导方程，运用古典显式有限差分格式，也就是时间导数采用向前差分，空间二阶导数采用中心差分来研究Couette流。进而通过C+编程来进行Y方向上的速度

2、分布模拟，通过EXCEL做成表格。分别考察了t等于20,200,2000的情况，得出相应结果。最后，就古典显式有限差分格式的收敛性给出了讨论，当系数u大于0.5时，求解将不再收敛，进而验证了数值分析中Von Neumann得出的结论。关键字：Couette流（库艾特流） C+ 古典显式有限差分 收敛性1 引言1.1基础知识1 下面，我们先引入流体力学的一些基本假设及控制方程。 我们认为，流体质点是几何尺寸同流动空间相比是极小量，又含有大量分子的微元体。连续介质则是质点连续地充满所占空间的流体或固体。进而提出连续介质模型，即把流体视为没有间隙地充满它所占据的整个空间的一种连续介质，且其所有的物理

3、量都是空间坐标和时间的连续函数的一种假设模型：u =u(t,x,y,z)。 流体的压缩性是用体积压缩率k来量度 (1-1)式中：p为外部压强。 在研究流体流动过程中，若考虑到流体的压缩性，则称为可压缩流动，相应地称流体为可压缩流体，例如高速流动的气体。若不考虑流体的压缩性，则称为不可压缩流动，相应地称流体为不可压缩流体，如水、油等。本文考虑的流体均为不可压缩流体。 流体的粘性指在运动的状态下，流体所产生的抵抗剪切变形的性质。粘性大小由粘度来量度。流体的粘度是由流动流体的内聚力和分子的动量交换所引起的。粘度有动力粘度和运动粘度之分。动力粘度由牛顿内摩擦定律导出： (1-2)式中：为切应力，Pa；

4、为动力粘度，Pas；为流体的剪切变形速率。 根据流体是否满足牛顿内摩擦定律，将流体分为牛顿流体和非牛顿流体。牛顿流体严格满足牛顿内摩擦定律且保持为常数。非牛顿流体的切应力与速度梯度不成正比，一般又分为塑性流体、假塑性流体、胀塑性流体3种。本文考虑的流体均为牛顿流体。 描述流体物理量有两种方法，一种是拉格朗日描述；一种是欧拉描述。由于本文需要只介绍欧拉描述。欧拉描述，也称空间描述，它着眼于空间点，认为流体的物理量随空间点及时间而变化，即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为(q1,q2,q3)，用欧拉坐标表示的各空间点上的流体物理量如速度、压强等，在任一时刻t的值，可写为q1、q2

5、、q3及t的函数。从数学分析知道，当某时刻一个物理量在空间的分布一旦确定，该物理量在此空间形成一个场。因此，欧拉描述实际上描述了一个个物理量的场。 若以f表示流体的一个物理量，其欧拉描述的数学表达式是(设空间坐标取用直角坐标) (1-3) 进一步的，根据流体流动过程以及流动过程中的流体的物理参数是否与时间相关，可将流动分为定常流动与非定常流动。定常流动是流体流动过程中各物理量均与时间无关，这种流动称为定常流动。非定常流动是流体流动过程中某个或某些物理量与时间有关，则这种流动称为非定常流动。本文考虑的是非定常Couette流，即其流动过程中物理量与时间有关。 流体流动所遵循的物理定律，是建立流体

6、运动基本方程组的依据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考虑不同的流态，如层流与湍流。 在流体力学中，系统是指某一确定流体质点集合的总体。系统与外界无质量交换，但可以有力的相互作用，及能量(热和功)交换。控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动形状任意的空间体积。包围这个空间体积的边界面，称为控制面。控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能量交换，也可与控制体外的环境有力的相互作用。下面，我们引入计算流体力学中

7、的控制方程，我们仅介绍与本文有关的方程。 质量守恒方程(连续性方程)：质量的守恒写作： (1-4)其中是流体的密度。 在不可压缩流体的情况 不是时间或空间的函数。方程简化为： (1-5)也就是 (1-6) 动量守恒方程(运动方程)：动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或N-S方程，其微分形式表达如下： (1-7)式中：、分别是单位质量流体上的质量力在三个方向上的分量；是流体内应力张量的分量。若不采用简化书写的完整形式非常繁琐，分别为： 上组就是著名的N-S方程，以克劳德

8、-路易纳维和乔治加布里埃尔斯托克斯命名，这些方程建立了流体的粒子动量的改变率(加速度)和作用在液体内部的压力的变化和耗散粘滞力(类似于摩擦力)以及重力之间的关系。1.2问题背景 本文考虑的是这样的Couette流动：假设在两相距1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0，上平板突然以恒定的速度U=1m/s平动，求在任意时刻t两板间的速度分布。 该问题的控制方程最后可总结为热传导方程： (1-8) 上式是抛物型方程，因此可用有限差分法求解。并且收敛性由Von Neumann收敛性定理保证。 Couette流动是流体力学中基础而又十分重要的一种流动，讨论及模拟它的方法有很多

9、。在数值分析方面，Crank-Nicolson(克兰克-尼科尔森)格式(Crank-Nicolson差分在求解抛物线方程中广泛应用），显式古典有限差分格式，MacCormack格式2甚至是格子Boltzmann方法3等都可以很好的解决Couette流动。而在模型模拟方面，商业软件如Fluent,CFX均能很好的模拟并求解该问题。本文采用的显式古典有限差分格式，并且应用C+程序实现求解，是一种既简单有十分有效的方法，因此是有实际意义的。 另一方面，对于模型方程，采用一维热传导方程，即方程 (1-8)： 利用Von Neumann稳定性方法可得： (1-9)。 这一结果是以Fourier方法，Vo

10、n Neumann稳定性判据得出的。当然，也可以直接运用矩阵方法得出，得到(1-9)的方法有很多，本文只介绍以Fourier方法得到的Von Neumann稳定性。2 正文 让我们回顾一下本文研究的问题：假设在两相距1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0，上平板突然以恒定的速度U=1m/s平动，求在任意时刻t两板间的速度分布。如下图所示。图一 假设整个流场内部没有压力梯度，流动中运动平板对流体产生的剪力维持流体流动由于两板设为无限长，可忽略端部效应，这样可以知道每一个等x的截面速度分布相同，因此只需求x=0截面的速度分布即可。首先，我们来确定控制方程和定解条件。运用

11、非定常二维不可压黏性流的基本方程，即： (2-1) 根据此前提假设，令v=0，则由连续性方程得：(2-2) 。这样y方向的方程可以自动满足.而x方向的方程则可以简化为 (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)。为了更好的进行之后的有限差分离散方程组，先引入以下几种有限差分表达式及其模型

12、2：1X方向一阶导数的一阶向前差分 2.X方向一阶导数的一阶向后差分 3.X方向一阶导数的二阶中心差分 4. X方向二阶导数的二阶中心差分 由上面几种形式可知，对于方程 (2-3) 有显式古典差分格式，也就是时间导数采用向前差分，空间二阶导数采用中心差分的提法(2-4)设。下面介绍当取何值时解是稳定的。 利用有限差分格式进行计算时是按时间层逐层推进的。如果考虑二层差分格式，那么计算第层上的值时，要用到第层上计算出来的结果值，而计算，时的舍入误差（包括的情况，不过此时是由于初始数据不精确而引起的）必然会影响到的值。从而就要分析这种误差传播的情况.希望误差的影响不是越来越大，以致掩盖差分格式的解的

13、面貌，这就是所谓稳定性的问题。对于稳定性的分析，我们一般运用Fourier方法4，有时候我们会用更简便的方法。在前面我们将方程(2-3)化成了便于计算的差分形式(2-4)实际上取，代入相应的差分方程。现在代入上式中得 消去公因子有 这时出现一个引子，我们把它称为是增长因子，记作如果在实际的问题中出现矩阵，称为是增长矩阵。 下面我们给出判别准则von Neumann条件2：定理 差分格式稳定的必要条件是当，,对所有的有 ， ， (2-5)其中表示的特征值，为常数对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。 设,则式可写为 ,即 放大因子 所以

14、,为保证,应有 只要满足式，差分格式就是稳定的。因此，将上述结论用于方程 (2-3)可知当=1/2时，差分稳定。接下来，考虑程序的编写与调试。程序中网格节点数k作为变量，用U(i,1)表示，用U(i,2)表示，由 (2-4)计算出n+1时间层上的节点函数，然后赋值给U(i,1)数组，完成对n时间层上节点函数的更新，再重复循环操作，不断推进直到到达tmax。我们可以通过tmax取不同的值来提高结果的正确度。这种迭代计算是编程中常见的手段。在进入代码之前我们先补充边界条件的离散形式：初始条件的离散形式：；。完成了以上分析之后，列出代码如下：/*本段代码用C+语言编写5，并在Visual Stadi

15、o2013上运行。为了输入方便，用v代替，u代替*/#include math.h#include using namespace std;int main()float Y200, U2003;int i, j, k, nmax;double v, dt, tmax;float u, dy, n;scanf_s(%d%d%lf%f, &k, &nmax, &v, &u);/*GRID GENNERATION网格生成*/dy = 1 / (float)(k - 1);for (i = 1; i = k; i+)Yi = dy*(i - 1);/*dt = u*dy*dy / v;tmax =

16、nmax*dt;INITIAL CONDITIONS初始条件*/for (i = 1; i = k - 1; i+)Ui1 = 0;Uk1 = 1;n = 1;do/*BOUNDARY CONDITION边界条件*/for (j = 1; j = 2; j+)U1j = 0.0;Ukj = 1.0;/*差分计算*/for (i = 2; i = k- 1; i+)Ui2 = Ui1 + (float)(u*(Ui + 11 - 2.0*Ui1 + Ui - 11);for (i = 1; i = k; i+)Ui1 = Ui2;n+; while (n nmax);/*输出结果*/for (i

17、 = 1; i = k; i+)printf(Y%d=%f,U%d2=%fn, i, Yi, i, Ui2);system(pause); 该程序有4个输入值，即语句scanf_s(%d%d%lf%f, &k, &nmax, &v, &u); k为网格节点数，在本次运行中我们一律取为21，也就是y方向网格数为20。而nmax得值于tmax密切相关，即给定了时间推进的截止时间，我们下面做了三组以nmax为控制变量的运行，分别取nmax=20，200, 2000。而v,也就是，一律取为一律取为。最后u，也就是，我们将以其取值来讨论抛物线方程的稳定性定理。 首先，第一组试验取值取k=21,nmax=

18、20,v=,u=0.25。按此输出后得到如下数据。进而在EXCEL上可直观画出此时的速度分布图像，如图二所示。Y(m)U(m/s)000.0500.100.1500.200.2500.30.0000020.350.0000140.40.000070.450.0002940.50.0010650.550.0033780.6图二0.0094750.650.0237030.70.0532520.750.1081290.80.1995910.850.3367840.90.5223970.950.74925911第二组试验取值取k=21,nmax=200,v=,u=0.25。按此输出后得到如下数据。进而

19、在EXCEL上可直观画出此时的速度分布图像，如图三所示。Y(m)U(m/s)000.050.0216190.10.0438850.150.0674360.20.0928850.250.1208130.30.1517550.350.1861890.40.2245220.450.2670760.50.3140930.550.3656670.60.4218410.650.48250.70.5474180.75 图三0.6162520.80.6885470.850.7637460.90.8412050.950.9202111第三组试验取值取k=21,nmax=2000,v=,u=0.25。按此输出后得

20、到如下数据。进而在EXCEL上可直观画出此时的速度分布图像，如图四所示。Y(m)U(m/s)000.050.050.10.099990.150.1499990.20.1999980.250.2499980.30.2999970.350.3499970.40.3999970.450.4499960.50.4999960.550.5499960.60.5999960.650.6499960.7图四0.6999960.750.7499970.80.7999970.850.8499980.90.8999980.950.94999911 在上面的运行中，我们完成了三组输入值不同的操作，进一步的，将三张图

21、集中到一张图来，如图五。图五由图五不难看出，以时间为控制变量，随着时间的增加，上平板拖动的影响区将逐渐的增大，当经过很长时间（如图中的tmax=20833min)时，速度U沿y方向上的分布将趋近于线性。上三次操作最后都得到了收敛的稳定的解。那么，什么时候无法得到方程(2-3)稳定解呢？我们继续进行如下试验。第四组试验取值取k=21,nmax=2000,v=,u=0.51。此时运行程序，运行结果如下图六。图六如图六所示，此时得到的解不再收敛，而是没有规律的发散了，说明此时u=0.51的设置使方程(2-3)无法得到准确定解。再设置试验五，试验六分别为k=21,nmax=2000,v=,u=0.5和

22、k=21,nmax=2000,v=,u=0.49。给出两组试验运行截图,分别为图七，图八。图七图八 不难看出，试验五，六得出的结果与试验三一致。将其与试验四进行对比，不难看出，使得方程稳定的条件是u=1/2,也就是=1/2。说明稳定性与网格的空间间距，时间间距有关，证明了von Neumann条件的实际正确性。3 总结 通过前面的探讨，C+编写的程序可以较好地进行Y方向上的速度分布模拟，得出的结果通过EXCEL做成表格与实际情况相符合。方程 (2-3)是扩散式抛物线方程，运用古典显式有限差分格式，也就是时间导数采用向前差分，空间二阶导数采用中心差分来研究时，其具有一阶时间精度和二阶空间精度，是符合计算要求的。在文中所