利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言完整篇doc.docx
- 文档编号:4183576
- 上传时间:2022-11-28
- 格式:DOCX
- 页数:8
- 大小:16.07KB
利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言完整篇doc.docx
《利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言完整篇doc.docx》由会员分享,可在线阅读,更多相关《利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言完整篇doc.docx(8页珍藏版)》请在冰豆网上搜索。
利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言完整篇doc
利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言1
Forpersonaluseonlyinstudyandresearch;notforcommercialuseFortran77语言源程序!
————————————————————————————————————
!
利用有限体积算法三阶迎风型QUICK离散格式和
!
人工压缩算法求解方腔流动问题(Fortran77语言版本)!
————————————————————————————————————programQUICK_cavity
parameter(mx=101,my=101)
implicitdoubleprecision(a-h,o-z)
dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)
dimensionun(mx,my+1),vn(mx+1,my),pn(mx+1,my+1)
common/ini/u,v,p
c2=2.25
re=100.0
dt=0.0005
dx=1.0/float(mx-1)
dy=1.0/float(my-1)
!
----------------------------------------------------------------------------------------
!
u、v、p为t时刻值,un、vn、pn为t+1时刻值,
!
mx、my为最大网格数,c2为虚拟压缩系数的平方,re为雷诺数。
!
-----------------------------------------------------------------------------------------
num=0
err=100.00
!
nun,计数器;err,判断人工压缩法求解收敛的标准callinitial
!
调入初始条件,以下为人工压缩算法求解
dowhile(err.gt.1e-4.and.num.lt.1e6)
err=0.0
callquick(u,v,p,un,vn,mx,my,dx,dy,dt,re)
!
QUICK离散格式求解动量方程,得到un、vn
callcalp(p,un,vn,pn,mx,my,dx,dy,dt,c2)
!
求压强pn
callcheck(u,v,p,un,vn,pn,mx,my,dx,dy,dt,c2,err)!
校验流场信息,判断是否收敛,同时更新u、v、pwrite(*,*)'error=',err
num=num+1
write(*,*)num
!
屏幕跟踪输出
enddo
calloutput(u,v,p,mx,my,dx,dy)
!
输出结果文件
End
!
!
subroutineinitial
!
初始化流场
parameter(mx=101,my=101)
doubleprecisionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)common/ini/u,v,p
doi=1,mx+1
doj=1,my+1
p(i,j)=1.0
enddo
enddo
doi=1,mx
doj=1,my+1
u(i,j)=0.0
if(j.eq.my+1)u(i,j)=4.0/3.0
if(j.eq.my)u(i,j)=2.0/3.0
enddo
enddo
doi=1,mx+1
doj=1,my
v(i,j)=0.0
enddo
enddo
endsubroutine
!
!
subroutinequick(u,v,p,un,vn,mx,my,dx,dy,dt,re)
!
以QUICK格式离散动量方程
implicitdoubleprecision(a-h,o-z)
dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1),vn(mx+1,my)
doubleprecisionmiu
miu=1.0/re
!
以下求解x方向速度un----------------------------------------------------------------------------doi=3,mx-2
doj=3,my-1
fw=0.5*(u(i-1,j)+u(i,j))*dy
fe=0.5*(u(i,j)+u(i+1,j))*dy
fs=0.5*(v(i,j-1)+v(i+1,j-1))*dx
fn=0.5*(v(i,j)+v(i+1,j))*dx
df=fe-fw+fn-fs
aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fw
aww=-0.125*alfa(fw)*fw
ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fw
aee=0.125*(1.0-alfa(fe))*fe
as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fs
ass=-0.125*alfa(fs)*fs
an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fs
ann=0.125*(1.0-alfa(fn))*fn
ap=aw+ae+as+an+aww+aee+ass+ann+df
!
aw、ae、as、an...均为有限体积算法中各项系数,详见前文三阶迎风型QUICK离散格式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)+aww*u(i-2,j)+aee*u(i+2,j)&
+ass*u(i,j-2)+ann*u(i,j+2))-dt*(p(i+1,j)-p(i,j))/dx
enddo
enddo
!
-------------------------------------------------------------------------------------------
j=2
doi=3,mx-2
callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)
enddo
j=my
doi=3,mx-2
callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)
enddo
i=2
doj=2,my
利用有限体积算法三阶迎风型QUICK离散格式和人工压缩算法求解方腔流动问题Fortran语言1第2页
callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)
enddo
i=mx-1
doj=2,my
callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)
enddo
!
内层边界由一阶迎风型离散格式得到----------------------------------------------------!
-------------------------------------------------------------------------------------------doi=2,mx-1
un(i,1)=-un(i,2)
un(i,my+1)=2.0-un(i,my)
enddo
doj=1,my+1
un(1,j)=0.0
un(mx,j)=0.0
enddo
!
外层边界条件按物理边界条件给出
!
-----------------------------------------------------------------------------
!
以下求解y方向速度vn-----------------------------------------------
doi=3,mx-1
doj=3,my-2
fw=0.5*(u(i-1,j)+u(i-1,j+1))*dy
fe=0.5*(u(i,j)+u(i,j+1))*dy
fs=0.5*(v(i,j-1)+v(i,j))*dx
fn=0.5*(v(i,j)+v(i,j+1))*dx
df=fe-fw+fn-fs
aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fw
aww=-0.125*alfa(fw)*fw
ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fw
aee=0.125*(1.0-alfa(fe))*fe
as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fs
ass=-0.125*alfa(fs)*fs
an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fs
ann=0.125*(1.0-alfa(fn))*fn
ap=aw+ae+as+an+aww+aee+ass+ann+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)+aww*v(i-2,j)+aee*v(i+2,j)&
+ass*v(i,j-2)+ann*v(i,j+2))-dt*(p(i,j+1)-p(i,j))/dy
enddo
enddo
!
-----------------------------------------------------------------------------
j=2
doi=3,mx-1
callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)
enddo
j=my-1
doi=3,mx-1
callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddo
i=2
doj=2,my-1
callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddo
i=mx
doj=2,my-1
callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddo
!
----------------------------------------------------------------------------doi=2,mx
vn(i,1)=0.0
vn(i,my)=0.0
enddo
doj=1,my
vn(1,j)=-vn(2,j)
vn(mx+1,j)=-vn(mx,j)
enddo
!
----------------------------------------------------------------------------Endsubroutine
!
!
functionalfa(x)
!
函数,正1负0
doubleprecisionalfa,x
if(x.gt.0.d0)then
alfa=1.0
else
alfa=0.0
endif
end
!
!
subroutineupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)
!
以一阶迎风型离散格式得到内层边界un的值
implicitdoubleprecision(a-h,o-z)
dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1)
doubleprecisionmiu
miu=1.0/re
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-1,j))*dy+0.5*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1))*dxap=aw+ae+as+an+df
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*(p(i+1,j)-p(i,j))/dxEndsubroutine
!
!
subroutineupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)
!
以一阶迎风型离散格式得到内层边界vn值
implicitdoubleprecision(a-h,o-z)
dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),vn(mx+1,my)
doubleprecisionmiu
miu=1.0/re
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)-u(i-1,j)-u(i-1,j+1))*dy+0.5*(v(i,j+1)-v(i,j-1))*dxap=aw+ae+as+an+df
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 利用 有限 体积 算法 迎风 QUICK 离散 格式 人工 压缩 求解 流动 问题 Fortran 语言 完整 doc
![提示](https://static.bdocx.com/images/bang_tan.gif)
链接地址:https://www.bdocx.com/doc/4183576.html