| 查看: 298 | 回复: 0 | |||
[求助]
求教一个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与理论值相差太大，一直找不到原因。  涉及的公式 |
回复此楼» 猜你喜欢
如何从铁电相到顺电相。
已经有1人回复
-
各位大佬，求一份最新Wien2k版本作为学习用途
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有214人回复
-
南科大活性流体和软物质课题组诚招2027级博士、硕士研究生和博后
已经有42人回复
-
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有24人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
PVA溶解
已经有2人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
找到一些相关的精华帖子，希望有用哦~
请教 fortran 运行错误的原因
已经有13人回复
- 求教fortran默认值的问题
已经有8人回复
- 求高人帮忙看一下我这个Fortran程序的错误出在哪里？
已经有12人回复
- 物理模拟是用c语言好还是fortran好？
已经有30人回复
- fortran的一个错误
已经有9人回复
- 怎么使用fortran程序创建一个立方晶格模型
已经有5人回复
- 【求助】读文本数据出现段错误[fortran]
已经有14人回复
- 【求助】请教如何在超算上运行fortran程序
已经有4人回复
- 【求助】此FORTRAN90程序没有错误，但是运行不出来。
已经有4人回复
- 【求助】用fortran求解大型线性方程组时出现的错误【已解决】
已经有11人回复
- 【求助】请问有拟合分子力场参数的FORTRAN程序吗？
已经有13人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
