24小时热门版块排行榜

返回列表

zeng0927

金虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1022
帖子: 6
在线: 61.5小时
虫号: 1349199
注册: 2011-07-18
性别: GG
专业: 核技术及其应用

[求助] 求教一个FORTRAN一维稳态流体模拟的程序错误

PROGRAM MAIN
! Predictor Corrector Method for Coupled Exchange in Column
! Boundary conditions C(N+1)=C(N), C(0)=CHS or CDS
! IMPLICIT REAL*8 (A-H),REAL*8 (P-Z)
DIMENSION CHG(0:200),CHGM(0:200)
DIMENSION CDG(0:200),CDGM(0:200)
DIMENSION PA(0:200),VF(0:200),AVF(0:200)
REAL*8 K,M,NTIME,L
OPEN(UNIT=3,FILE='TEST1.OUT',FORM='FORMATTED',STATUS='UNKNOWN')
OPEN(UNIT=4,FILE='INPUT.TXT',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(*,*)'Input column lenghth (cm)'
READ(4,*)L
WRITE(*,*)'Input temperature (K)'
READ(4,*)T
WRITE(*,*)'Input let in pressure (Torr)'
READ(4,*)PIN
WRITE(*,*)'Input let out pressure (Torr)'
READ(4,*)POUT
WRITE(*,*)'Input porosity (1)'
READ(4,*)E
WRITE(*,*)'Input viscosity (10-6pa.s)'
READ(4,*)M
WRITE(*,*)'Input permeability (10-12m2)'
READ(4,*)K
WRITE(*,*)'Input total exchange time(sec)'
READ(4,*)TIME
L=L*0.01
!M=M*0.000001
!K=K*0.000000000001
N=199
NPO=N+1
DELL=L/(NPO)
DELT=0.1
DELTDT=DELT/2
NTSTOP=TIME/DELT
CHGR=1
CDGR=0
DO 10 J=1,NPO
CHG(J)=0
10 CDG(J)=0
CHG(0)=CHGR
CDG(0)=CDGR
DO 12 NT=1,NTSTOP
! Perform Predictor Operation (half time step)
CALL COEFFV(PA,VF,AVF,PIN,POUT,DELL,DELTDT,L,K,E,M,N)
CALL UPTWOHG(CHG,CHG,CHGM,AVF,DELTDT,E,N)
CALL UPTWODG(CDG,CDG,CDGM,AVF,DELTDT,E,N)
! Perform Corrector Operation (full time step)
CALL COEFFV(PA,VF,AVF,PIN,POUT,DELL,DELT,L,K,E,M,N)
CALL UPTWOHG(CHG,CHGM,CHGM,AVF,DELT,E,N)
CALL UPTWODG(CDG,CDGM,CDGM,AVF,DELT,E,N)
DO 11 J=0,NPO
CHGM(J)=(CHGM(J)+CHG(J))/2
11 CDGM(J)=(CDGM(J)+CDG(J))/2
CALL COEFFV(PA,VF,AVF,PIN,POUT,DELL,DELT,L,K,E,M,N)
CALL UPTWOHG(CHG,CHGM,CHG,AVF,DELT,E,N)
CALL UPTWODG(CDG,CDGM,CDG,AVF,DELT,E,N)
NTIME=NT*DELT
WRITE(3,202) NTIME,' ',CHG(1),' ',CHG(10),' ',CHG(100),' ',CHG(150)
12 CONTINUE
202 FORMAT (1F7.2,A1,1F12.5,A1,1F12.5,A1,1F12.5,A1,1F12.5)
STOP
END

SUBROUTINE COEFFV(PA,VF,AVF,PIN,POUT,DELL,DELT,L,K,E,M,N)
DIMENSION PA(0:200),VF(0:200),AVF(0:200)
REAL*8 K,M,L
DO 10 J=1,N
10 PA(J)=(PIN**2-((PIN**2)-(POUT**2))*(J*DELL/L))**0.5
PA(0)=PIN
PA(N+1)=POUT
DO 11 J=1,N
11 VF(J)=(PA(J-1)-PA(J+1))*K/(2*DELL*E*M)
DO 12 J=1,N
12 AVF(J)=VF(J)*DELT/(4*DELL)
AVF(0)=AVF(1)
AVF(N+1)=AVF(N)
RETURN
END

SUBROUTINE UPTWOHG(CHG,CHGM,CHGN,AVF,DELT,E,N)
! This is version for gas phase H concentration change in column coordinates
! IMPLICIT REAL*8 (A-H),REAL*8 (P-Z)
DIMENSION CHG(0:200),CHGM(0:200),CHGN(0:200)
DIMENSION AVF(0:200)
DIMENSION A(199),B(199),C(199),QH(199),UDAG(199)
DO 10 J=1,N
A(J)=-AVF(J)
B(J)=1
10 C(J)=AVF(J)
! Calculate H array at beginning of time step
DO 11 J=1,N
QH(J)=CHG(J)+AVF(J)*(CHGM(J-1)-CHGM(J+1))+2*(AVF(J-1)-AVF(J+1))*CHG(J)
11 CONTINUE
QH(N)=QH(N)+AVF(N)*CHGM(N-1)
CALL TRIDAG (A,B,C,QH,UDAG,N)
DO 12 J=1,N
12 CHGN(J)=UDAG(J)
CHGN(0)=CHG(0)
CHGN(N+1)=CHG(N)
RETURN
END

SUBROUTINE UPTWODG(CDG,CDGM,CDGN,AVF,DELT,E,N)
! This is version for gas phase D concentration change in column coordinates
! IMPLICIT REAL*8 (A-H),REAL*8 (P-Z)
DIMENSION CDG(0:200),CDGM(0:200),CDGN(0:200)
DIMENSION AVF(0:200)
DIMENSION A(199),B(199),C(199),QD(199),UDAG(199)
DO 10 J=1,N
A(J)=-AVF(J)
B(J)=1
10 C(J)=AVF(J)
! Calculate D array at beginning of time step
DO 11 J=1,N
QD(J)=CDG(J)+AVF(J)*(CDGM(J-1)-CDGM(J+1))+2*(AVF(J-1)-AVF(J+1))*CDG(J)
11 CONTINUE
CALL TRIDAG (A,B,C,QD,UDAG,N)
DO 12 J=1,N
12 CDGN(J)=UDAG(J)
CDGN(0)=CDG(0)
CDGN(N+1)=CDG(N)
RETURN
END

SUBROUTINE TRIDAG (A,B,C,Q,U,N)
! This is from Numerical Recipes in Fortran
! IMPLICIT REAL*8 (A-H),REAL*8 (P-Z)
PARAMETER (NMAX=200)
DIMENSION GAM(NMAX),A(N),B(N),C(N),Q(N),U(N)
! IF(B(1).EQ.0.)PAUSE
BET=B(1)
U(1)=Q(1)/BET
DO 11 J=2,N
GAM(J)=C(J-1)/BET
BET=B(J)-A(J)*GAM(J)
! IF(BET.EQ.0.)PAUSE
U(J)=(Q(J)-A(J)*U(J-1))/BET
11 CONTINUE
DO 12 J=N-1,1,-1
U(J)=U(J)-GAM(J+1)*U(J+1)
12 CONTINUE
RETURN
END

具体的问题是：
速度VF是线性递增的，我估计是有压力降的原因，不知道对不对；

还有就是SUBROUTINE UPTWOHG这个子程序算出的CHG与理论值相差太大，一直找不到原因。

涉及的公式

回复此楼