附录D 二维泊萧叶黏性流动 (D).doc

附录D 二维不可压缩黏性流体流动问题的数值解与计算程序二维流动问题是有解析解的二维不可压缩黏性流动。对它采用算法和压力迭代解法进行数值求解。同时，为了初学者入门和练习方便,这里给出了由语言和语言编写的计算二维不可压缩黏性流体流动问题的计算程序，供大家学习参考。D-1 算法和压力迭代解法求解二维不可压缩黏性流动问题1. 二维不可压缩黏性流体流动问题设在水平方向上有两块无限长固定不动的平行平板，它们的间距为，两平板间充满不可压缩黏性流体。平板间两个横截面和上压力分别为和，当和不同时，平板间不可压缩黏性流体就会产生流动,并在平板间形成一个速度分布剖面，这就是著名的二维不可压缩黏性流体流动问题，如图D.1和图D.2所示。假定忽略质量力，且认为流动是定常层流。在平板间的横截面上黏性流体的速度分布精确解为： （D.1）图D.2二维流动横截面速度分布示意图图D.1二维流动示意图2. 基本方程组、初始条件和边界条件二维流动问题在数学上可以用二维不可压缩黏性流动方程组来描述：连续方程： （D.2）动量方程： （D.3）其中和分别为和方向的速度分量，是压力，是雷诺数。初始条件：。边界条件：左边界（）和右边界（）均为出口边界，速度应满足输出边界条件：；左边界上压力为：，右边界上压力为：；上边界（）和下边界（）均为刚性壁面，应满足无滑移边界条件：。3. 算法和压力迭代解法差分格式网格采用交错网格。计算采用算法。节点和变量设置如图D.3所示。图D.3 算法交错网格的节点坐标位置关系示意图表D.1在算法交错网格中变量下标范围设置方向方向在算法中由于采用交错网格，速度和压力分别设置在三套不同的网格节点上。压力值设置在主控体积单元的中心；在方向，速度设置在主控体积单元东、西两个侧面的辅助体积单元上，而在方向，速度设置在主控体积单元南、北两个侧面的辅助体积单元上。为了便于处理和的边界条件，在边界四周加设一层虚拟网格。三个变量的下标范围如表D.1所示：对方程（D.3）中对流项和扩散项采用中心差分格式。动量方程中的差分格式为： （D.4）其中 （D.5）动量方程中的差分格式为： （D.6）其中 （D.7） 压力差分格式为： （D.8）其中 （D.9）当在所有网格点中满足下面条件时，则可得到定常解： （D.10）式中是一小量，一般取。压力和速度之间迭代过程采用压力迭代解法。具体迭代过程如下：首先由式（D.1）差分格式计算出值当m = 1时，、用、值来代替，求出值，其中m为迭代次数；然后利用压力修正公式： （D.11）求出当m = 1时，用值来代替，计算出值，其中是松弛因子，必须满足如下稳定性条件： （D.12）为了减少迭代次数，通常可取最佳值；最后利用动量方程、差分格式求出和值，并利用散度差分格式，重新计算，检验散度是否满足，如果不满足，则需要重复计算，直到收敛。4. 计算结果分析计算分别采用标准的语言和语言编写程序。计算中网格点数为10121。取。在的情况下，通过计算可以求出速度分布为： （D.14）图D.4 采用语言程序计算得到的x = 0.5处速度剖面与精确解比较的结果图D.5 采用语言程序计算得到的x = 0.5处的速度剖面与精确解比较结果 图D.4给出了采用语言程序计算得到的处速度剖面上计算结果与精确解比较的结果；图D.5给出了采用语言程序计算得到处速度剖面上计算结果与精确解比较的结果。从图D.4和图D.5中可以看出，采用算法和压力迭代解法所得到的计算结果和精确解是十分吻合的。图D.6和图D.8分别给出了采用语言和语言编制的程序，计算得到的水平速度云图；图D.7和图D.9分别给出了采用语言和语言编制的程序，计算得到的水平速度矢量图。这些计算结果表明计算已经达到定常状态，计算结果真实地反映了二维不可压缩黏性流体流动。图D.6 采用语言程序计算得到的水平速度云图图D.7 采用语言程序计算得到的水平速度矢量图图D.8 采用语言程序计算得到的水平速度云图RAN程序计算的水平速度云图图D.9 采用语言程序计算得到的水平速度矢量图D-2 二维不可压缩黏性流体问题的数值计算源程序1. 语言源程序/ MAC-Chorin2D.cpp : 定义控制台应用程序的入口点。/*-*利用算法和压力迭代解法求解二维不可压缩黏性流动问题(语言版本) *-*/采用Visual Studio 2003写的标准C程序。如果在其他平台编译有问题，请利用Visual Studio 2003建立一个控制台工程进行编译。-*/#include stdafx.h#include stdio.h#include math.h#define Jx 101 /x方向网格点数#define Jy 21 /y方向网格点数#define RE 100 /Reynolds数#define DT 0.0002 /时间步长#define EPS 1e-6 /收敛限/全局变量double uJx+2Jy+2,vJx+2Jy+2,pJx+2Jy+2,tuJx+2Jy+2,tvJx+2Jy+2;/*-速度边界条件入口: 无；出口: u、v 为已经给定的边界，u、v的上下边界为无滑移固壁，u、v的左右边界为出口边界。-*/void bounduv(double uJx+2Jy+2,double vJx+2Jy+2) int i,j; for(i=0;i=Jx+1;i+) ui1=0; uiJy=0; vi0=-vi1; viJy=-viJy-1; for(j=1;j=Jy;j+) u0j=u1j; uJxj=uJx-1j; v1j=v1j; vJxj=vJx-1j; /*-压力边界条件入口: 无；出口: p为已经给定的边界，只需要给定左右压力，左侧1.0,右侧0.0。-*/ void boundp(double pJx+2Jy+2) int j; for(j=0;j=Jy+1;j+) p1j=1.0; pJxj=0.0; /*-初始化入口: 无；出口: u、v、p，已经给定的初始值。-*/void init(double uJx+2Jy+2,double vJx+2Jy+2,double pJx+2Jy+2,double dx,double dy,double Re) int i,j; for(i=0;i=Jx+1;i+) for(j=0;j=Jy+1;j+) uij=0; vij=0; pij=0; bounduv(u,v); boundp(p);/*-根据n时刻u、v、p采用中心差分格式计算n+1时刻u、v入口: u、v、p，n时刻u、v、p，tu、tv，临时变量，dx、dy，网格宽度，dt，时间步长，Re，雷诺数；出口: u、v，n+1时刻u、v。-*/void solveuv(double uJx+2Jy+2,double vJx+2Jy+2,double pJx+2Jy+2,double dx,double dy,double dt,double Re,double tuJx+2Jy+2, double tvJx+2Jy+2) double vav,uav,adv,vis,prs; int i,j; for(i=1;i=Jx-1;i+) for(j=2;j=Jy-1;j+) vav=(vij+vi+1j+vij-1+vi+1j-1)/4.0; adv=-uij*(ui+1j-ui-1j)/dx/2.0-vav*(uij+1-uij-1)/dy/2.0; vis=(ui+1j-2*uij+ui-1j)/dx/dx+(uij+1-2*uij+uij-1)/dy/dy)/Re; prs=-(pi+1j-pij)/dx; tuij=(adv+vis+prs)*dt+uij; for(i=2;i=Jx-1;i+) for(j=1;j=Jy-1;j+) uav=(ui-1j+1+uij+1+ui-1j+uij)/4.0; adv=-uav*(vi+1j-vi-1j)/dx/2.0-vij*(vij+1-vij-1)/dy/2.0; vis=(vi+1j-2*vij+vi-1j)/dx/dx+(vij+1-2*vij+vij-1)/dy/dy)/Re; prs=-(pij+1-pij)/dy; tvij=(adv+vis+prs)*dt+vij; for(i=1;i=Jx;i+) for(j=1;j=Jy;j+) uij=tuij; vij=tvij; bounduv(u,v);/*-根据n时刻u、v、p计算n+1时刻p入口: u、v、p，n时刻u、v、p，lambda，常数；出口: p，n+1时刻p。-*/void solvep(double pJx+2Jy+2,double uJx+2Jy+2,double vJx+2Jy+2,double dx,double dy,double lambda) int i,j; for(i=1;i=Jx;i+) for(j=1;j=Jy;j+) pij=pij-lambda*(uij-ui-1j)/dx+(vij-vij-1)/dy); boundp(p);/*-判断收敛入口: u、v，当前时刻u、v，dx、dy，网格宽度，eps，收敛限制；出口: dm， 当前网格点上速度散度最大值；返回: 1达到收敛，0尚未收敛。-*/int conv(double uJx+2Jy+2,double vJx+2Jy+2,double dx,double dy,double eps,double &dm) int i,j; double er; dm=0; for(i=1;i=Jx;i+) for(j=1;j=Jy;j+) er=fabs(uij-ui-1j)/dx+(vij-vij-1)/dy); if(dmer)dm=er; if(dmeps)return 1; else return 0;/*-输出全场计算结果和x=0.5处的速度u剖面入口: u、v、p，需要输出的u、v、p，dx、dy，网格宽度；出口: 无。-*/void output(char *sN,double uJx+2Jy+2,double vJx+2Jy+2,double pJx+2Jy+2,double dx,double dy) int i,j; FILE *fp; fp=fopen(sN,w); fprintf(fp,Title = computational Results n); fprintf(fp,VARIABLES = x, y, u,v,p n); fprintf(fp,ZONE T= Zone 1, I=%d, J=%d,F=POINTn,Jy,Jx); for(i=1;i=Jx;i+) for(j=1;j=Jy;j+) fprintf(fp,%16f%16f%20e%20e%20en,(i-1)*dx,(j-1)*dy,(ui-1j+uij)/2.0,(vij-1+vij)/2.0,pij); fclose(fp); fp=fopen(result.txt,w); i=(Jx+1)/2; for(j=1;j=Jy;j+) fprintf(fp,%16f%20e%20en,(j-1)*dy-0.1,uij,pij); fclose(fp);void main() int n,isconv; double dx,dy,lambda; double dm; dx=1.0/(Jx-1); /x方向网格间距 dy=0.2/(Jy-1); /y方向网格间距 init(u,v,p,dx,dy,RE); /*- lambda常数与时间步长、网格间距和Reynolds数都有关。 关系到收敛性，应仔细选择！ -*/ lambda=2e-3; n=0; do solvep(p,u,v,dx,dy,lambda); solveuv(u,v,p,dx,dy,DT,RE,tu,tv); isconv=conv(u,v,dx,dy,EPS,dm); n+; if(n%100=0) printf(%5d Step Max vel divergence=%16en,n,dm); while(!isconv); output(result.plt,u,v,p,dx,dy);-2. 语言源程序! MAC-Chorin2D.for /*-*利用算法和压力迭代解法求解二维不可压缩黏性流动问题(语言版本) *-*/ program MAC_Chorin implicit double precision (a-h,o-z) parameter (M1=101,M2=21) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension u(0:M1+1,0:M2+1),v(0:M1+1,0:M2+1),p(0:M1+1,0:M2+1) dimension tu(0:M1+1,0:M2+1),tv(0:M1+1,0:M2+1)!x方向网格点数 Jx=M1!y方向网格点数 Jy=M2!Reynolds数 Re=100!x方向网格间距 dx=1.0/(Jx-1)!y方向网格间距 dy=0.2/(Jy-1)!时间步长 dt=0.0002!收敛限 EPS=1d-6 call init(u,v,p)!lambda常数与时间步长、网格间距和Reynolds数都有关。关系到收敛性，应仔细选择！ dlambda=2d-3; n=0!计算压力1 call solvep(p,u,v,dlambda)!计算速度 call solveuv(u,v,p,tu,tv)!判断是否满足连续性条件 isconv=iconv(u,v,dm) n=n+1 if(mod(n,100).eq.0)write(*,*)n,Step Max Vel Divergence=,dm if(isconv.eq.0)goto 1 call output(u,v,p) end!-!速度边界条件!入口: 无；!出口: u、v，已经给定边界，u、v上下边界为无滑移固壁，u、v左右边界为!出口边界。!- subroutine bounduv(u,v) implicit double precision (a-h,o-z) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension u(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1) do 10 i=0,Jx+1 u(i,1)=0; u(i,Jy)=0; v(i,0)=-v(i,1); v(i,Jy)=-v(i,Jy-1);10 continue do 11 j=1,Jy u(0,j)=u(1,j); u(Jx,j)=u(Jx-1,j); v(1,j)=v(1,j); v(Jx,j)=v(Jx-1,j);11 continue end!-!压力边界条件!入口: 无；!出口: p，为已经给定边界，只需要给定左右压力，左侧1.0,右侧0.0。!- subroutine boundp(p) implicit double precision (a-h,o-z) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension p(0:Jx+1,0:Jy+1) do 20 j=0,Jy+1 p(1,j)=1.0; p(Jx,j)=0.0;20 continue end !-!初始化!入口: 无；!出口: u、v、p，已经给定的初始值。!- subroutine init(u,v,p) implicit double precision (a-h,o-z) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension u(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1) do 30 i=0,Jx+1 do 30 j=0,Jy+1 u(i,j)=0; v(i,j)=0; p(i,j)=0;30 continue call bounduv(u,v) call boundp(p); end!-!根据n时刻u、v、p.采用中心差分格式计算n+1时刻u、v!入口: u、v、p，n时刻u、v、p，tu、tv临时变量；!出口: u、v，n+1时刻u、v。!- subroutine solveuv(u,v,p,tu,tv) implicit double precision (a-h,o-z) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension u(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1) dimension tu(0:Jx+1,0:Jy+1),tv(0:Jx+1,0:Jy+1) do 40 i=1,Jx-1 do 40 j=2,Jy-1 vav=(v(i,j)+v(i+1,j)+v(i,j-1)+v(i+1,j-1)/4.0; adv=-u(i,j)*(u(i+1,j)-u(i-1,j)/dx/2.0 $ -vav*(u(i,j+1)-u(i,j-1)/dy/2.0; vis=(u(i+1,j)-2*u(i,j)+u(i-1,j)/dx/dx $ +(u(i,j+1)-2*u(i,j)+u(i,j-1)/dy/dy)/Re; prs=-(p(i+1,j)-p(i,j)/dx; tu(i,j)=(adv+vis+prs)*dt+u(i,j);40 continue do 41 i=2,Jx-1 do 41 j=1,Jy-1 uav=(u(i-1,j+1)+u(i,j+1)+u(i-1,j)+u(i,j)/4.0; adv=-uav*(v(i+1,j)-v(i-1,j)/dx/2.0 $ -v(i,j)*(v(i,j+1)-v(i,j-1)/dy/2.0; vis=(v(i+1,j)-2*v(i,j)+v(i-1,j)/dx/dx $ +(v(i,j+1)-2*v(i,j)+v(i,j-1)/dy/dy)/Re; prs=-(p(i,j+1)-p(i,j)/dy; tv(i,j)=(adv+vis+prs)*dt+v(i,j);41 continue do 42 i=1,Jx do 42 j=1,Jy u(i,j)=tu(i,j); v(i,j)=tv(i,j);42 continue call bounduv(u,v) end!-!根据n时刻u、v、p计算n+1时刻p!入口: u、v、p，n时刻u、v、p；dlambda，常数；!出口: p，n+1时刻p。!- subroutine solvep(p,u,v,dlambda) implicit double precision (a-h,o-z) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension u(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0:Jy+1) do 50 i=1,Jx do 50 j=1,Jy p(i,j)=p(i,j) $ -dlambda*(u(i,j)-u(i-1,j)/dx+(v(i,j)-v(i,j-1)/dy);50 continue call boundp(p) end!-!判断收敛!入口: u、v，当前时刻u、v、p；!出口: dm，当前网格点上速度散度最大值，!返回: 1达到收敛，0尚未收敛。!- integer function iconv(u,v,dm) implicit double precision (a-h,o-z) common /G_def/ Jx,Jy,dx,dy,dt,Re,EPS dimension u(0:Jx+1,0:Jy+1),v(0:Jx+1,0:Jy+1),p(0:Jx+1,0: