2022年附录E有限体积法求解库塔流 .pdf

《2022年附录E有限体积法求解库塔流 .pdf》由会员分享，可在线阅读，更多相关《2022年附录E有限体积法求解库塔流 .pdf（15页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、-E.1- 附录 E 二维不可压缩黏性流体Couette流动问题的有限体积算法与计算程序二维 Couette流动是一个不可压缩黏性流动的典型流动，并有解析解，可以用来检验数值算法计算精度和可靠性。对它采用有限体积算法一阶迎风型离散格式进行数值求解。同时，为了初学者入门和练习方便, 这里给出了用 C 语言和Fortran77 语言编写的计算二维不可压缩黏性流体Couette流动问题计算程序，供大家学习参考。E-1 利用有限体积算法一阶迎风型离散格式求解二维不可压缩黏性流体Couette流动问题1. 二维不可压缩黏性流体Couette流动问题的提法二维不可压缩黏性流体Couette流动：有两个无限

2、长平板，间距为L ，两板之间充满密度为 1、静止的不可压缩黏性流体。无限长平板组成的通道两端压力相等，下板固定不动，上板以量纲为一的速度1 自左向右平移运动（图E.1） 。2. 基本方程组、初始条件和边界条件设流体是黏性流体。二维Couette流动问题在数学上可以由二维不可压缩黏性流动 N-S方程组来表示，把它写成通用变量的微分方程组形式，有：uvSxyxxyy（E.1）其中u为变量在水平x方向的流速，v为在垂直 y 方向的流速，为黏度，S为源项。源项中不仅包含压力梯度项，也包含时间导数项。初始条件：上板以量纲为一的速度1 沿着上壁面方向自左向右运动。图 E.1 二维不可压缩黏性流体Couet

3、te流动问题示意图精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 15 页-E.2- 边界条件：流动速度uv、在上下壁面可采用无滑移边界条件，在左右两端采用自由输出边界条件；压强p 采用自由输出边界条件。3. 计算网格划分和控制体单元与节点定义采用交错网格。图E.2 和图 E.3 是计算网格、控制体单元和节点示意图。节点 P 所在主控制体单元如图E.2 中有阴影部分所示。 在x方向与节点 P 相邻的节点为 W 和 E ，在 y 方向与节点 P 相邻的节点为 S和 N ，主控制单元界面分别为 w s e n、 。压力 p 和速度 u v、

4、 分别在三套不同网格中，如图E.3 中有阴影部分所示。4. 有限体积算法离散格式对方程（ E.1）在图 E.2 所示节点 P主控制单元内积分，有：VVVVVudVvdVdVdVS dVxyxxyy（E.2）由于不可压缩黏性流体Couette流动是二维问题， 因此，控制体单元体积V 仅是面积，而它的边界是长度。设,wesnAAy AAx，利用 Gauss定理，可将方程（ E.2）改写成如下有限体积离散格式：图 E.2 方腔流动计算网格、控制体单元和节点示意图图 E.3 计算采用的交错网格示意图精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共

5、 15 页-E.3- ewnsewnsu Au Av Av AAAAASx yxxyy（E.3）对上式中ewnsxxyy、采用一阶向前差分近似，则有：,PWEPewNPPSnsxxxxyyyy（E.4）同时记：,eewwewnnssnsFuAFuAFvAFvA（E.5）,eewwewPEPWnnssnsPNPSAADDxxAADDyy（E.6）则式（ E.2）写成：eewwnnsseEPwPWnNPsPSFFFFDDDDSx y（E.7）式中PEWNSewnsDDDD、都是主控制单元内节点上的已知量，如果利用差分计算得到主控制单元边界上的流通量eewwnnssFFFF、，就可以求出节点上的未知

6、量PEWNS、。 、为了便于讨论， 现对一维对流扩散方程的一阶迎风型离散格式进行分析。首先讨论无源项一维对流扩散方程的一阶迎风型离散格式。当流动为正向时，即0,0,0,0ewewuuFF，主控制单元界面取值：,wWeP（E.8）则方程（ E.7）离散为：wweewPwwWeEDFDFFDFD（E.9）当流动为负向时，即0,0,0,0ewewuuFF，控制体单元界面取值为：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 15 页-E.4- ,wPeE（E.10）则方程（ E.7）离散为：weeewPwwWeEDDFFFDFD（E.11）将

7、两种流动方向离散格式（E.9）和（ E.11）合并后，得到统一的一维对流扩散方程一阶迎风型离散格式表达式：PPWWEEaaa（E.12）式中,00,WwwPWEewEeeaDmax FaaaFFaDmaxF（E.12a）同理，可以得到带有源项的二维对流扩散方程的一阶迎风型离散格式：PPWWEESSNNaaaaaSx y（E.13）式中,0 ,0,0 ,0,WwwEeeSssNnnPWESNewnsaDmax FaDmaxFaDmax FaDmaxFaaaaaFFFF（E.13a）源项 S 为：upStx（E.14）若把()nu表示为nt时刻的动量，1()nu表示1nt时刻的动量，则可以得到源项

8、 S离散格式为：1nnPPewVuuS dVx yppyt（E.15）最后，得到带有源项的二维对流扩散方程有限体积算法一阶迎风型显式离散格式：1nnnnnnnnnPPPPWWEESSNNewuua ua ua ua ua ux yppyt（E.16）式中系数ka为一阶迎风格式中各对应系数。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 15 页-E.5- 5. 计算结果分析利用上述有限体积算法一阶迎风型离散格式和相应的初始条件和边界条件，采用 MAC算法中压力耦合方程，求解二维不可压缩黏性流体Couette流动问题。图 E.4 给出了二

9、维不可压缩黏性流体Couette流动沿 y 方向的速度分布，并和精确解进行了比较， 十分吻合。图 E.5 是二维不可压缩黏性流体Couette流动水平 x方向的速度云图。 图 E.6 是二维不可压缩黏性流体Couette流动速度矢量分布图。图 E.4 二维不可压缩黏性流体Couette流动沿 y 方向的速度分布图 E.5 二维不可压缩黏性流体Couette流动水平x 方向速度云图精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 15 页-E.6- E-2 二维不可压缩黏性流体Couette流动问题的数值计算源程序1. C 语言源程序/ f

10、vm_upwind_MAC_couette.cpp /*- 以一阶迎风型格式和 Chorin压力迭代求解二维 Couette流动问题（ C语言版本）-*/ #include #include #include #define MX 100 /最大网格数#define MY 20 /最大网格数#define Lambda 0.002 /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; /全局变量d

11、ouble max(double a,double b) / 定义 max 函数 double c; if(ab) c=b; else c=a; return c; /*- 初始化图 E.6 二维不可压缩黏性流体Couette流动速度矢量分布图精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 15 页-E.7- 输入 : 无；输出： u、 v、p、dx、dy ，初始速度、压强和空间步长。-*/ void init(double uMX+1MY+2,double vMX+2MY+1,double pMX+2MY+2,double& dx,d

12、ouble& dy) int i,j; dx=1.0/MX; dy=0.2/MY; for(i=0;i=MX;i+) for(j=0;j=MY+1;j+) 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

13、 vMX+2MY+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; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 15 页-E.8- piMY+1=piM

