流体力学论文

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

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

3、流体视为没有间隙地充满它所 占据的整个空间的一种连续介质，且其所有的物理量都是空间坐标和 时间的连续函数的一种假设模型：u =u(t,x,y,z)。流体的压缩性是用体积压缩率 k来量度, dv v d ：7 k =(1-1)dp dp式中：p为外部压强。在研究流体流动过程中，若考虑到流体的压缩性，则称为可压缩 流动，相应地称流体为可压缩流体，例如高速流动的气体。若不考虑 流体的压缩性，则称为不可压缩流动，相应地称流体为不可压缩流体， 如水、油等。本文考虑的流体均为不可压缩流体。流体的粘性指在运动的状态下，流体所产生的抵抗剪切变形的性 质。粘性大小由粘度来量度。流体的粘度是由流动流体的内聚力和分

4、 子的动量交换所引起的。粘度有动力粘度 n和运动粘度y之分。动力 粘度由牛顿内摩擦定律导出：dy (1-2)式中：工为切应力，pa; r为动力粘度，pa s； du/dy为流体的剪切变 形速率。根据流体是否满足牛顿内摩擦定律，将流体分为牛顿流体和非牛 顿流体。牛顿流体严格满足牛顿内摩擦定律且 n保持为常数。非牛顿 流体的切应力与速度梯度不成正比， 一般又分为塑性流体、假塑性流 体、胀塑性流体3种。本文考虑的流体均为牛顿流体。描述流体物理量有两种方法，一种是拉格朗日描述；一种是欧拉 描述。由于本文需要只介绍欧拉描述。欧拉描述，也称空间描述，它 着眼于空间点，认为流体的物理量随空间点及时间而变化，

5、即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为(q1,q2,q3),用欧拉坐标表示的各空间点上的流体物理量如速度、压强等，在任一时刻t的值，可写为q1、q2、q3及t的函数。从数学分析知道，当某时 刻一个物理量在空间的分布一旦确定，该物理量在此空间形成一个 场。因此，欧拉描述实际上描述了一个个物理量的场。若以f表示流体的一个物理量，其欧拉描述的数学表达式是 （设 空间坐标取用直角坐标）f = f(x, y,z,t) = f(r,t)(1-3)进一步的，根据流体流动过程以及流动过程中的流体的物理参数 是否与时间相关，可将流动分为定常流动与非定常流动。 定常流动是 流体流动过程中各物理量均

6、与时间无关，这种流动称为定常流动。非定常流动是流体流动过程中某个或某些物理量与时间有关，则这种流 动称为非定常流动。本文考虑的是非定常 couette流，即其流动过程 中物理量与时间有关。流体流动所遵循的物理定律，是建立流体运动基本方程组的依 据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、 热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考 虑不同的流态，如层流与湍流。在流体力学中，系统是指某一确定流体质点集合的总体。 系统与 外界无质量交换，但可以有力的相互作用，及能量 （热和功）交换。控 制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不 动形状任意的空

