| 查看: 511 | 回复: 0 | |||
[交流]
【求助】求大侠看看这个程序怎么修改
|
|
PROGRAM EMP CHARACTER*10 FIL1,KWORD DIMENSION A(10,500),X(10,3) !定义两个数组 REAL X !TNT的参数 VD=6930 !m/s DEN=1.630 !g/cm3 ! PH=21*10E9 !pa ! write(*,*)' AN or other EXPLOSIVE? (AN=1)' ! r=(vd*vd/1608/1608+1.)**.5 !铵油多方指数的计算 !DEN,TIE分别为炸药的密度和厚度 DEN=1.63 TIE=20 !DEP,TIP分别为平板的密度和厚度 DEP=1.7 TIP=8 !Q表示质量比 Q=(DEN*TIE)/DEP/TIP !R即为等熵指数K 单位分别为kg/m3,m/s,pa !或者为g/cm3,km/s,Gpa R=2.72 !R=DEN0*VD**2/PH-1.0 !角度 CTA=30 !插值点数 N=20 EPS=0.001 WRITE(*,3)Q,R,TIE !Q质量比R多方指数及GUMMA WRITE(*,7)VD,DEN !炸药的爆速和密度 7 FORMAT(6H Vd=,F7.2,5H(m/s),10H DENSITYe=,F5.3,9H(g/cm**3)) 3 FORMAT(4H R=,F7.4,9H GAMMA=,F5.3,22H Explosive thickness, & F8.4,4H(mm)) WRITE(*,5) 5 FORMAT(' Angle(DEG) Vp(m/s) H(mm) P(GPa) X(mm)') DEM=DEN/2. SS=1.0172 ! CALL MMM(DEM,VD,R,EPS,CTA,SS,CTAB) !CTA(指的是角度)子例行程序 CALL PMFRI(N,2.2292,4.2945,R,A,0.,0.,1.) !务必注意这两步的作用 DO 10 I=N,1,-1 DO 10 J=1,10 10 A(J,I+1)=A(J,I) DO 20 J=1,10 20 A(J,1)=0. A(7,1)=3.1415926/2. DD=-180./3.1415926*A(1,1) !角度转化 DV=-VD*2.*SIN(A(1,1)/2.) !泰勒公式计算的抛掷速度 DHH=-TIE*A(9,1) !炸药的厚度乘以角度 PH=DEN*(VD/1000.)**2/(R+1.) !爆轰波头的压力 WRITE(*,4)DD,DV,DHH,PH,A(8,1)*Tie 4 FORMAT(10(1X,F9.4)) 200 DO 30 J=1,10 X(J,1)=A(J,1) 30 X(J,2)=A(J,2) CALL RIM(X,R,-1,Q,EPS,P) !调用边界插值函数 DO 40 I=1,10 40 A(I,1)=X(I,3) DD=-180./3.1415926*A(1,1) DV=-VD*2.*SIN(A(1,1)/2.) DHH=-TIE*A(9,1) ! WRITE(*,4)P P的数值 WRITE(*,4)DD,DV,DHH,P*PH,A(8,1)*Tie,A(1,1),A(4,1),A(2,1),A(6,1) N=N-1 DO 50 I=1,N DO 60 J=1,10 X(J,1)=A(J,I) 60 X(J,2)=A(J,I+2) CALL INSIDE(X,EPS) !调用内部插值函数 DO 70 J=1,10 70 A(J,I+1)=X(J,3) 50 CONTINUE IF(N.EQ.1) GOTO 100 GOTO 200 100 WRITE(*,*)' Anther Calculation(Y/N)?' ! READ(*,'(A)')KWORD ! IF(KWORD.EQ.'Y') GOTO 909 ! WRITE(*,77)FIL1 STOP END C plane-fluid subrountine !******************************* !等熵的压力关系 !******************************* FUNCTION PS(X)!密度我变量 V=1.63/X A=373.8 !A,B,C GPA B=2.747 CC=0.734 R1=4.15 R2=0.90 W=0.3 PS=(A*EXP(-R1*V)+B*EXP(-R2*V)+CC/(V**(W+1)))*1.0 RETURN END !******************************* !等熵的比内能函数 !******************************* FUNCTION ES(X) !X指的是密度 V=1.63/X A=373.8 !A,B,C GPA B=2.747 CC=0.734 R1=4.15 R2=0.90 W=0.3 ES=A*EXP(-R1*V)/R1+B*EXP(-R2*V)/R2+CC/(V**W)/W RETURN END !******************************* !等熵的声速函数 !******************************* !X指的是密度 FUNCTION CS(X)! 声速的函数 V=1.63/X A=373.8 !A,B,C GPA B=2.747 CC=0.734 R1=4.15 R2=0.90 W=0.3 DPV=-A*R1*EXP(-R1*V)-B*R2*EXP(-R2*V)-(W+1)*CC/V**(W+1) CS2=-V**2/X*DPV CS=SQRT(CS2) RETURN END !******************************* !普朗特-迈耶函数 !******************************* !X指的是密度 FUNCTION F(X) !普朗特函数需要修改CS时含有变量的函数 ! XI=ES(X)+PS(X)/X ! XI0=ES(DENH)+PS(DENH)/DENH+0.5*qH**2 XX=LOG(SQRT(30-XI(X)))!XI0需要带入 F=-SQRT((2.0*EXP(2.0*XX)-CS(X)**2)/CS(X)**2) END !******************************* !定义了焓的函数 !******************************* FUNCTION XI(X) XI=ES(X)+PS(X)/X END !******************************* !辛普生积分函数 !******************************* FUNCTION SIMPSON(F,A,B,N) H=(B-A)/(2.0*N) S=F(A)-F(B) DO I=1,N S1=F(A+(2.0*I-1)*H) S2=F(A+2.0*I*H) S=S+4.0*S1+2.0*S2 END DO SIMPSON=H/3.0*S END !********************************* !利用辛普生积分求得马赫角 !********************************* FUNCTION V(X) V=SIMPSON(F(X),30.0,XI(X),20) RETURN END !******************************* !普迈函数的反函数 !******************************* FUNCTION VN(Q,EPS) !反函数!Q角度R等熵指数需要修改 Y0=0. Y1=2. IF(V(Y1).LE.Q)THEN Y0=Y1 Y1=5. END IF IF(V(Y1).LE.Q)THEN Y0=Y1 Y1=10. END IF IF(V(Y1).LE.Q)THEN Y0=Y1 Y1=15. END IF IF(V(Y1).LE.Q)THEN Y0=Y1 Y1=20. END IF IF(V(Y1).LE.Q)THEN Y0=Y1 Y1=30. END IF 1 YY=(Y0+Y1)/2. IF(V(YY).LE.Q)THEN Y0=YY ELSE Y1=YY !VN得出来的是a=log(i0-i) END IF IF(ABS(Y1-Y0).LT.(.01*EPS)) GOTO 2 GOTO 1 2 VN=YY RETURN END !******************************* !正弦函数的反函数 !******************************* FUNCTION ASINN(X) ASINN=ATAN(1./X) RETURN END !******************************* !焓的反函数 !******************************* FUNCTION VI(Q,EPS) Y0=0. Y1=2. IF(XI(Y1).LE.Q)THEN Y0=Y1 Y1=5. END IF IF(XI(Y1).LE.Q)THEN Y0=Y1 Y1=10. END IF IF(XI(Y1).LE.Q)THEN Y0=Y1 Y1=15. END IF IF(XI(Y1).LE.Q)THEN Y0=Y1 Y1=20. END IF IF(XI(Y1).LE.Q)THEN Y0=Y1 Y1=30. END IF 3 YY=(Y0+Y1)/2. IF(XI(YY).LE.Q)THEN Y0=YY ELSE Y1=YY !VN得出来的是a=log(i0-i) END IF IF(ABS(Y1-Y0).LT.(.01*EPS)) GOTO 4 GOTO 3 4 VI=YY RETURN END !******************************* !内部插值函数 已修改 还需要修改焓0，VN !******************************* SUBROUTINE INSIDE(X,EPS) DIMENSION X(10,3) X(1,3)=(X(10,2)-X(10,1)+X(1,1)+X(1,2))/2. X(10,3)=X(1,3)-X(1,1)+X(10,1) X(5,3)=VN(X(10,3),EPS) X(3,3)=30.0-EXP(2.0*X(5,3))!X(3,0)指的是I0 X(2,3)=VI(X(3,3),EPS) X(6,3)=SQRT(2.0*EXP(2.0*X(5,3))/CS(X(2,3))**2) X(7,3)=ASINN(X(6,3)) A=1./TAN(X(1,1)+X(7,1)) B=1./TAN(X(1,2)-X(7,2)) X(9,3)=(X(8,2)-X(8,1)+A*X(4,1)-B*X(9,2))/(A-B) X(8,3)=X(8,1)+A*(X(9,3)-X(9,1)) RETURN END !******************************* !边界插值函数 !******************************* SUBROUTINE RIM(X,R,K,Q,EPS,PRESS) DIMENSION X(10,3) A=1./TAN(X(1,2)+K*X(7,2)) F=X(1,2)+K*X(7,2) PK1=P(X(2,1))*Q/(R+1.) PK=PK1 DO 10 I=1,2 SINF=SIN(X(1,1)-F) X(1,3)=X(1,1)-K*(X(3,2)-X(3,1)-A*(X(4,2)-X(4,1)))*SIN(F)/SINF & *PK X(10,3)=FLOAT(K)*(X(1,3)-X(1,1))+X(10,2) X(5,3)=VN(X(10,3),EPS) X(3,3)=X(3,0)-EXP(2.0*X(5,3)) ! X(2,3)=VI(X(3,3)) !VI是焓的反函数 X(6,3)=SQRT(2.0*EXP(2.0*X(5,3))/(CS(X(2,3))**2)) X(7,3)=ASINN(X(6,3)) PK3=PK1+(X(3,3)-X(3,1))/0.5/(X(2,1)+X(2,3))*Q/(R+1.) PK=.5*PK1+.5*PK3 10 CONTINUE X(8,3)=(SIN(X(1,3))-SIN(X(1,1)))*K/PK+X(8,1) X(9,3)=(COS(X(1,1))-COS(X(1,3)))*K/PK+X(9,1) PRESS=PK/(Q/(R+1.))!计算P的表达式 RETURN END !******************************* !密度流初值选取 已修改 !******************************* SUBROUTINE PDENFRI(N,DEN0,DENM,R,X,Q0,X0,Y0) DIMENSION X(10,N) DO 1 I=1,N X(2,I)=DEN0+(DENM-DEN0)/(N-1)*(I-1) IF(X(2,I).LE..1E-6)THEN X(2,I)=DEN0+0.1E-6 END IF X(3,I)=ES(X(2,I))+PS(X(2,I))/X(2,I) ! X(4,I)=VI(X(3,I)) X(5,I)=LOG(SQRT(30.0-X(3,I))) X(6,I)=SQRT(2.0*EXP(2.0*X(5,I))) X(7,I)=ASIN(X(5,I)) X(8,I)=X0 X(9,I)=Y0 X(10,I)=V(X(5,1)) !需要考虑一下 X(1,I)=V(X(5,I))+Q0 1 CONTINUE RETURN END |
回复此楼» 猜你喜欢
080500求调剂
已经有5人回复
-
269电子信息求调剂，可转专业
已经有5人回复
-
285求调剂
已经有15人回复
-
265求调剂
已经有8人回复
-
274求调剂求调剂
已经有5人回复
-
284求调剂
已经有4人回复
-
土木水利专硕276分求调剂
已经有4人回复
-
材料考研求调剂总分280
已经有9人回复
-
285求调剂
已经有6人回复
-
材料调剂
已经有9人回复
找到一些相关的精华帖子，希望有用哦~
大侠们帮忙看看这是什么意思呢？不懂
已经有5人回复
- 修改意见 各位大侠帮忙看一下什么意思
已经有4人回复
- 急求：各位大侠，哪位知道文章修改重投后ID号会改变吗？
已经有9人回复
- MFC 中单文档应用程序框架 如何调用其他功能的对话框？
已经有13人回复
- 一段中文的翻译，各位大侠能否帮忙修改下
已经有5人回复
- perl程序求助，请大神们赐教~~~
已经有9人回复
- 求帮忙看一下IV curve
已经有3人回复
- 大家爱帮我看看这段MATLAB程序为啥出错了
已经有5人回复
- 请帮忙编写一个perl程序处理文本
已经有26人回复
- 各位大侠帮我解个小困惑啊
已经有16人回复
- 求大侠帮忙解决一个B类错误！！谢谢
已经有5人回复
- 急求一个matlab编码，各位如果知道一定要帮我一把啊
已经有10人回复
- C语言变量问题，求助大侠！！
已经有17人回复
- MC一个简单的程序，运行不了，哪位大侠帮看看怎么回事。多谢！
已经有7人回复
- 请各位大侠帮忙看看如何修改
已经有6人回复
- 请各位大侠帮忙看看：煤堆自燃边界条件设置问题。谢谢！
已经有6人回复
- 求帮忙啊~ matlab一多项式结果太长，怎么写到txt里？
已经有14人回复
- 求各位前辈帮忙！matlab解复数方程，出现图中的警告，怎来修改呢？谢谢了啊
已经有5人回复
- 帮我看看这样明年有没有必要继续申报？
已经有6人回复
- 【求助】求助哪位大侠帮我分析下电子密度差分图
已经有16人回复
- 【求助】帮忙看一个程序
已经有16人回复
- 【求助】有没有哪位大侠报考的定向博士考上后改为非定向的？
已经有10人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