14、Y; for(j=1;j=MY;j+) p0j=p1j; pMX+1j=pMXj; /*- 一阶迎风格式输入： u、v、p、dx、dy ，n 时刻速度、 n+1 时刻压强 ,空间步长；输出： 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-

15、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); 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*ui

16、j-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;j+) un0j=un1j; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 15 页-E.9- 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+ui

17、j+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; vnij=vij+dt/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;

18、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,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); 精

19、选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 15 页-E.10- for(j=0;j=MY;j+) for(i=0;i1e-6 & step1e6) /err，收敛限制； step，迭代次数限制 err=0.0; getp(u,v,p,dx,dy); /更新压强upwind(u,v,p,un,vn,dx,dy); /更新速度for(i=0;i=MX;i+) for(j=0;jerr) err=value; uij=unij; for(i=0;i=MX+1;i+) for(j=0;jerr) err=value; vij=vnij; 精

20、选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 15 页-E.11- printf(err=%len,err); output(u,v,p,dx,dy); /输出结果- - 2.Fortran77 语言源程序program fvm_upwind_MAC_couette !- 以一阶迎风型离散格式和Chorin压力迭代求解二维 Couette流动问题（Fortran77 语言版本）!- parameter(mx=101,my=21) implicit double precision(a-h,o-z) dimension u(mx,my+

21、1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1),vn(mx+1,my) double precision lambda integer step re=100. !雷诺数dt=0.0005 !时间步长lambda=0.002 !Chorin 压力迭代常数，直接影响收敛性dx=1.0d0/(mx-1) dy=0.2d0/(my-1) step=0 err=100. call initial(u,v,p) !初始化do while(err1e-6.and.steperr) err=temp u(i,j)=un(i,j) enddo enddo do i=1,mx+1

22、do j=1,my temp=abs(v(i,j)-vn(i,j)/dt if(temperr) err=temp v(i,j)=vn(i,j) enddo enddo write(*,*) err=,err enddo call output(u,v,p,mx,my,dx,dy) !输出End !- !初始化!输入：无；!输出： u、v、p，初始速度、压强。!- subroutine initial(u,v,p) parameter(mx=101,my=21) double precision u(mx,my+1),v(mx+1,my),p(mx+1,my+1) do i=1,mx+1 do

23、 j=1,my+1 p(i,j)=1.0 enddo enddo do i=1,mx do j=1,my+1 u(i,j)=0.0 if(j=my+1) u(i,j)=2.0 enddo enddo do i=1,mx+1 do j=1,my v(i,j)=0.0 enddo enddo endsubroutine !- ! Chorin压力迭代!输入： u、v、p、mx、my、 dx、dy、lambda ，n 时刻速、压强、其他各量；精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 15 页-E.13- !输出： p ，n+1 时刻压

24、强。!- subroutine calp(u,v,p,mx,my,dx,dy,lambda) implicit double precision(a-h,o-z) dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),d(mx+1,my+1) double precision lambda do i=2,mx do j=2,my d(i,j)=(u(i,j)-u(i-1,j)/dx+(v(i,j)-v(i,j-1)/dy enddo enddo do i=2,mx do j=2,my p(i,j)=p(i,j)-lambda*d(i,j) enddo end

25、do do i=2,mx p(i,1)=p(i,2) p(i,my+1)=p(i,my) enddo do j=1,my+1 p(1,j)=p(2,j) p(mx+1,j)=p(mx,j) enddo endsubroutine !- !一阶迎风型离散格式!输入： u、v、pn、mx、my、dx、dy、dt、re ，n 时刻速度、 n+1 时刻压强、其他各量；!输出： un,vn，n+1 时刻速度。!- subroutine upwind(u,v,pn,un,vn,mx,my,dx,dy,dt,re) implicit double precision(a-h,o-z) dimension u

26、(mx,my+1),v(mx+1,my),pn(mx+1,my+1),un(mx,my+1),vn(mx+1,my) double precision miu miu=1.0/re do i=2,mx-1 do j=2,my aw=miu+max(0.5*(u(i-1,j)+u(i,j)*dy,0.0) ae=miu+max(0.0,-0.5*(u(i,j)+u(i+1,j)*dy) as=miu+max(0.5*(v(i,j-1)+v(i+1,j-1)*dx,0.0) an=miu+max(0.0,-0.5*(v(i,j)+v(i+1,j)*dx) df=0.5*(u(i+1,j)-u(i-

27、1,j)*dy+0.5*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1)*dx ap=aw+ae+as+an+df 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 15 页-E.14- un(i,j)=u(i,j)+dt/dx/dy*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j)+as*u(i,j-1)+an*u(i,j+1) & -dt*(pn(i+1,j)-pn(i,j)/dx enddo enddo !- do i=2,mx-1 un(i,1)=-un(i,2) un(i,my+1

28、)=2.0-un(i,my) enddo do j=1,my+1 un(1,j)=un(2,j) un(mx,j)=un(mx-1,j) enddo !u 边界条件!- !- do i=2,mx do j=2,my-1 aw=miu+max(0.5*(u(i-1,j)+u(i-1,j+1)*dy,0.0) ae=miu+max(0.0,-0.5*(u(i,j)+u(i,j+1)*dy) as=miu+max(0.5*(v(i,j-1)+v(i,j)*dx,0.0) an=miu+max(0.0,-0.5*(v(i,j)+v(i,j+1)*dx) df=0.5*(u(i,j)+u(i,j+1)-

29、u(i-1,j)-u(i-1,j+1)*dy+0.5*(v(i,j+1)-v(i,j-1)*dx ap=aw+ae+as+an+df vn(i,j)=v(i,j)+dt/dx/dy*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j)+as*v(i,j-1)+an*v(i,j+1) & -dt*(pn(i,j+1)-pn(i,j)/dy enddo enddo !- do i=2,mx vn(i,1)=0.0 vn(i,my)=0.0 enddo do j=1,my vn(1,j)=-vn(2,j) vn(mx+1,j)=-vn(mx,j) enddo !v 边界条件!- E

30、ndsubroutine !- !格式输出，采用Tecplot 数据格式画图!输入： u、v、p、mx、my、 dx、dy，最后结果；精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 15 页-E.15- !输出：无。!- subroutine output(u,v,p,mx,my,dx,dy) implicit double precision(a-h,o-z) dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1) dimension uc(mx,my),vc(mx,my),pc(mx,my),x(

31、mx),y(my) do i=1,mx x(i)=(i-1)*dx enddo do j=1,my y(j)=(j-1)*dy enddo do i=1,mx do j=1,my uc(i,j)=0.5*(u(i,j)+u(i,j+1) vc(i,j)=0.5*(v(i,j)+v(i+1,j) pc(i,j)=0.25*(p(i,j)+p(i+1,j)+p(i,j+1)+p(i+1,j+1) enddo enddo open(1,file=result.plt) write(1,*) TITLE = 2D Results write(1,*) V ARIABLES = X, Y, U, V, P write(1,*) ZONE I=,mx, J=,my, F=POINT write(1,10) (x(i),y(j),uc(i,j),vc(i,j),pc(i,j),i=1,mx),j=1,my) close(1) 10 format(5(e15.8,5x) endsubroutine - - - 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 15 页