MC DRAG - 발사체의 항력 계수를 추정하는 컴퓨터 프로그램

1981년

 

제2차 세계 대전 이후로 항공기, 미사일 및 무기 발사체의 공기역학적 특성을 추정하는 더 빠르고 정확한 방법에 대한 필요성이 점점 더 커지고 있다. 지난 10년 동안 이러한 필요성은 이용 가능한 데이터의 체계적인 컴파일, 이론적 유동장 솔루션에 기반한 계산, 그리고 위의 조합을 통해 충족되었다.

최근 몇 년 동안 크고 강력한 컴퓨팅 기계의 확산으로 더 빠르고 균일하며 정확한 공기역학적 추정치를 구현하는 것에 대한 광범위한 관심이 높아지고 있다. 유동장 계산을 기반으로 한 접근 방식은 임의의 발사체 형태에 대해 정확도와 근사의 균일성을 향상시킬 수 있는 장기적인 전망을 제공한다. 그러나 더 발전된 컴퓨터에서도 이 접근 방식은 일반적으로 상당히 길고, 지정된 범위의 마하 수, 레이놀즈 수 및 요 레벨에만 적용될 수 있으며, 실제 매끄럽지 않은 형태의 발사체 형태에는 적용하기 어렵다.

공기역학 데이터는 항상 다항식에 맞출 수 있다. 이 과정은 중간 크기의 컴퓨터에서도 빠르게 진행되며, 종종 매우 좋은 적합도를 제공한다. 그러나 이러한 다항식 적합도를 원래 데이터베이스 이상으로 외삽하는 것은 본질적으로 위험하다. 외삽이 필요한 경우, 데이터는 이론에 기초한 방정식에 맞춰 외삽된 영역 전반에 걸쳐 유효해야 한다.

이 보고서에서는 특정 공기역학적 유사성 규칙을 통해 제로 요 항력 계수와 마하 수 사이의 관계를 도출한다. 이 관계는 (a) 특정 형상 및 크기 매개변수와 (b) 최소 제곱에 의해 값이 결정된 추가 매개변수를 포함한다. 이러한 최소 제곱 값은 마하 수 범위 0.5에서 5, 발사체 직경 범위 4에서 400mm에 걸쳐 유효하다. 따라서 이러한 범위 내에서 항력 계수는 주어진 크기 및 형상 매개변수 집합에 대해 추가적인 피팅 과정 없이 직접 계산할 수 있다. 프로그램 MC DRAG가 이 계산을 수행한다. 이 프로그램은 소형 무기 탄환, 재진입체 모델, 포탄 등 세 가지 예시에 적용될 예정이다.

 