7、间体积。包围这个空间体积的边界面，称为控制面。 控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的 流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能 量交换，也可与控制体外的环境有力的相互作用。下面，我们引入计算流体力学中的控制方程，我们仅介绍与本文 有关的方程。质量守恒方程（连续性方程）：质量的守恒写作:号 + v . ("v)= 0(1-4)其中p是流体的密度。在不可压缩流体的情况 v = 0p不是时间或空间的函数。方程简化为:（1-5）也就是.:u :vjw一 十 一 十(1-6):x：:y；z动量守恒方程（运动方程）：动量守恒是流体运动时应遵循的另一 个

8、普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等 于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或n-s方程，其微分形式表达如下:dup=pfbx +/p cpyx十闻zxdt次立触= pfby+也+印yy+ hzydt次y包pdw=pfbz +至hyz -t"+闻zzi dty包(1-7)式中：fbx、fby、fbz分别是单位质量流体上的质量力在三个方向上的分量；源是流体内应力张量的分量。若不采用简化书写的完整形式非 常繁琐，分别为：/ dududu 1率副?馈-川卜晦-劫+余(dv a,加。 dp d l%2 r_ al a r始丫上a /加上加

9、丫上组就是著名的n-s方程，以克劳德-路易纳维和乔治加布 里埃尔斯托克斯命名，这些方程建立了流体的粒子动量的改变率（加 速度）和作用在液体内部的压力的变化和耗散粘滞力（类似于摩擦力） 以及重力之间的关系。1.2问题背景本文考虑的是这样的couette流动：假设在两相距1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0,上平板突然以恒定的速度u=1m/s平动，求在任意时刻t两板间的速度分布该问题的控制方程最后可总结为热传导方程:(1-8)-u：2u-7 -2t x上式是抛物型方程，因此可用有限差分法求解。并且收敛性由 von neumann攵敛性定理保证。couette 流动是流

10、体力学中基础而又十分重要的一种流动，讨论 及模拟它的方法有很多。在数值分析方面，crank-nicolson（克兰克- 尼科尔森）格式（crank-nicolson 差分在求解抛物线方程中广泛应 用），显式古典有限差分格式，maccormack格式2甚至是格子boltzmann方法3等都可以很好的解决 couette流动。而在模型模 拟方面，商业软件如fluent,cfx均能很好的模拟并求解该问题。本 文采用的显式古典有限差分格式，并且应用c+播序实现求解，是一种既简单有十分有效的方法，因此是有实际意义的。另一方面，对于模型方程，采用一维热传导方程，即方程（1-8）:2 u二 u=ot -x-

11、 2t xt 1 s 二o -利用von neumann!急te性方法可得：取 2 （1-9）。这一结果是以fourier方法，von neumann稳定性判据得出的。 当然，也可以直接运用矩阵方法得出，得到 （1-9）的方法有很多，本 文只介绍以fourier方法得到的von neumann!急定性。二.正文让我们回顾一下本文研究的问题：假设在两相距 1m的无限大平板间充满水，平板原来都处于静止状态，在某一时刻t=0,上平板突然以恒定的速度u=1m/s平动，求在任意时刻t两板间的速度分布。 如下图所示。(2-1)假设整个流场内部没有压力梯度，流动中运动平板对流体产生的 剪力维持流体流动由于两

12、板设为无限长， 可忽略端部效应，这样可以 知道每一个等x的截面速度分布相同，因此只需求x=0截面的速度分 布即可。首先，我们来确定控制方程和定解条件。运用非定常二维不可压十-2二u-2二 y黏性流的基本方程，即:个nn二u二u二 u u v ：t：x：y口l二u二u二u u v t-x-y根据此前提假设，令v=0,则由连续性方程得：ix = 0(2-2)。这样y方向的方程可以自动满足.而x方向的方程则可以简化为一一 2u u-t-=-2t y (2-3)上式其实就是(1-8)。此时，一个方程解一个未知数，方程封闭可 解。近似的，取n为十的负六次方。边界条件为 y=0,u=0;y=1,u=u=1

13、。而初始条件(t=0)为 y<1,u=0;y=1,u=1。接着，开始进行网格的划分。这是 cfd中的关键步骤。简单的， 采用均匀网格，网格点数设为k,即边界点分别为1和k,因此，网格 空间尺度上步长为 x=1/(k-1)。为了更好的进行之后的有限差分离散方程组，先引入以下几种有限差分表达式及其模型2：1. x方向一阶导数的一阶向前差分i , ji +1, j2.x方向一阶导数的一阶向后差分包ui,j -ui4,ji改,j&x3.x方向一阶导数的二阶中心差分：u _ ui 1,j - ui,j i,j2 x4. x方向二阶导数的二阶中心差分2 .2 y,u uui书,j 2ui,

14、j + ui,jtt =21出八j(ax)由上面几种形式可知，对于方程一.2u 1 二 u=r . 2t y （2-3）导数采用中心差分的提法n 1 nui-ui二 tn 1 nuiuinn ni ui 1 -2ui ui4=j2:(汉)日5 -2陋5)(x)(2-4)有显式古典差分格式，也就是时间导数采用向前差分，空间二阶设儿=鼻。下面介绍当取何值时解是稳定的（m利用有限差分格式进行计算时是按时间层逐层推进的如果考虑n ! 1二层差分格式，那么计算第n+1层上的值uj时，要用到第n层上计算 nnnnnn出来的结果值uj，uja,，uj+.而计算uj，uj中，uj卡时的舍入误差(包括n=0的情

15、况，不过此时是由于初始数据不精确而引起n : 1的)必然会影响到uj的值。从而就要分析这种误差传播的情况.希望 误差的影响不是越来越大，以致掩盖差分格式的解的面貌，这就是所 谓稳定性的问题。对于稳定性的分析，我们一般运用 fourier方法4,有时候我们 会用更简便的方法。在前面我们将方程(2-3)化成了便于计算的差分形 式(2-4)n 1 n i-t z n n n .ui =ui ,2(ui 1 -2ui - ui j)(以)n n ikjh实际上取uj =v e ，代入相应的差分方程。现在代入上式中得n t ikjh n ikjhn ik( j41)hn ikjh n ik( j)hv

