logo资料库

流过平面的超声速流动数值模拟.docx

发布时间：2022-06-08 发布人：admin 分类：说明书 资料大小：0.02M 资料格式：docx 举报 版权申诉
第1页 / 共17页
第2页 / 共17页
第3页 / 共17页
第4页 / 共17页
第5页 / 共17页
第6页 / 共17页
第7页 / 共17页
第8页 / 共17页
资料共17页，剩余部分请下载后查看
    !平板上的有粘超声速流动：完整的 N-S 方程的 MacCormack 数值解法 !变量说明： !IMAX:网格的 X 方向个数 !JMAX:网格的 Y 方向个数 !MAXIT:循环推进限制的最大次数 !ITER：循环推进控制变量 !MACH_FAR:远场的马赫数 !V_FAR:远场的来流速度 !TEMP_FAR:远场的温度 !A_FAR:远场的声速 !P_FAR:远场的压力 !ROU_FAR:远场的密度 !MIU_FAR:远场的动力粘性系数 !K_FAR:远场的热传导系数 !E_FAR:远场的单位微团内能 !LHORI:平板长度 !TW:物面温度（假设物面为常壁面温度条件） !GAMMA:比热比 !PR:普朗特数 !MIU_SEA:标准海平面的动力粘性系数 !TEMP_SEA:标准海平面的温度 !R:气体常数 !CV：定容比热容 !CP:定压比热容 !REL:初始雷诺数 !DELTA:层流边界层厚度 !LVERT：计算域 Y 方向的高度 !DX,DY:X,Y 方向步长 !IE:标示调用 E/F 中粘性应力项及热传导项的时刻 !U、E、F：N-S 程原始的通量项的列向量形式 !U1:U1 项初始时刻值及 DT 时间后校正值(U/E/F 中的其他各分量同)，作为循环推进不断更 新 !U1_ES:U1 项 T+DT 时刻的预测值(U/E/F 中的其他各分量同) !VX,VY,ROU,P,E,TEMP,MIU,K:原始变量，作为循环不断推进更新变化 !ROU_PRE:用来存储前一时刻的密度值，与当前时刻的密度值比较后来作为判断是否达到 收敛的条件 !定义常量模块 MODULE GLOBAL IMPLICIT NONE REAL::TW0,GAMMA,PR,MIU_SEA,R,TEMP_SEA,CV,CP PARAMETER(TW0=1.0,GAMMA=1.4,PR=0.71,MIU_SEA=1.7894E-5,R=287.,TEMP_SEA=28 8.16,CV=717.5,CP=1004.5)
    END MODULE !主程序 PROGRAM MAIN USE GLOBAL IMPLICIT NONE INTEGER::IMAX,JMAX,MAXIT,ITER,I,J REAL::MACH_FAR,TEMP_FAR,A_FAR,P_FAR,ROU_FAR,V_FAR,LHORI,TW,REL,E_FAR, MIU_FAR,K_FAR,DELTA,LVERT,DX,DY,DT REAL,ALLOCATABLE::VX(:,:),VY(:,:),A(:,:),ROU(:,:),P(:,:),E(:,:),TEMP(:,:),MIU(:,:),K(:,:),U1 (:,:),U2(:,:),U3(:,:),U5(:,:),& E1(:,:),E2(:,:),E3(:,:),E5(:,:),F1(:,:),F2(:,:),F3(:,:),F5(:,:),ROU_PRE(:,:) REAL,EXTERNAL::DYNVIS,THERMC INTEGER,EXTERNAL::CONVER OPEN(6,FILE='OUTPUT.DAT') OPEN(5,FILE='INPUT.DAT') !WRITE(*,"('请输入划分网格 X 方向个数：')") READ(5,*) IMAX !WRITE(*,"('请输入划分网格 Y 方向个数：')") READ(5,*) JMAX WRITE(6,101) IMAX,JMAX 101 FORMAT(5X,'***** !WRITE(*,"('请输入循环推进最大次数限制：')") READ(5,*) MAXIT WRITE(6,102) MAXIT 102 FORMAT(5X,'***** 循环推进最大次数限制 *****'/,8x,'MAXIT=',I4/) !WRITE(*,"('请输入来流马赫数：')") READ(5,*) MACH_FAR !WRITE(*,"('请输入来流温度：')") READ(5,*) TEMP_FAR !WRITE(*,"('请输入来流（海平面）声速：')") READ(5,*) A_FAR !WRITE(*,"('请输入来流压力：')") READ(5,*) P_FAR !WRITE(*,"('请输入平板的长度：')") READ(5,*) LHORI ALLOCATE(VX(IMAX,JMAX),VY(IMAX,JMAX),A(IMAX,JMAX),ROU(IMAX,JMAX),P(I MAX,JMAX),E(IMAX,JMAX),TEMP(IMAX,JMAX),MIU(IMAX,JMAX),K(IMAX,JMAX),& U1(IMAX,JMAX),U2(IMAX,JMAX),U3(IMAX,JMAX),U5(IMAX,JMAX),E1(IMAX,JMAX),E 2(IMAX,JMAX),E3(IMAX,JMAX),E5(IMAX,JMAX),F1(IMAX,JMAX),& F2(IMAX,JMAX),F3(IMAX,JMAX),F5(IMAX,JMAX),ROU_PRE(IMAX,JMAX)) TW=TW0*TEMP_FAR MIU_FAR=DYNVIS(TEMP_FAR) !Sutherlan 定律 *****'/,8x,'IMAX=',I4/,8X,'JMAX=',I4/) 输入的网格信息
    要 参 数 来 流 主 计 主 要 算 域 FORMAT(5X,'***** FORMAT(5X,'***** ROU_FAR=P_FAR/(R*TEMP_FAR) V_FAR=MACH_FAR*A_FAR REL=ROU_FAR*V_FAR*LHORI/MIU_FAR E_FAR=CV*TEMP_FAR K_FAR=THERMC(MIU_FAR) !根据普朗特常数的公式 DELTA=5*LHORI/SQRT(REL) LVERT=5*DELTA DX=LHORI/(IMAX-1) DY=LVERT/(JMAX-1) WRITE(6,103) MACH_FAR,A_FAR,V_FAR,TEMP_FAR,P_FAR,ROU_FAR,MIU_FAR,K_FAR,E_FAR 103 数 *****'/,8x,'MACH_FAR=',F16.12/,8X,'A_FAR=',F16.12/,8X,'V_FAR=',F16.8/,8X,'TEMP_FAR=', F16.12/,8X,& 'P_FAR=',F16.4/,8X,'ROU_FAR=',F16.12/,8X,'MIU_FAR=',F16.12/,8X,'K_FAR=',F16.12/,8X,'E_ FAR=',F16.8/) WRITE(6,104) LHORI,LVERT,DX,DY 104 参 *****'/,8x,'LHORT=',F16.12/,8X,'LVERT=',F16.12/,8X,'DX=',F16.12/,8X,'DY=',F16.12/) DO ITER=1,MAXIT IF(ITER.EQ.1) THEN !初始化流场所有参数操作 DO J=1,JMAX DO I=1,IMAX IF(J.EQ.1) THEN VX(I,J)=0;VY(I,J)=0;TEMP(I,J)=TW ELSE VX(I,J)=V_FAR;VY(I,J)=0;TEMP(I,J)=TEMP_FAR END IF ROU(I,J)=ROU_FAR;P(I,J)=P_FAR;E(I,J)=E_FAR;MIU(I,J)=MIU_FAR;K(I,J)=K_FAR;A(I,J)=( GAMMA*R*TEMP(I,J))**0.5 END DO END DO END IF DO I=1,IMAX DO J=1,JMAX A(I,J)=(GAMMA*R*TEMP(I,J))**0.5 END DO END DO CALL TSTEP(VX,VY,A,MIU,ROU,DX,DY,IMAX,JMAX,DT) !计算 DT 的子程序 CALL MAC(VX,VY,ROU,P,E,TEMP,MIU,K,U1,U2,U3,U5,E1,E2,E3,E5,F1,F2,F3,F5,IMAX,JMAX,D X,DY,DT,TEMP_FAR,P_FAR,V_FAR,TW,ROU_PRE) !调用 MacCormack 方法的子程序 IF(CONVER(ITER,MAXIT,IMAX,JMAX,ROU_PRE,ROU).EQ.1) EXIT !检验收敛的子程序 END DO
    CALL MODT(VX,ROU,DY,IMAX,JMAX) !确认连续性的子程序 CALL OUTPUT(VX,VY,ROU,P,E,TEMP,MIU,K,IMAX,JMAX,P_FAR,V_FAR,TEMP_FAR,A_FAR) 生成用以绘图的数据文件子程序 STOP END PROGRAM MAIN ! !TSTEP(计算时间步长的子程序) SUBROUTINE TSTEP(VX,VY,A,MIU,ROU,DX,DY,IMAX,JMAX,DT) USE GLOBAL IMPLICIT NONE INTEGER::I,J,IMAX,JMAX REAL::VX(IMAX,JMAX),VY(IMAX,JMAX),A(IMAX,JMAX),MIU(IMAX,JMAX),ROU(IMA X,JMAX),VPIE(IMAX,JMAX),DT_CFL(IMAX,JMAX),VPIE_MAX,DX,DY,DT REAL,PARAMETER::K_COU=0.6 !为保证解的稳定性而给出的"虚拟因子" REAL::REX(IMAX,JMAX),REY(IMAX,JMAX),REXMAX,REYMAX VPIE_MAX=0;DT=100.;VPIE=0;DT_CFL=0;REXMAX=0;REYMAX=0 !初始化参数，为了避 免 DT 初值对后续产生不好的影响，预设一个很大的 DT 值 DO I=2,IMAX-1 DO J=2,JMAX-1 REX(I,J)=ROU(I,J)*VX(I,J)*DX/MIU(I,J) REY(I,J)=ROU(I,J)*VY(I,J)*DY/MIU(I,J) VPIE(I,J)=4./3.*GAMMA*MIU(I,J)**2/(PR*ROU(I,J)) IF(VPIE(I,J).GT.VPIE_MAX) THEN VPIE_MAX=VPIE(I,J) END IF IF(REXMAX.LT.REX(I,J)) THEN REXMAX=REX(I,J) END IF IF(REYMAX.LT.REY(I,J)) THEN REYMAX=REY(I,J) END IF END DO END DO WRITE(6,97) REXMAX,REYMAX 97 FORMAT(5X,'REXMAX:',F12.6,5X,'REYMAX:',F20.18/) IF(REXMAX.GT.40.OR.REYMAX.GT.4) THEN WRITE(6,98) 98 FORMAT(5X,'步长选取出错') !STOP END IF DO I=2,IMAX-1 DO J=2,JMAX-1
    DT_CFL(I,J)=(ABS(VX(I,J))/DX+ABS(VY(I,J))/DY+A(I,J)*SQRT(1./(DX*DX)+1./(DY*DY))+ 2*VPIE_MAX*(1/(DX*DX)+1/(DY*DY)))**(-1.) IF(DT_CFL(I,J).LT.DT) THEN DT=DT_CFL(I,J) END IF END DO END DO DT=K_COU*DT WRITE(6,99) DT 99 FORMAT(5X,'DT',F30.20/) RETURN END !MAC(调用 MacCormack 方法的子程序) SUBROUTINE MAC(VX,VY,ROU,P,E,TEMP,MIU,K,U1,U2,U3,U5,E1,E2,E3,E5,F1,F2,F3,F5,IMAX,JMAX,D X,DY,DT,TEMP_FAR,P_FAR,V_FAR,TW,ROU_PRE) USE GLOBAL IMPLICIT NONE INTEGER::IE,I,J,IMAX,JMAX REAL::VX(IMAX,JMAX),VY(IMAX,JMAX),ROU(IMAX,JMAX),P(IMAX,JMAX),E(IMAX,J MAX),ET(IMAX,JMAX),TEMP(IMAX,JMAX),MIU(IMAX,JMAX),K(IMAX,JMAX),U1(IMA X,JMAX),U2(IMAX,JMAX),U3(IMAX,JMAX),U5(IMAX,JMAX),& U1_ES(IMAX,JMAX),U2_ES(IMAX,JMAX),U3_ES(IMAX,JMAX),U5_ES(IMAX,JMAX),E1( IMAX,JMAX),E2(IMAX,JMAX),E3(IMAX,JMAX),E5(IMAX,JMAX),F1(IMAX,JMAX),F2(IM AX,JMAX),F3(IMAX,JMAX),F5(IMAX,JMAX),ROU_PRE(IMAX,JMAX),& DX,DY,DT,TEMP_FAR,P_FAR,V_FAR,TW REAL,EXTERNAL::DYNVIS,THERMC,TAUXX,TAUXY,TAUYY,QX,QY DYNVIS,THERMC,TAUXX,TAUXY,TAUYY,QX,QY 是函数，不是变量 DO I=1,IMAX !将原始变量转换为 U 通量 DO J=1,JMAX ET(I,J)=ROU(I,J)*(E(I,J)+((VX(I,J))**2+(VY(I,J))**2)/2.) U1(I,J)=ROU(I,J) U2(I,J)=ROU(I,J)*VX(I,J) U3(I,J)=ROU(I,J)*VY(I,J) U5(I,J)=ET(I,J) ROU_PRE(I,J)=ROU(I,J) !存储当前时刻每一网格点的的密度值 END DO END DO DO I=1,IMAX !计算 t 时刻的 E/F 列向量 DO J=1,JMAX IE=1 !代表当前调用的是求预测步 E 向量中粘性应力项及热传导项的 CASE E1(I,J)=ROU(I,J)*VX(I,J) ! 声 明
    E2(I,J)=ROU(I,J)*VX(I,J)**2+P(I,J)-TAUXX(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) E3(I,J)=ROU(I,J)*VX(I,J)*VY(I,J)-TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) E5(I,J)=(ET(I,J)+P(I,J))*VX(I,J)-VX(I,J)*TAUXX(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)-V Y(I,J)*TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)+QX(TEMP,K,DX,DY,IMAX,JMAX,I, J,IE) IE=2 !代表当前调用的是求预测步 F 向量中粘性应力项及热传导项的 CASE F1(I,J)=ROU(I,J)*VY(I,J) F2(I,J)=ROU(I,J)*VX(I,J)*VY(I,J)-TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) F3(I,J)=ROU(I,J)*VY(I,J)**2+P(I,J)-TAUYY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) F5(I,J)=(ET(I,J)+P(I,J))*VY(I,J)-VX(I,J)*TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)-V Y(I,J)*TAUYY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)+QY(TEMP,K,DX,DY,IMAX,JMAX,I, J,IE) END DO END DO DO I=2,IMAX-1 !预测 T+DT 时刻的 U1/U2/U3/U5 DO J=2,JMAX-1 U1_ES(I,J)=U1(I,J)-(E1(I+1,J)-E1(I,J))/DX*DT-(F1(I,J+1)-F1(I,J))/DY*DT U2_ES(I,J)=U2(I,J)-(E2(I+1,J)-E2(I,J))/DX*DT-(F2(I,J+1)-F2(I,J))/DY*DT U3_ES(I,J)=U3(I,J)-(E3(I+1,J)-E3(I,J))/DX*DT-(F3(I,J+1)-F3(I,J))/DY*DT U5_ES(I,J)=U5(I,J)-(E5(I+1,J)-E5(I,J))/DX*DT-(F5(I,J+1)-F5(I,J))/DY*DT END DO END DO DO I=2,IMAX-1 !将通量预测值解耦 DO J=2,JMAX-1 CALL DECODE(U1_ES(I,J),U2_ES(I,J),U3_ES(I,J),U5_ES(I,J),ROU(I,J),VX(I,J),VY(I,J),ET(I,J),E(I,J), TEMP(I,J),P(I,J)) !注意，此时原始变量 VX,VY,ROU,ET,E,TEMP,P 均已经更新为 T+DT 时刻 的预测值 MIU(I,J)=DYNVIS(TEMP(I,J)) !注意，此时原始变量 MIU 已更新为 T+DT 时刻的预测值 K(I,J)=THERMC(MIU(I,J)) !注意，此时原始变量 K 已更新为 T+DT 时刻的预测值 END DO END DO CALL BC(VX,VY,TEMP,P,ROU,E,ET,MIU,K,V_FAR,TEMP_FAR,P_FAR,TW,IMAX,JMAX) ! 预测步设置边界条件 DO I=1,IMAX !计算 t+dt 时刻的预测 E/F 列向量 DO J=1,JMAX IE=3 !代表当前调用的是求校正步 E 向量中粘性应力项及热传导项的 CASE E1(I,J)=ROU(I,J)*VX(I,J) !注意，此处及以下中的 E/F 向量中每一项均已经更新为 t+dt 时刻 的预测值 E2(I,J)=ROU(I,J)*VX(I,J)**2+P(I,J)-TAUXX(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) E3(I,J)=ROU(I,J)*VX(I,J)*VY(I,J)-TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) E5(I,J)=(ET(I,J)+P(I,J))*VX(I,J)-VX(I,J)*TAUXX(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)-V Y(I,J)*TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)+QX(TEMP,K,DX,DY,IMAX,JMAX,I, J,IE)
    IE=4 !代表当前调用的是求校正步 F 向量中粘性应力项及热传导项的 CASE F1(I,J)=ROU(I,J)*VY(I,J) F2(I,J)=ROU(I,J)*VX(I,J)*VY(I,J)-TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) F3(I,J)=ROU(I,J)*VY(I,J)**2+P(I,J)-TAUYY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE) F5(I,J)=(ET(I,J)+P(I,J))*VY(I,J)-VX(I,J)*TAUXY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)-V Y(I,J)*TAUYY(VX,VY,MIU,DX,DY,IMAX,JMAX,I,J,IE)+QY(TEMP,K,DX,DY,IMAX,JMAX,I, J,IE) END DO END DO !WRITE(*,"('NO PROBLEM')") DO I=2,IMAX-1 !校正步，对所有内点进行校正，得到校正后的 T+DT 时刻的 U1/U2/U3/U5 DO J=2,JMAX-1 U1(I,J)=0.5*(U1(I,J)+U1_ES(I,J)-(E1(I,J)-E1(I-1,J))/DX*DT-(F1(I,J)-F1(I,J-1))/DY*DT) U2(I,J)=0.5*(U2(I,J)+U2_ES(I,J)-(E2(I,J)-E2(I-1,J))/DX*DT-(F2(I,J)-F2(I,J-1))/DY*DT) U3(I,J)=0.5*(U3(I,J)+U3_ES(I,J)-(E3(I,J)-E3(I-1,J))/DX*DT-(F3(I,J)-F3(I,J-1))/DY*DT) U5(I,J)=0.5*(U5(I,J)+U5_ES(I,J)-(E5(I,J)-E5(I-1,J))/DX*DT-(F5(I,J)-F5(I,J-1))/DY*DT) END DO END DO DO I=2,IMAX-1 !将通量预测值解耦 DO J=2,JMAX-1 CALL DECODE(U1(I,J),U2(I,J),U3(I,J),U5(I,J),ROU(I,J),VX(I,J),VY(I,J),ET(I,J),E(I,J),TEMP(I,J),P(I,J) ) !注意，此时原始变量 VX,VY,ROU,ET,E,TEMP,P 均已经更新为 T+DT 时刻的校正值 MIU(I,J)=DYNVIS(TEMP(I,J)) !注意，此时原始变量 MIU 已更新为 T+DT 时刻的校正值 K(I,J)=THERMC(MIU(I,J)) !注意，此时原始变量 K 已更新为 T+DT 时刻的校正值 END DO END DO CALL BC(VX,VY,TEMP,P,ROU,E,ET,MIU,K,V_FAR,TEMP_FAR,P_FAR,TW,IMAX,JMAX) ! 校正步设置边界条件 END !将 U 通量解耦为原始变量的子程序（DECODE） SUBROUTINE DECODE(A1,A2,A3,A5,ROU,VX,VY,ET,E,TEMP,P) USE GLOBAL IMPLICIT NONE REAL::A1,A2,A3,A5,ROU,VX,VY,ET,E,TEMP,P ROU=A1 VX=A2/A1 VY=A3/A1 ET=A5 E=A5/A1-(VX**2+VY**2)/2. TEMP=E/CV P=ROU*R*TEMP
    RETURN END !设置边界条件的子程序 BC SUBROUTINE BC(VX,VY,TEMP,P,ROU,E,ET,MIU,K,V_FAR,TEMP_FAR,P_FAR,TW,IMAX,JMAX) USE GLOBAL IMPLICIT NONE INTEGER::I,J,IMAX,JMAX REAL::VX(IMAX,JMAX),VY(IMAX,JMAX),TEMP(IMAX,JMAX),P(IMAX,JMAX),ROU(IM AX,JMAX),E(IMAX,JMAX),ET(IMAX,JMAX),MIU(IMAX,JMAX),K(IMAX,JMAX),& V_FAR,TEMP_FAR,P_FAR,TW REAL,EXTERNAL::DYNVIS,THERMC VX(1,1)=0;VY(1,1)=0 !前缘点无滑移边界条件 TEMP(1,1)=TEMP_FAR;P(1,1)=P_FAR !前缘点温度和压力分别为自由来流值 ROU(1,1)=P(1,1)/(R*TEMP(1,1)) E(1,1)=CV*TEMP(1,1) ET(1,1)=ROU(1,1)*(E(1,1)+((VX(1,1))**2+(VY(1,1))**2)/2) MIU(1,1)=DYNVIS(TEMP(1,1)) K(1,1)=THERMC(MIU(1,1)) DO J=2,JMAX !左边界 VX(1,J)=V_FAR !左边界处除前缘点外 x 方向速度为来流速度 VY(1,J)=0 !左边界处 y 方向速度为 0 TEMP(1,J)=TEMP_FAR !左边界处温度为来流自由值 P(1,J)=P_FAR !左边界处压力为来流自由值 ROU(1,J)=P(1,J)/(R*TEMP(1,J)) E(1,J)=CV*TEMP(1,J) ET(1,J)=ROU(1,J)*(E(1,J)+((VX(1,J))**2+(VY(1,J))**2)/2) MIU(1,J)=DYNVIS(TEMP(1,J)) K(1,J)=THERMC(MIU(1,J)) END DO DO I=1,IMAX !上边界 VX(I,JMAX)=V_FAR !上边界处 x 方向速度为来流速度 VY(I,JMAX)=0 !上边界处 y 方向速度为 0 TEMP(I,JMAX)=TEMP_FAR !上边界处温度为来流自由值 P(I,JMAX)=P_FAR !上边界处压力为来流自由值 ROU(I,JMAX)=P(I,JMAX)/(R*TEMP(I,JMAX)) E(I,JMAX)=CV*TEMP(I,JMAX) ET(I,JMAX)=ROU(I,JMAX)*(E(I,JMAX)+((VX(I,JMAX))**2+(VY(I,JMAX))**2)/2.) MIU(I,JMAX)=DYNVIS(TEMP(I,JMAX)) K(I,JMAX)=THERMC(MIU(I,JMAX)) END DO DO I=2,IMAX !平板边界
    分享到：
    赞 收藏

    相关推荐

    资料库

    课程资源

    共收录17145份资料，累计13个分类，关注成员有19位，主要包括：PHP,网络管理,网页制作,Java,.Net,数据库,3G/移动开发,C/C++,游戏开发,嵌入式,讲义,软件测试,专业指导

    热门标签

    PHP 网络管理 网页制作 Java .Net 数据库 3G/移动开发 C/C++ 游戏开发 嵌入式 讲义 软件测试 专业指导

    最新资料