高等计算流体力学-08

• real :: basel,basev,baser,baset,basem ！无量刚基准， 长度-速度-密度-温度-马赫数 • Real:: ro(iq,jq),p(iq,jq),roe(iq,jq),tp(iq,jq), • rovx(iq,jq),rovy(iq,jq),vx(iq,jq),vy(iq,jq),ho(iq,jq), & • ！ ρ，p， ρE（单位体积总能），T（温度）， ρu， ρv，u， v，总焓 • &, hro(iq,jq),hre(iq,jq),hrx(iq,jq),hry(iq,jq)& • ！ 通量代数和（方程右端项）连续-能量-动量 • &,rom1(iq,jq),roem1(iq,jq),rovxm1(iq,jq),rovym1(iq,jq) ！ N时刻守恒变量
• real :: vib(iq,jq),vjb(iq,jq) ！控制体界面中点运动的法向速度（动网格 计算中有用） • real :: & • &uil (iq,jq),uir(iq,jq),& 控制体i-界面左右两侧 • &vil (iq,jq),vir(iq,jq),& 的速度、压力和密度 • &pil (iq,jq),pir(iq,jq),& • &ril (iq,jq),rir(iq,jq),&
integer ::nmax,n,irest, ioup,iinp, icgl, iimplicit, idgstart,nts1p, nunsloop,ninnerloop, ip,iprint,istart,nb ！ nmax: maximum advancing steps； ! irest: start from very beginning(0) or start from a previous computation(1); ! iprint: output every iprint steps ！ icgl=1: gloal time step; =0 local time step ！ip=1输出计算结果； n：current time-step（累积）； istart： current time-step（当前计算）
R φL φ φi 1/2, j ,k i 1/2, j , k i 1/2, j , k
MUSCL插值
计算左右状态 ˆ 针对每个控制体进行限制,限制后物理量计为 φ
ˆ iL1/2, j ,k φi , j ,k 0.5B(2φiL1/2, j ,k ,2φiR1/2, j ,k ) φ ˆ iR1/2, j ,k φi , j ,k 0.5B(2φiR1/2, j ,k ,2φiL1/2, j ,k ) φ
问题：双马赫反射
激波前
ro(i,j)=1.4 vx(i,j)=0. vy(i,j)=0. p(i,j)=1.0 激波后 ro(i,j)=8.0 vx(i,j)=7.1447 vy(i,j)=-4.125 p(i,j)=116.5
1.2
1
0.8
y
0.6
0.4
0.2
0.5
1
1.5
2
2.5
3
3.5
x
数值解
φ [u, v, w, p, T ]T φi 1/2, j ,k ri 3/2, j ,k ri 1/2, j ,k φi , j ,k ri 1/2, j ,k ri 1/2, j ,k φi 1, j ,k ri 1/2, j ,k ri 1/2, j ,k ri 3/2, j ,k ri 1/2, j ,k
im,jm
数据结构
ib=3,jb=3
i+1/2,j+1/2
i+1,j+1
i+1/2,j i-1/2,j+1-/2,j
i,j
i,j-1/2 i-1/2,j-1/2 i,j i,j
变量定义
• module main • parameter(iq=405,jq=105,ngmax=max(iq,jq)+5) !数组维数 • parameter(small=1.0d-20,small1=1.0d-10,alarge=1.0d20) ! 小量和大量 • parameter(one=1.0,zero=0.0) ! 1和0 • logical :: steady1 ！steady1: logical variable, ==.true.
!
!
虚拟网格设置
• bc_ghostcell_value
– 除固体壁面边界条件外,其 他类型边界通过线性外推 得到虚拟网格点上的值;
i=ib-1 ro(i,j)=max(2.*ro(i+1,j)-ro(i+2,j),1.0e-10) vx(i,j)=2.*vx(i+1,j)-vx(i+2,j) vy(i,j)=2.*vy(i+1,j)-vy(i+2,j) p(i,j)=max(2.*p(i+1,j)-p(i+2,j),1.0e-10) vmu(i,j)=2.*vmu(i+1,j)-vmu(i+2,j) vmul(i,j)=2.*vmul(i+1,j)-vmul(i+2,j) alagmx(i,j)=alagmx(i+1,j) j=jb xhalf=0.5*(x(i,j)+x(i+1,j)) if(xhalf>=1./6.) then !x<1/6,均匀来流条件，否则固壁条件 vx(i,j-1)=vx(i,j) vy(i,j-1)=-vy(i,j) p (i,j-1)=p (i,j) ro(i,j-1)=ro(i,j) vmu(i,j-1)=vmu(i,j) vmul(i,j-1)=vmul(i,j) alagmy(i,j-1)=alagmy(i,j)
muscl_interpolation
• • • • • dr=(uil(i+1,j)- vx(i,j))*2. dl=(vx(i,j) -uir(i,j))*2. uir(i,j) =vx(i,j)-0.5*alimiter(dl,dr) uil(i+1,j)=vx(i,j)+0.5*alimiter(dr,dl) • • • • • dr=(ujl(i,j+1)- vx(i,j))*2. dl=(vx(i,j) -ujr(i,j))*2. ujr(i,j) =vx(i,j)-0.5*alimiter(dl,dr) ujl(i,j+1)=vx(i,j)+0.5*alimiter(dr,dl)
– 固体壁面外虚拟网格上物 理量的设置要使经 interpolation_to_interface处理 后的边界上的值满足固壁 边界条件 (其实在无粘计 算中也可简单外推, 但如果 求解的是NS方程，粘性通 量计算在MUSCL插值之前 进行，要求壁面上的值满 足边界条件)
计算界面状态
• interpolation_to_interface • alsinv=1./(alagmx (i,j)+alagmx (i-1,j)) – 采用线性插值计算界面状 态, 考虑网格尺度影响。 • uil (i,j)=(alagmx (i,j)*vx – 左右状态同为 φi1/2, j,k (i-1,j)+alagmx (i-1,j)*vx (i,j)) *alsinv
• • • •
&ujl (iq,jq),ujr(iq,jq),& &vjl (iq,jq),vjr(iq,jq),& &pjl (iq,jq),pjr(iq,jq),& &rjl (iq,jq),rjr(iq,jq)
控制体j-界面左右两侧 的速度、压力和密度
• • • • • • • • • • •
1 x( y 2 2 2 ) y(2 x 2 2 ) B( x, y ) sign( x ) sign( y ) 2 2 x 2 xy 2 y 2 3 2
φiL1/2, j ,k φiL1/2, j,k φi, j, k φiR1/2, j ,k φi , j ,k φiR1/2, j,k
compute steady flow; ==.false. the flow is unsteady • integer :: ib,jb,im,jm,irsolver ！ib,jb,im,jm网格指标； irsolver， irsolver=0 Roe; =1 Lax -riemann solver
• real :: cfl,ttime,tstep,epsall • ! Cfl数，计算时间，时间步长，收敛标准 • real ::cp,ga,pi,fga,ga1,rcp,rcpcv,rfga,cv, prl,prt,prte • ！Cp，γ，π，ga1=ga-1，fga=(ga-1.0)/ga，rfga=1.0/fga， rcp=1.0/cp， cv=cp/ga， rcpcv=cp-cv，prl,prt,prte：层流、湍流、 有效普朗特数 • real :: vmul(iq,jq),vmu(iq,jq) ！层流及有效粘性系数 • real :: vxin,vyin,tpin,pin,pout！远场（进出口）参数 • real :: rms1,rms2,rms3,rms4 ！连续、动量、能量方程的最大残 差。 • integer :: imax1,jmax1, imax2,jmax2, imax3,jmax3, imax4,jmax4 ！ 连续、动量、能量方程的最大残差位置。 • end module main
数值解
程序结构
• 程序结构
核心算法
• • • • • • • • • • • subroutine solver(nrk) use main CALL CPU_TIME ( time1 ) ! set ghost-cell at boundary for various boundary conditions call bc_ghostcell_value ! interpolate the primitive variables to cell interfaces call interpolation_to_interface ! MUSCL interpolation at cell interface except boundaries call muscl_interpolation ! set muscl_interpolated value at the cell interfaces on the boundary call bc_muscl_interpolation • • • • • • • • • • • • • • • ! evaluate the inviscid fluxes if(irsolver==0) then call inviscid_fluxes_roe else call inviscid_fluxes_lax end if Runge-Kutta time stepping call Runge_Kutta(nrk) update variables call update_variables CALL CPU_TIME ( time2 ) write(*,*) 'time spent' write(*,*) time2-time1 return end
其中为沿的质点导数沿有riemann不变量其中为沿的质点导数均熵假定沿有riemann不变量二固壁边界法向量指向求解域内部有限差分具体实施过程设满足边界条件的边界点数值解为它和之间通过三条指向域外的特征线对应的相容关系相联系具体实施过程应用边界条件特征关系边界条件求解有限体积
第八讲 求解二维Euler方程的计算程序
• real :: x(iq,jq),y(iq,jq)&,xc(iq,jq),yc(iq,jq),step(iq,jq),vol(iq,jq) • ！控制体角点和中心点坐标，时间步长，控制体面 积 • real :: aix(iq,jq),aiy(iq,jq),ajx(iq,jq),ajy(iq,jq),& • ！控制体面积矢量 Snx，Sny • &aim(iq,jq),ajm(iq,jq) & • ! 控制体面积大小 • &,alagmx(iq,jq),alagmy(iq,jq) • 控制体长度和宽度
二十二、程序介绍
基本特征
• • • • • 结构网格二阶有限体积方法 重构:准一维重构+TVD-MUSCL插值 Riemann Solver: Lax, Roe 时间推进: 二阶TVD Runge-Kutta 边界处理:基于Roe格式的简化特征边界条件
定义数组维数时imax=im+s,jmax=jm+s, s为虚拟网格层数