16、e =ve a(ve - 2v e v e )vn1eikjh =vn1 a，(eikh -2 e，eikjh消去公因子eikjh有n 1 nikh . ikhv v 1 a , (e e -2) ikh-ikh这时出现一个引子化ame +e 一2),我们把它称为是增长因子，记作g(7).如果在实际的问题中出现矩阵，称为是增长矩阵。下面我们给出判别准则von neumann条件2：定理差分格式稳定的必要条件是当飞0, n丁 ,对所有的k乏r有九j(gg,k) w1+mwj=12 . p, (2-5)其中上j(gck)表示g(&k)的特征值，m为常数.对于适定的线性微分 方程，格式如果差

17、分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了n i jr设uj=rnel |则错误！未找到引用源。式可写为rn川叽srnei"p+(1-2s)rnei(巧+$赭9，即rn 1. =rn sei? 1 -2s se，丁放大因子g = ” =st 1 -2s se 巧 rn=1 -2s 2scos?=1 -2s 1 -cos/. 2 v=1 -4ssin 一2所以1-公必为保证仃工1,应有二 t . 1s 二2 -x 2只要满足错误！未找到引用源。式，差分格式就是稳定的。因此，将上述结论用于方程 (2-3)可知当九=1/2时，差分稳定。接下来，考虑程序的编写与调试。程序中

18、网格节点数k作为变量, 用u(i,1)表示u：,用u(i,2)表示由uin+ = uin + n :t 2 (uin+ - 2un + ul)(x)(2-4)计算出n+1时间层上的节点函数，然后赋值给u(i,1)数组，完成对n 时间层上节点函数的更新，再重复循环操作，不断推进直到到达tmax。我们可以通过tmax取不同的值来提高结果的正确度。这种迭代计算 是编程中常见的手段。在进入代码之前我们先补充边界条件的离散形式：nnu 1 =0,uk =1,n =1,2,3 初始条件的离散形式：u1 =0,i =1,2,3k1 ; u； =1。完成了以上分析之后，列出代码如下：/*本段代码用c+语言编写

19、5,并在visual stadio2013上运行。为了输入方便，用v代替n , u代替九*/#include "math.h"#include <iostream> 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

