1、-E.1-附录 E 二维不可压缩黏性流体 流动问题Couet的有限体积算法与计算程序二维 流动是一个不可压缩黏性流动的典型流动，并有解析解，可以Couet用来检验数值算法计算精度和可靠性。对它采用有限体积算法一阶迎风型离散格式进行数值求解。同时，为了初学者入门和练习方便,这里给出了用 语言和C语言编写的计算二维不可压缩黏性流体 流动问题计算程序，供Fortan7 Couet大家学习参考。E-1 利用有限体积算法一阶迎风型离散格式求解二维不可压缩黏性流体 流动问题Couet1.二维不可压缩黏性流体 流动问题的提法二维不可压缩黏性流体 流动：有两个无限长平板，间距为 ，两板t L之间充满密度为 1
2、、静止的不可压缩黏性流体。无限长平板组成的通道两端压力相等，下板固定不动，上板以量纲为一的速度 1 自左向右平移运动（图 E.1） 。2. 基本方程组、初始条件和边界条件设流体是黏性流体。二维 流动问题在数学上可以由二维不可压缩黏Couet性流动 方程组来表示，把它写成通用变量的微分方程组形式，有：N-S（E.1）vSxyxy其中 为变量 在水平 方向的流速， 为 在垂直 方向的流速， 为黏度，u 为源项。源项中不仅包含压力梯度项，也包含时间导数项。S图 E.1 二维不可压缩黏性流体 流动问题示意图ouet-E.2-初始条件：上板以量纲为一的速度 1 沿着上壁面方向自左向右运动。边界条件：流动
3、速度 在上下壁面可采用无滑移边界条件，在左右两端uv、采用自由输出边界条件；压强 采用自由输出边界条件。p3. 计算网格划分和控制体单元与节点定义采用交错网格。图 E.2 和图 E.3 是计算网格、控制体单元和节点示意图。节点 所在主控制体单元如图 E.2 中有阴影部分所示。在 方向与节点P x相邻的节点为 和 ，在 方向与节点 相邻的节点为 和 ，主控制单元WEyPSN界面分别为 。压力 和速度 分别在三套不同网格中，如图 E.3 中有wsen、 、 、 puv、阴影部分所示。4. 有限体积算法离散格式对方程（E.1）在图 E.2 所示节点 主控制单元内积分，有：P（E.2VVVVVudvd
4、ddSxyxy）由于不可压缩黏性流体 流动是二维问题，因此，控制体单元体积 仅Couet是面积，而它的边界是长度。设 ，利用 定理，,wesnAyAxGaus可将方程（E.2）改写成如下有限体积离散格式：图 E.2 方腔流动计算网格、控制体单元和节点示意图 图 E.3 计算采用的交错网格示意图-E.3-ewnse nsuAvAASxyxxyy （E.3）对上式中 采用一阶向前差分近似，则有： ewnsxy、 、 、,PWEPewNSnsxxyy（E.4）同时记： ,eewwnnssFuAuAvv（E.5）,ewPEPWnsNSDxxAyy（E.6）则式（E.2）写成： ewnseEPwPWNP
5、SFFDxy（E.7）式中 都是主控制单元内节点上的已知量，如PEWSewns、 、 、 、 、 、 、 、果利用差分计算得到主控制单元边界上的流通量 ，就可ewnsFF、 、 、以求出节点上的未知量 。 、PEWNS、 、 、 、为了便于讨论，现对一维对流扩散方程的一阶迎风型离散格式进行分析。首先讨论无源项一维对流扩散方程的一阶迎风型离散格式。当流动为正向时，即 ，主控制单元界面取值：0,0,ewewuFWP-E.4-（E.8）则方程（E.7）离散为：wewPwWeEDFDF（E.9）当流动为负向时，即 ，控制体单元界面取值为：0,0,ewewu,wPeE（E.10）则方程（E.7）离散为：
6、weewPwWeEDFDF（E.11）将两种流动方向离散格式（E.9）和（E.11）合并后，得到统一的一维对流扩散方程一阶迎风型离散格式表达式：PWEaa（E.12）式中 ,0,WwPEewEeaDmxF（E.12a）同理，可以得到带有源项的二维对流扩散方程的一阶迎风型离散格式：PWESNaaaxy（E.13）式中 ,00,WwEeeSssNnnPESewsaDmxFaDmxF（E.13a）-E.5-源项 为： SupStx（E.14）若把 表示为 时刻的动量， 表示 时刻的动量，则可以得到源项()nunt 1()nu1nt离散格式为： S1nnPPewVSdxypyt（E.15）最后，得到带
7、有源项的二维对流扩散方程有限体积算法一阶迎风型显式离散格式：（E.161nnnnn nPPPWESN ewuauauaxypyt）式中系数 为一阶迎风格式中各对应系数。k5. 计算结果分析利用上述有限体积算法一阶迎风型离散格式和相应的初始条件和边界条件，采用 MAC 算法中压力耦合方程，求解二维不可压缩黏性流体 流动问题。Couet图 E.4 给出了二维不可压缩黏性流体 流动沿 y 方向的速度分布，并和精Couet确解进行了比较，十分吻合。图 E.5 是二维不可压缩黏性流体 流动水平tx 方向的速度云图。图 E.6 是二维不可压缩黏性流体 流动速度矢量分布ouet图。-E.6-E-2 二维不可
8、压缩黏性流体 流动问题的数值计算源程序Couet图 E.4 二维不可压缩黏性流体 流动沿 y 方向的速度分Couet布图 E.5 二维不可压缩黏性流体 流动水平 x 方向速度云图Couet图 E.6 二维不可压缩黏性流体 流动速度矢量分布图Couet-E.7-1. 语言源程序C/ fvm_upwind_MAC_couette.cpp/*-以一阶迎风型格式和 压力迭代求解二维 流动问题（C 语言版本）Chorinouet-*/#include #include #include #define MX 100 /最大网格数#define MY 20 /最大网格数#define Lambda 0.0
9、02 /Chorin 压力迭代常数，直接影响收敛性#define Re 100.00 /雷诺数#define dt 0.0005 /时间步长double uMX+1MY+2,vMX+2MY+1,pMX+2MY+2,unMX+1MY+2,vnMX+2MY+1;/全局变量double max(double a,double b) /定义 max 函数double c;if(ab)c=b;elsec=a;return c;/*-初始化输入: 无；输出：u、v、p、dx、dy ，初始速度、压强和空间步长。-*/void init(double uMX+1MY+2,double vMX+2MY+1,do
10、uble pMX+2MY+2,doubledx=1.0/MX;dy=0.2/MY;for(i=0;i=MX;i+)for(j=0;j=MY+1;j+)-E.8-uij=0.0;if(j=MY+1) uij=2.0;for(i=0;i=MX+1;i+)for(j=0;j=MY;j+)vij=0.0;for(i=0;i=MX+1;i+)for(j=0;j=MY+1;j+) pij=1.0;/*-迭代更新压强Chorin输入：u、v、p、dx、dy ，n 时刻速度、压强、空间步长；输出：p ，n+1 时刻压强。-*/void getp(double uMX+1MY+2,double vMX+2MY+
11、1,double pMX+2MY+2,double dx,double dy)double dMX+2MY+2;int i,j;for(i=1;i=MX;i+)for(j=1;j=MY;j+)dij=(uij-ui-1j)/dx + (vij-vij-1)/dy;for(i=1;i=MX;i+)for(j=1;j=MY;j+)pij=pij-Lambda*dij;for(i=1;i=MX;i+)pi0=pi1;piMY+1=piMY;for(j=1;j=MY;j+)p0j=p1j;pMX+1j=pMXj;/*-一阶迎风格式输入：u、v、p、dx、dy ，n 时刻速度、n+1 时刻压强,空间步长
12、；-E.9-输出：un、vn ，n+1 时刻速度。-*/void upwind(double uMX+1MY+2,double vMX+2MY+1,double pMX+2MY+2,double unMX+1MY+2,double vnMX+2MY+1,double dx,double dy)int i,j;double aw,ae,as,an,df,ap,miu;miu=1.0/Re; for(i=1;i=MX-1;i+)for(j=1;j=MY;j+)aw=miu+max(0.5*(ui-1j+uij)*dy,0.0);ae=miu+max(0.0,-0.5*(uij+ui+1j)*dy)
13、;as=miu+max(0.5*(vij-1+vi+1j-1)*dx,0.0);an=miu+max(0.0,-0.5*(vij+vi+1j)*dx);df=0.5*(ui+1j-ui-1j)*dy+0.5*(vij+vi+1j-vij-1-vi+1j-1)*dx;ap=aw+ae+as+an+df;unij=uij+dt/dx/dy*(-ap*uij+aw*ui-1j+ae*ui+1j+as*uij-1+an*uij+1)-dt*(pi+1j-pij)/dx;/边界条件for(i=1;i=MX-1;i+)uni0=-uni1;uniMY+1=2.0-uniMY;for(j=0;j=MY+1;
14、j+)un0j=un1j;unMXj=unMX-1j;for(i=1;i=MX;i+)for(j=1;j=MY-1;j+)aw=miu+max(0.5*(ui-1j+ui-1j+1)*dy,0.0);ae=miu+max(0.0,-0.5*(uij+uij+1)*dy);as=miu+max(0.5*(vij-1+vij)*dx,0.0);an=miu+max(0.0,-0.5*(vij+vij+1)*dx);df=0.5*(uij+uij+1-ui-1j-ui-1j+1)*dy+0.5*(vij+1-vij-1)*dx;ap=aw+ae+as+an+df;-E.10-vnij=vij+dt/
15、dx/dy*(-ap*vij+aw*vi-1j+ae*vi+1j+as*vij-1+an*vij+1)-dt*(pij+1-pij)/dy;/边界条件for(i=1;i=MX;i+)vni0=0.0;vniMY=0.0;for(j=0;j=MY;j+)vn0j=-vn1j;vnMX+1j=-vnMXj;/*-格式输出，采用 数据格式画图Tecplot输入：u、v、p、dx、dy ，最后结果；输出：无。-*/void output(double uMX+1MY+2,double vMX+2MY+1,double pMX+2MY+2,double dx,double dy)double uo,vo
16、,po;int i,j;FILE *fp;fp=fopen(“Result.plt“,“w“);fprintf(fp,“variables=x,y,u,v,pnzone i=%d,j=%d,f=pointrn“,MX+1,MY+1);for(j=0;j=MY;j+)for(i=0;i=MX;i+)uo=0.5*(uij+uij+1);vo=0.5*(vij+vi+1j);po=0.25*(pij+pi+1j+pij+1+pi+1j+1);fprintf(fp,“%20.10e %20.10e %20.10e %20.10e %20.10ern“,i*dx,j*dy,uo,vo,po);fclose(fp);
