库埃特流

不可压Couette流的数值模拟姓名：安彬学号：201101003032专业：空间工程1、问题的提出库埃特（Couette）流动定义如下：设有两个相距为D的平行平板，上面的平板以速度,eu运动，下面的平板静止，速度0u。考虑这两个平板之间的粘性流动。在xy平面内，流动图1所示。两平板间产生流动的驱动力只有一项：由上平板运动引起的、作用于流体上的切应力。由此产生了横截面上的速度剖面()uuy。图1这种流动的控制方程是x方向的动量方程对于库埃物流动，这个方程还可以大大的简化。为此，考察图1，我们注意到库埃特的模型在x轴正、负方向上都无限延伸。既然这种流yxxxzxxDupDtxxyz动没有起点和终点，那么流场的变化必定与x无关，即，所有量的0x。对于连续性方程0tV，将其应用于定常流动，有0uxy1既然库埃特流动中0ux，那么方程1变为0yyy2在下平板即0y处取值，有0，则方程2化为00yy3如果我们在0y点将展开为泰勒级数，就有2220002yyyyyyy4在上平板处取值，式4成为2220002yyDDDyy5由于0D，00，而且由方程3有0yy＝0，那么惟一合理的结果只能是：对于所有的n，00nnyy，由此可知在整个流场中06这是库埃特流动的物理特征。也就是说，流场中任何一点的垂直速度分量都为零。这表明库埃特流线是平行的直线，这一结果在观察图1时凭直觉就可以想到。最后，将y方向动量方程xyyyzyyDpDtyxyz应用于没有体积力的库埃特流动，就有0yypyy7由2yyVy得20yyuvxyy8因此，方程7简化为0py9于是可以得出这样的结论：对于库埃特流动，在x方向和y方向上都没有压力梯度。此时让我们回到x方向的动量方程。对于无体积力的定常二维流动，我们有yxxxuuuxyxxy10根据方程2xxuVx和xyyxuxy，对于库埃特流动20xxuvuxyx11yxuuuxyy12将式11和式12代入方程10，得0uyy13如果假设流动是不可压的、恒温的，则常数。于是13变为220uy14方程14就是不可压恒温库埃特流动的控制方程。求方程14的解析解很简单。对y进行两次积分，我们有12ucyc15式中，1c和2c是积分常数，它们的值由边界条件确定。因为在下平板0y处，0u（图1），从式15可推出20c。在上平板yD处，已知euu。于是在式15中1/ecuD。代入这些值，式15变为euyuD16式16是不可压库埃特流动速度分布的解析解。从中我们注意到，解析解u仅随变化，而且是线性的。图1画出了这种线性的速度分布。2、求解方法简介在数值求解中，我们这样来设置：假定速度剖面不同于式16给出的解析解，而不是线性。例如，假设速度剖面为0,对0yD；u=u,当euyD17将它作为计算的初始速度剖面，即0t时刻的初始条件。我们从这个初始条件出发，建立流场的时间推进解法。在推进过程中，速度剖面的形状会逐渐趋于一条直线。当速度剖面的形状趋于直线时，我们就认为此已达到定常状态，所以在数值计算中，我们使用的是用非定常库埃特流来逼近定常库埃特流，最终得到其稳定的速度分布。非定常的库埃特流的控制方程易得为22uuty18方程18是一个抛物线型偏微分方程，因此时间推进的结果代表了一个适定的解。3、差分格式的建立首先把方程18表示为无量纲的形式。定义如下无量纲变量'euuu'yyD'/ettDu把方程18无量纲化'2'''2euuuDty19在方程19中我们发现1eeDuDR其中eDR是按照两平板之间距离D计算的雷诺数，于是方程19变为'2'''21eDuuRty20式20就是要用数值求解的方程。为了得到数值解，我们分别应用显式和隐式的有限差分格式，显式格式使用时间前差空间中心差（FTCS）格式，隐式格式则使用克拉克－尼科尔森（Crank-Nicolson）格式。为了简化记号，我们省略方程20中所有变量上的“’”号，得到221eDuutRy21在FCFS格式中，用向前差分代替了ut，用二阶导数的二队精度中心差分格式代替了22uy，得到了方程21这种特定形式的差分方程，即111221nnnnniiiiieDuuuuutRy22整理后可写成11122nnnnniiiiieDtuuuuuRy23方程23即为方程21的显式差分格式。在Crank-Nicolson格式中，将方程22中右边的空间差分表示成第n个时间层上的量与第n+1个时间层上的量的平均值，就是用11111111222221nnnnnniiiiiinniieDuuuuuuuutRy或111111112222nnnnnnnniiiiiiiieDtuuuuuuuuyR24来表示方程21，把等式24中所有n+1层的项移至等式左边并整理，方程24改写为1111122211221221()2nnniiieDeDeDnnniiieDeDtttuuuyRyRyRttuuuyRyR25方程25可表示为如下形式11111nnniiiiAuBuAuK26其中22eDtAyR2721eDtByR28112212nnniiiieDeDttKuuuyRyR29求解方程26的网格如图2所示。图中，平板间的垂直距离D（y方向）被N等分（共N+1个网格点），每段长度为y，即DyN30图21u和1Nu可以通过边界条件得到10u3111Nu32式26代表了N－1个未知数2u，3u，，Nu形成的N－1个方程，形成一个方程组。将这个方程组具体写出来，其第一个方程是1111232nnnAuBuAuK33因为10u，因此方程33变为N+1111232nnBuAuK34式26所代表的最后一个方程为11111nnnNNNNAuBuAuK35因为11Nu。因此方程35变为111nnNNNAuBuKA37于是，由式26所代表的方程组可以用矩阵形式表示成122133144155111100000000000000000000000000000000000000nnnnnNNnNNBAKuABAKuABAKuABAKuKuABAKAuAB38解这个方程组，将得到解12nu，13nu，，1nNu。它们是n+1时刻的速度值。反复进行这个过程，直到速度剖面收敛到一个稳定的状态。4、差分格式的性质分析下面分析显式格式和隐式格式的性质。显式格式的建立相对比较简单，其格式也使编程实现较为容易。根据冯·诺伊曼（VonNeumann）稳定性判据得，212eDtRy。一般情况下，对于网格步长y是一定的，所以就要求，在数值求解过程中，t不能过大。隐式格式的建立与求解虽然较为复杂，每一步的求解也需花费更多的时间，但是在算法稳定性上则远远优于显式格式，有些隐式格式是无条件稳定的，即对于任意的t都能得到最后的定常解，这就能大大减小计算的次数，从而提高计算速度。但若t无限扩大，则会使得每一步的计算精度得不到满足。5、程序求解对于显式格式，根据VonNeumann稳定性判据得到212eDtRy,根据式23则可以求出每个时间步骤的整个速度剖面的所有速度，使用MABLAB画出其图形。对于隐式格式，式24为无条件稳定，在讨论过程中定义2eDtcRy，然后对38进行求解，那可以得到每个时间步骤的整个速度剖面的所有速度，同样使用MATLAB画出其图形。改变c或dt观察图形及循环次数的变化。6、结果分析应用显式格式求解，得到如果所示图形。Delta=0.3R=875Delta=0.4R=625Delta=0.4R=521Delta=0.6结果发散从图中容易得到，速度分布最终趋向于一条直线，即速度分布为正比例分布。与式16一致，因此我们可以说这种方法是可行的。当Delta增大时达到相同精度的循环次数减少，当Delta大于0.5时因出现发散而无法得到数值解。应用隐式格式求解，得到如下图形dt=0.2c=0.5R=16dt=0.1c=1R=32dt=0.05c=2R=65dt=0.05c=1R=131dt=0.1c=3R=23容易得到，使用隐式格式求解得到的解与显式格式得到的解是一致的，也与式16得到的结果一致。观察图像可知，从总体上讲达到相同的精度隐式格式所用的循环数大大减少。循环次数随着dt、c的减小而增大，dt和c越大循环次数越少，但当c和dt过大时计算过程中会造成严重截断误差和舍入误差。因此合理选择dt和c的值对于实现高效运算十分重要。附录显式格式的MATLAB源程序clcclearalldelta=0.4;u(1:21)=0;u(22)=1;y=[0:0.05:1];he=[01];zon=[01];plot(he,zon,'r')holdallc=[0:0.05:1];i=0;while10forn=2:1:21u2(n)=delta*(u(n+1)-2*u(n)+u(n-1))+u(n);endu2(22)=1;u2(1)=0;u2(2)=0;u21=u2(2:22);if(max(abs(u21-c))0.001)break;endu=u2;u3=u2(2:22);plot(y,u3,'b-')holdalli=i+1;endxlabel('y')ylabel('u')稳式格式的MATLAB源程序clc;clearall;dt=0.2;dy=0.2;r=5000;%c=dt/(2*dy*dy*r);c=0.6;a=-c;b=1+2*c;u(1:(1/dt))=0;u(1/dt+1)=1;A(1:(1/dt-2))=a;B(1:(1/dt-1))=b;x0=diag(B,0);x1=diag(A,1);x2=diag(A,-1);x=x0+x1+x2;y=[0:dy:1];i=0;while10forn=2:1:1/dtk(n)=(1-2*c)*u(n)+c*(u(n+1)+u(n-1));endk(n)=k(n)-a;k1=k(2:n);u(2:1/dt)=k1/x;u(1)=0;u(1/dt+1)=1;plot(y,u,'b-');holdall;i