20、/ (float)(k - 1);for (i = 1; i <= k; i+)yi = dy*(i - 1);/*dt = u*dy*dy / v;tmax = 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

21、 + (float)(u*(ui + 11 - 2.0*ui1 + ui - 11);for (i = 1; i <= k; i+)ui1 = ui2;n+; while (n < nmax);/*输出结果*/for (i = 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为网格节点

22、数，在本次运行中我们一律取为21,也就是y方向网格数为20。而nmax得值于tmax密切相关，即给定了时间 推进的截止时间，我们下面做了三组以 nmax为控制变量的运行，分 别取nmax=20, 200, 2000。而 v,也就是r , 一律取为一律取为 10'm/s2。最后u,也就是九，我们将以其取值来讨论抛物线方程的稳 定性定理。首先，第一组试验取值取 k=21, nmax=20,v=10m/s2 ,u=0.25。按此输出后得到如下数据。进而在excelt可直观画出此时的速度分布图像,如图一所不y(m)u(m/s)000.0500.100.1500.200.2500.30.0000

23、020.350.0000140.40.000070.450.000294nmax=20, tamx=208min0.50.0010650.550.0033780.60.0094750.650.0237030.70.0532520.750.1081290.80.1995910.850.3367840.90.5223970.950.74925911第二组试验取值取k=21, nmax=200,v=i0-6m/s2 ,u=0.25o按此输出后得到如下数据。进而在 excelt可直观画出此时的速度分布图像,如图二所不y(m)u(m/s)000.050.0216190.10.0438850.150.06

24、74360.20.0928850.250.1208130.30.1517550.350.1861890.40.2245220.450.2670760.50.3140930.550.3656670.60.4218410.650.48250.70.5474180.750.6162520.80.6885470.850.7637460.90.8412050.950.9202111i >hmax=2001 tmax=2083min1n ajtb a a/n trn ?j7u (0020.40.160.r611.2图三第三组试验取值取k=21, nmax=2000,v=io"m/s2 ,u

25、=0.25。按此输出后得到如下数据。进而在excel±可直观画出此时的速度分布图像，如图四所示。y(m)u(m/s)00runax=2000j tmak =208330.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.649996图四0.70.6999960.750.7499970.80.7999970.850.8499980.90.8999980.950.94

26、999911在上面的运行中，我们完成了三组输入值不同的操作，进一步的,将三张图集中到一张图来，如图五nmax=20时，不同的tmax由图五不难看出，1m射间为年制变量,随着时间的增加，上平板 拖动的影响区将逐渐的增大，当经过很长时间(如图中的 tmax=20833min)时，速度u沿y方向上的分布将趋近于线性。上三次操作最后都得到了收敛的稳定的解。那么，什么时候无法得到方程(2-3)稳定解呢？我们继续进行如下试验。第四组试验取值取 k=21, nmax=2000,v=io'm/s2 ,u=0.51。此时运行程序，运行结果如下图六/luaxri隼piiikfct 憎机-li, u 强口m

27、u值力般峋麻舒.a芸0日期 uibi wjcoi nkim' rail>«1j - dtij lel_wliie2 i-hi. mmhhhhvi. u ix h2 i-19¥mz!i4-44r42l4z4¥ i：轴l 篦.uijh21-3b71 黯 qud j31mqi 旗榄*.mrhh“。口 tl. u ” bia i吗右白鹏嗝也用”的版e摊总2wornsuisna 1咏hbaf16hf.uwawi,iina i rrsbe蒙曲"knnfwra ffmtmt 7i7h1-3wui4h.uiii2 i- i m i14hmuhhniahi.

28、jbwiw.uiriii iiiibhlmzeifimutih涮_mwhhh “51 1 总. h配.u ” h21 719：1 田2<is却q"的e® _魄ihhh “:i b蚓dfl. u* i«h g斗阚越两叮之b1的整4.蚂as砸 f m p-fl. uhilcz -1eismw6s54mme.»»»«»«t 1 1 ? pt ssrhhhj ui li lies -is?th-watwai?mnsi. hhhhhh 7i11 i>«i-i'.|amkmh .uhi

29、ie2 i" :h t1 45272-4：l!i274tt<iwv»». wmiihhfl fl1* 口. uu i i i ln j-11 kelih'/zimm 44jh4«iw_mwuhhh hit:4一?豌&hh.uls i12】=一i白1西igqdseq?士»»hhhh hie i-s-tseeaa仙i« i ca】一m“际谷日酒4晒 g 附缸国邺的 口 ？ t. uh71c2 j “7城耳门41 两£65葭皓g<u |i if jehhhnvm*)c3 eh珞hikwlg&q

30、uot;?，mamm 11 !f ifjiwhhh .uiit he3l»4lb：l：iim4d tkmu.hmmuhhflxh 1-41-hhhmn .u92hil2 j-ltf!jtia44421443h«mh.»mmi4hh f23 l-i.hem ,u1±i je2 -1.的*g ifwils-iffittei. -.再取 wproji&<l4d图六如图六所示，此时得到的解不再收敛，而是没有规律的发散了,说明此时u=0.51的设置使方程（2-3）无法得到准确定解再设置试验五，试验六分别为k=21, nmax=2000,v=i0-6m

31、/s2 ,u=0.5和k=21, nmax=2000,v=io'm/s2 ,u=0.49。给出两组试验运行截图，分别为图七，图八/tuaxii|1 划幻 |专pbd>cm-等廿计|（上国1«1£“21】3“打一日.nc2 3-b 卜口】一日k!£13 j-mh141th15 3-bpiuhfcrlstl-m, o m 田-mm u1bajwkwujcijwhmirli i i z iiz k4m9nlpu 114«ml.u14lt±-a.fi45i»n畏盹4e. u heu1btwhvhjjei? 11 n嗝股 nnmr

32、iihibixi4.i4iivftkmm.u 2b lri±j-a.95bfl« bwm«, u£iin±i-i. enseei h.sn»u>l na hh4mu-u 1 1 11 t_wu4mlm 的询虬u ra1±。.函就i «e#eou3ti-«.ie«999ru b 111 »wnariiil iiz i-2-hvvvv 3mms_u 7 fe_29994# s«e«.u6£2k-e.34ff#»«m»n. u e¥ kl-fl.39¥¥9ixigw图七血日*5，5 ：131.王2