24小时热门版块排行榜

返回列表

20033657

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 640.6
红花: 1
帖子: 70
在线: 111.1小时
虫号: 1000797
注册: 2010-04-19
专业: 结构工程

[交流] 【求助】求大侠看看这个程序怎么修改

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

回复此楼