| 查看: 264 | 回复: 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人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
物理学I论文润色/翻译怎么收费?
已经有265人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有11人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有13人回复
-
桂林理工大学物理学专业招收调剂
已经有13人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
美国顶级物理期刊《应用物理快报》（APL）的编辑欺骗和歧视作者及AIP的官僚主义傲慢
已经有6人回复
-
津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生一名
已经有0人回复
找到一些相关的精华帖子，希望有用哦~
请教 fortran 运行错误的原因
已经有13人回复
- 求教fortran默认值的问题
已经有8人回复
- 求高人帮忙看一下我这个Fortran程序的错误出在哪里？
已经有12人回复
- 物理模拟是用c语言好还是fortran好？
已经有30人回复
- fortran的一个错误
已经有9人回复
- 怎么使用fortran程序创建一个立方晶格模型
已经有5人回复
- 【求助】读文本数据出现段错误[fortran]
已经有14人回复
- 【求助】请教如何在超算上运行fortran程序
已经有4人回复
- 【求助】此FORTRAN90程序没有错误，但是运行不出来。
已经有4人回复
- 【求助】用fortran求解大型线性方程组时出现的错误【已解决】
已经有11人回复
- 【求助】请问有拟合分子力场参数的FORTRAN程序吗？
已经有13人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
5