C     ******************************************************************
C     MC DRAG
C
C     ESTIMATE OF ZERO-YAW DRAG COEFFICIENT FOR A BODY OF REVOLUITION.
C     INPUTS
C     REFERENCE DIAMETER(MM)
C     TOTAL LENGTH(CAL)
C     NOSE LENGTH(CAL)
C     RATIO OF TANGENT RADIUS TO ACTUAL NOSE RADIUS(HEADSHAPE PARAMFTER)
C     BOATTAIL LENGTH((CAL)
C     BASE DIAMETER(CAL)
C     MEPLAT DIAMETER(CAL)
C     BAND DLAMETER(CAL)
C     CENTER OF GRAVITY(CAL.FROM NOSE)
C     BOUNDARY LAYER CDDE(L/T OR T/T)
C     PROJFCTILE IDENTIFICATION.
C     THE STANDARD DEVIATION OF THE DRAG ESTIMATE IS 10 PECENT AT
C     SUBSONIC AND TRANSONIC SPEEDS, AND 4 PERCENT AT SUPER SONIC SPEEDS.
C     ******************************************************************
      CHARACTER(LEN=100) :: CODE
      DIMENSION CD(24),CDH(24),CDSF(24),CDBND(24),CDBT(24),CDB(24)
      DIMENSION PBP1(24)
      REAL M(24),LT,LN,RTR,LBT,DB,DM,DBND,XCGN,M2
      REAL TA
      INTEGER BLC
      DATA(M(I),I=1,24)/.5,.6,.7,.8,.85,.9,.925,.95,.975,1.,1.1,1.2,1.3,
     11.4,1.5,1.6,1.7,1.8,2.,2.2,2.5,3.,3.5,4./
C    1 CONTINUE
C     READ(5,501) DREF,LT,LN,RTR,LBT,DB,DM,DRB,XCGN,BLC,CODEA,CODEB
      OPEN (UNIT=6,FILE='RESULTS.OUT')
C     ******************************************************************
C     INPUT
C     ******************************************************************
      DREF=155.0              ! PROJFCTILE REFERENCE DIAMETER (MM)
      LT=5.65                 ! PROJECTILE TOTAL LENGTH(CAL)
      LN=3.01                 ! NOSE LENGTH(CAL)
      RTR=0.50                ! HEADSHAPE PARAMFTER
      LBT=0.58                ! BOATTAIL LENGTH((CAL)
      DB=0.848                ! PROJFCTILE BASE DIAMETER(CAL)
      DM=0.090                ! PROJFCTILE MEPLAT DIAMETER(CAL)
      DBND=1.020              ! ROTATING BAND DLAMETER(CAL)
      XCGN=3.53               ! CENTER OF GRAVITY(CAL. FROM NOSE)
      BLC=2                   ! BOUNDARY LAYER CDDE(1: L/T OR 2: T/T)
      CODE='155MM M549'       ! PROJFCTILE IDENTIFICATION
C     ******************************************************************
      WRITE(6,*) CODE
      WRITE(6,"('PROJFCTILE REFERENCE DIAMETER (MM)   = ',F6.1)") DREF
      WRITE(6,"('PROJECTILE TOTAL LENGTH(CAL)         = ',F6.3)") LT
      WRITE(6,"('NOSE LENGTH(CAL)                     = ',F6.3)") LN
      WRITE(6,"('HEADSHAPE PARAMFTER                  = ',F6.2)") RTR
      WRITE(6,"('BOATTAIL LENGTH((CAL)                = ',F6.3)") LBT
      WRITE(6,"('BASE LENGTH((CAL)                    = ',F6.3)") DB
      WRITE(6,"('MEPLAT DIAMETER(CAL)                 = ',F6.3)") DM
      WRITE(6,"('ROTATING BAND DLAMETER(CAL)          = ',F6.3)") DBND
      WRITE(6,"('CENTER OF GRAVITY(CAL. FROM NOSE)    = ',F6.2)") XCGN
      WRITE(6,"('BOUNDARY LAYER CDDE(1: L/T, 2: T/T)  = ',I6)") BLC
      WRITE(6,"('*********************************************')")
      WRITE(6,"('     M       CD0       CDH      CDSF     CDBND
     1CDBT    CDB     PB/P1')")
C     ******************************************************************
      DO 300 I=1,24
      TA=(1.0-DM)/LN  ! THICKNESS RATIO
      M2=M(I)**2.0
      RE=23296.3*M(I)*LT*DREF
      RET=0.4343*(ALOG(RE))
      CFT=(0.455/(RET**2.58))*((1.0+0.21*M2)**(-0.32))    ! PRANDTL'S NUMBER
C   WETTED SURFACE AREA OF THE PROJECTILE NOSE
      DUM=1.0+((0.333+(0.02/(LN**2)))*RTR)
      SWN=1.5708*LN*(1.0+1.0/(8.0*(LN**2)))*DUM
      SWCYL=3.1416*(LT-LN)
      SW=SWN+SWCYL
C   ------------------------------------------------------------------
C     LAMINAR BOUNDARY LAYER
      IF (BLC.EQ.1) THEN
      CFL=(1.328/(SQRT(RE)))*((1.0+0.12*M2)**(-0.12))
      END IF
C   TURBULENT BOUNDARY LAYER
      IF (BLC.EQ.2) THEN
      CFL=CFT
      END IF
      CDSFL=1.2732*SW*CFL
      CDSFT=1.2732*SW*CFT
      CDSF(I)=(CDSFL*SWN+CDSFT*SWCYL)/SW
C     ******************************************************************
C     PRESSURE DRAG COEFFICIENT FOR A ROTATING BAND
      CHI=(M2-1.0)/(2.4*M2)
      IF(M(I).LE.1.0) PTP=(1.0+0.2*M2)**3.5
      IF(M(I).GT.1.0) PTP=((1.2*M2)**3.5)*((6.0/(7.0*M2-1.0))**2.5)
      CMEP=(1.122*(PTP-1.0)*(DM*DM))/M2
      IF(M(I).LE.0.91) CDHM=0.0
      IF(M(I).GE.1.41) CDHM=0.85*CMEP
      IF(M(I).GT.0.91.AND.M(I).LT.1.41) CDHM=(0.254+2.88*CHI)*CMEP
      IF(M(I).LT.1.00) PB2=1.0/(1.0+0.1875*M2+0.0531*M2*M2)
      IF(M(I).GE.1.00) PB2=1.0/(1.0+0.2477*M2+0.0345*M2*M2)
      PB4=(1.0+0.09*M2*(1.0-EXP(LN-LT)))*(1.0+0.25*M2*(1.0-DB))
      PBP1(I)=PB2*PB4
      CDB(I)=(1.4286*(1.0-PBP1(I))*(DB*DB))/M2
      IF(M(I).LT.0.95) CDBND(I)=(M(I)**12.5)*(DBND-1.0)
      IF(M(I).GE.0.95) CDBND(I)=(0.21+0.28/M2)*(DBND-1.0)
C     ******************************************************************
C     TOTAL DRAG COEFFICIENT AT ZERO ANGLE OF ATTACK
      IF(M(I)-1.0) 100,100,200
C     SUBSONIC-TRANSONIC SPEEDS.
  100 CONTINUE
C     LBT: BOATTAIL LENGTH
      IF(LBT) 102,101,102
  101 CONTINUE
C     WHEN THERE IS NO BOATTAIL
      CDBT(I)=0.0
      CD(I)=CDH(I)+CDSF(I)+CDBND(I)+CDBT(I)+CDB(I)
      WRITE(6,1509) M(I),CD(I),CDH(I),CDSF(I),CDBND(I),CDBT(I),CDB(I),
     1 PBP1(I)
      GOTO 105
  102 CONTINUE
C   WHEN THERE IS A BOATTAIL
      IF(M(I).LE.0.85) GO TO 101
      TB=(1.0-DB)/(2.0*LBT)
      TB23=2.0*TB*TB+(TB**3)
      EBT=EXP((-2.0)*LBT)
      BBT=1.0-EBT+2.0*TB*((EBT*(LBT+0.5))-0.5)
      CDBT(I)=2.0*TB23*BBT*(1.0/(0.564+1250.0*CHI*CHI))
  105 CONTINUE
      XMC=(1.0+0.552*(TA**0.8))**(-0.5)
      IF(M(I).LE.XMC) CDHT=0.0
      IF(M(I).GT.XMC) CDHT=0.368*(TA**1.8)+1.6*TA*CHI
      CDH(I)=CDHT+CDHM
      CD(I)=CDH(I)+CDSF(I)+CDBND(I)+CDBT(I)+CDB(I)
      WRITE(6,1509) M(I),CD(I),CDH(I),CDSF(I),CDBND(I),CDBT(I),CDB(I),
     1 PBP1(I)
      GOTO 300
C   ------------------------------------------------------------------
C     SUPERSONIC SPEEDS.
C   ------------------------------------------------------------------
  200 CONTINUE
      RE2=M2-1.0
      RE=SQRT(RE2)
      ZE=RE
      SSMC=1.0+0.368*(TA**1.85)
      IF(M(I).LT.SSMC) ZE=SQRT(SSMC*SSMC-1.0)
      C1=0.7156-0.5313*RTR+0.595*(RTR**2)
      C2=0.0796+0.0779*RTR
      C3=1.5870+0.0490*RTR
      C4=0.1122+0.1658*RTR
      RZ2=1.0/(ZE*ZE)
      CDHT=(C1-C2*(TA**2))*RZ2*((TA*ZE)**(C3+C4*TA))
      CDH(I)=CDHT+CDHM
C     ******************************************************************
C     PRESSURE DRAG COEFFICIENT DUE TO BOATTAIL (OR FLARE)
      IF(LBT) 202,201,202
  201 CONTINUE
C     WHEN THERE IS NO BOATTAIL
      CDBT(I)=0.0
      CD(I)=CDH(I)+CDSF(I)+CDBND(I)+CDBT(I)+CDB(I)
      WRITE(6,1509) M(I),CD(I),CDH(I),CDSF(I),CDBND(I),CDBT(I),CDB(I),
     1 PBP1(I)
      GOTO 300
  202 CONTINUE
C   WHEN THERE IS A BOATTAIL
      TB=(1.0-DB)/(2.0*LBT)
C     ******************************************************************
      IF(M(I)-1.1) 205,205,207
  205 CONTINUE
C     SUBSONIC-TRANSONIC SPEEDS.
      TB23=2.0*TB*TB+(TB**3)
      EBT=EXP((-2.0)*LBT)
      BBT=1.0-EBT+2.0*TB*((EBT*(LBT+0.5))-0.5)
      CDBT(I)=2.0*TB23*BBT*(1.774-9.3*CHI)
      CD(I)=CDH(I)+CDSF(I)+CDBND(I)+CDBT(I)+CDB(I)
      WRITE(6,1509) M(I),CD(I),CDH(I),CDSF(I),CDBND(I),CDBT(I),CDB(I),
     1 PBP1(I)
      GOTO 300
  207 CONTINUE
C     SUPERSONIC SPEEDS.
      BR=0.85/RE       ! K
      AA2=((5.0*TA)/(6.0*RE))*(0.5*TA)**2-(0.7435/M2)*((TA*M(I))**1.6)
      AA1=(1.0-(3*RTR)/(5*M(I)))*AA2
      EXL=EXP(((-1.1952)/M(I))*(LT-LN-LBT))
      XXM=((2.4*M2*M2-4.0*RE2)*(TB*TB))/(2.0*RE2*RE2)
      AA=AA1*EXL+((2.0*TB)/RE)-XXM
      RB=1.0/BR        ! 1/K
      EXRT=EXP(-BR*LBT)
      AAB=(1.0-EXRT)+(2.0*TB*(EXRT*(LBT+RB)-RB))
      CDBT(I)=4.0*AA*TB*RB*AAB
      CD(I)=CDH(I)+CDSF(I)+CDBND(I)+CDBT(I)+CDB(I)
      WRITE(6,1509) M(I),CD(I),CDH(I),CDSF(I),CDBND(I),CDBT(I),CDB(I),
     1 PBP1(I)
  300 CONTINUE
C     ******************************************************************
 1509 FORMAT(F6.3,3X,7(F7.3,3X))
      IF(LN.LT.1.0.OR.DM.GT.0.5) GOTO 698
      GOTO 1510
  698 WRITE(6,1512)
 1512 FORMAT("WARNING... NOSE TOO SHORT OR TOO BLUNT. CHECK CDH.")
 1510 IF(LBT.GT.1.5.OR.DB.LT.0.65) GOTO 699
      GOTO 9999
  699 WRITE(6,1513)
 1513 FORMAT("WARNING... BOATTAIL TOO LONG OR TOO STEEP. CHECKC CDBT.")
 9999 CLOSE(6)
      END