Télécharger funcj.eso

Retour à la liste

funcj
C FUNCJ     SOURCE    BECC      09/12/07    21:15:28     6579
      subroutine funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
     &     k0, q,
     &     T_ag, T_bg,
     &     p_r, T_r, u_r,
     &     logdef,
     &     r_acj, T_acj, p_acj, u_acj, d_acj,
     &     r_bcj, T_bcj, p_bcj, u_bcj, d_bcj,
     &     k0cjdt, logan)
C
C     INPUT
C
C     nordpo           = order of polynomial for cp and cv (see also
C                        Tmaxcv)
C
C     Tmaxcv           = maximum temperature for cv polynomial expansion
C                        cv(T) = cv(Tmaxcv) if T > Tmaxcv
C
C     acv_r            = to compute cv and ether for r;
C                        vector such that
C                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
C
C     Rgas_r           = gas constant for r
C
C     acv_b            = to compute cv and ether for b (rss) ;
C
C     Rgas_b           = gas constant for (rss)
C
C     k0               = fundamental flame speed
C
C     q                = released heat rs->rss
C
C     T_ag, T_bg       = guess temperature for a (rs) and b (rss)
C                        used to evaluate gamma and linearisation of
C                        enthalpies
C
C     p_r, T_r, u_r    = pressure, temperature and velocities for r
C
C
C
C     OUTPUT
C
C     logdef           = logical. If .true., CJDF is admissible.
C                        Otherwise, it is not and the output is CJDT
C                        state
C
C     d_bcj              = reactive shock speed
C
C     d_acj              = non-reactive shock speed in CJDF case
C                          reactive shock speed in CJDT case
C
C     r_acj, T_acj,
C     p_acj, u_acj     = unburnt state behind the CJDF or CJDT
C
C     r_bcj, T_bcj,
C     p_bcj, u_bcj     = unburnt state behind the CJDF or CJDT
C
C     logan            = if .true., an anomaly has been detected
C
c      implicit none
      integer nordpo
      real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
     &     , acv_b(1:(nordpo+1)), Rgas_b, k0, q
     &     , T_ag, T_bg, p_r, T_r, u_r
     &     , r_acj, T_acj, p_acj, u_acj, d_acj
     &     , r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
      real*8 zero
      parameter(zero = 0.0D0)
      real*8 et_a, cv_a, et_b, cv_b, gam_a, gam_b
     &     , qti, qtiti, a, b, c
     &     , ht_acj, m
     &     , cv_r, et_r, k0cjdt, kcjdt
     &     , r_r
      logical logan, logdef, logdeb
      parameter(logdeb = .false.)
      if (logdeb) then
         write(*,*) 'Subroutine funcj.f'
      endif
C
      logan = .false.
C
C     Properties (gam_a, gam_b)
C
      call prith1(nordpo, acv_r, Tmaxcv, T_ag, et_a, cv_a)
      gam_a = (cv_a + Rgas_r) / cv_a
      call prith1(nordpo, acv_b, Tmaxcv, T_bg, et_b, cv_b)
      gam_b = (cv_b + Rgas_b) / cv_b
C
C**** Correction of q as it is zero
C
      q = max(q, (zero * Rgas_r * T_r))
      qti = q + ((cv_b * T_bg) - et_b)
C     qtiti used in case of CJDF
      qtiti = qti - ((cv_a * T_ag) - et_a)
      qtiti = max(qtiti, 0.0D0)
      if (logdeb )then
         write(*,*) 'gam_a, gam_b', gam_a, gam_b
         write(*,*) 'q, qti, qtiti'
         write(*,*) q, qti, qtiti
      endif
C
C**** Computation of the CJDT
C     It should be done here, in order to know K0,CJDT
C
      call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)
      a = 1.0D0 / gam_b
      a = a * a
      c = Rgas_r * T_r
      c = c * c
      b = 1.0D0 - a
      b = b * ((et_r + (Rgas_r * T_r)) + qti)
      b = b - (Rgas_r * T_r)
      kcjdt = (b*b) - (a*c)
      if ((kcjdt .lt. 0.0D0))then
         kcjdt = 0.0D0
         if (logdeb) then
            write(*,*) 'subroutine funcj'
            write(*,*) 'Computation of kcjdt'
            write(*,*) 'Delta < 0', kcjdt
            write(*,*) 'Delta = ', kcjdt
            write(*,*) 'b*b = ', b*b
            write(*,*) 'sqrt(|Delta|)/b = ', (sqrt(abs(kcjdt))/b)
         endif
C         logan = .true.
C         goto 9999
      endif
      kcjdt = b + sqrt(kcjdt)
      kcjdt = kcjdt / a
      if (kcjdt .lt. 0.0D0)then
         write(*,*) 'subroutine funcj'
         write(*,*) 'Computation of kcjdt'
         write(*,*) 'kcjdt*kcjdt < 0'
         logan = .true.
         goto 9999
      endif
      kcjdt = sqrt(kcjdt)
      d_bcj = kcjdt + u_r
C
C**** Computation of T_bcj (case of cjdt detonation)
C
      r_r = P_r / (Rgas_r * T_r)
      r_bcj = (gam_b - 1.0D0) / (gam_b + 1.0D0)
      r_bcj = r_bcj *
     &     ((2.0D0 * ((et_r + (Rgas_r * T_r)) + qti) / (kcjdt*kcjdt))
     &     + 1.0D0)
      if (r_bcj .lt. 0.0D0) then
         write(*,*) 'subroutine funcj'
         write(*,*) 'Computation of r_bcj'
         write(*,*) 'r_bcj < 0'
         logan = .true.
         goto 9999
      endif
      r_bcj = sqrt(r_bcj)
      r_bcj = r_r / r_bcj
      m = r_r * kcjdt
      p_bcj = p_r -
     &     ((m * m) * ((1.0D0 / r_bcj) - (1.0D0 / r_r)))
      if (p_bcj .lt. 0.0D0) then
         write(*,*) 'subroutine funcj'
         write(*,*) 'Computation of p_bcj'
         write(*,*) 'p_bcj < 0'
         logan = .true.
         goto 9999
      endif
      T_bcj = p_bcj / (Rgas_b * r_bcj)
      u_bcj = d_bcj - (m / r_bcj)
C
C******* Computation of VN state
C
      b = gam_a / (gam_a + 1.0D0)
      b = b * (1.0D0 + (P_r / (r_r * kcjdt * kcjdt)))
      c = 2.0D0 * ((et_r + (Rgas_r * T_r)) - (et_a - (cv_a * T_ag)))
      c = c / (kcjdt * kcjdt)
      c = c + 1.0D0
      c = c * (gam_a - 1.0D0) / (gam_a + 1.0D0)
      r_acj = b*b - c
      if (r_acj .lt. 0.0D0) then
         r_acj = 0.0D0
         if (logdeb) then
            write(*,*) 'subroutine funcj'
            write(*,*) 'Computation of r_acj'
            write(*,*) 'delta < 0', r_acj
         endif
C         logan = .true.
C         goto 9999
      endif
      r_acj = b - sqrt( r_acj )
      r_acj = r_R / r_acj
      m = kcjdt *  r_R
      P_acj = (1.0D0 / r_acj) - (1.0D0 / r_r)
      P_acj = P_R - (m**2 * P_acj)
      if (P_acj .lt. 0.0D0) then
         write(*,*) 'subroutine funcj'
         write(*,*) 'Computation of P_acj'
         write(*,*) 'P_acj < 0'
         logan = .true.
         goto 9999
      endif
      T_acj = P_acj / (Rgas_R * r_acj)
      u_acj = D_BCJ - (m / r_acj)
      k0cjdt = d_bcj - u_acj
CC
CC******* Test. We write the values thus determined
CC
C      write(*,*) 'T_acj,  p_acj,  r_acj,  u_acj '
C      write(*,*) T_acj,  p_acj,  r_acj,  u_acj
C      write(*,*) 'T_bcj, p_bcj, r_bcj, u_bcj, d_bcj'
C      write(*,*) T_bcj, p_bcj, r_bcj, u_bcj, d_bcj
C      stop
C
      if (k0 .ge. k0cjdt) then
         logdef = .false.
         d_acj = d_bcj
      else
C
C        CJDF
C
         logdef = .true.
         a = 1.0D0 / gam_b
         a = a * a
         b = 1.0D0 - a
         b = b * (gam_a / (gam_a - 1.0D0))
         b = b - 1.0D0
         b = b * (k0 * k0)
         c = 2.0D0 * (k0 * k0)
         c = c * (1.0D0 - a)
         c = c * qtiti
         c = c - ((k0**4) * a)
         T_acj = b * b + c
         if (T_acj .lt. 0.0D0) then
            T_acj = 0.0D0
            if (logdeb) then
               write(*,*) 'subroutine funcj'
               write(*,*) 'computation of T_acj'
               write(*,*) 'Delta < 0', T_acj
            endif
C            logan = .true.
C            goto 9999
         endif
         T_acj = sqrt(t_acj)
         T_acj = T_acj + b
         T_acj = T_acj / Rgas_r
C
C        CJDF state in case of right-travelling shock
C
         call funsho(nordpo, Tmaxcv, acv_r, Rgas_r,
     &        P_r, T_r, u_r, .true.,
     &        T_acj, p_acj, r_acj, u_acj, d_acj)
C
         r_bcj = (gam_b - 1.0D0) / (gam_b + 1.0D0)
         ht_acj = gam_a / (gam_a - 1.0D0) * Rgas_r * T_acj
         r_bcj = r_bcj *
     &        ((2.0D0 * (ht_acj + qtiti) / (k0*k0))
     &        + 1.0D0)
         if (r_bcj .lt. 0.0D0) then
            write(*,*) 'subroutine funcj'
            write(*,*) 'computation of r_bcj'
            write(*,*) 'Delta < 0'
            logan = .true.
            goto 9999
         endif
         r_bcj = sqrt(r_bcj)
         r_bcj = r_acj / r_bcj
         m = r_acj * k0
         p_bcj = p_acj -
     &        ((m * m) * ((1.0D0 / r_bcj) - (1.0D0 / r_acj)))
         if (p_bcj .lt. 0.0D0) then
            write(*,*) 'subroutine funcj'
            write(*,*) 'Computation of P_bcj'
            write(*,*) 'p_bcj < 0'
            logan = .true.
            goto 9999
         endif
         T_bcj = p_bcj / (Rgas_b * r_bcj)
         d_bcj = k0 + u_acj
         u_bcj = d_bcj - (m/r_bcj)
CC
CC****    Test. We write the values thus determined
CC
C         write(*,*) 'T_acj,  p_acj,  r_acj,  u_acj, d_acj '
C         write(*,*) T_acj,  p_acj,  r_acj,  u_acj, d_acj
C         write(*,*) 'T_bcj, p_bcj, r_bcj, u_bcj, d_bcj'
C         write(*,*) T_bcj, p_bcj, r_bcj, u_bcj, d_bcj
C         stop
C
         if (d_bcj .gt. d_acj) then
            if (logdeb) then
               write(*,*) 'subroutine funcj'
               write(*,*) 'Deflagration regime'
               write(*,*) 'd_bcj, d_acj ', d_bcj, d_acj
               write(*,*)
     $              'Reactive shock faster than the non-reactive one'
            endif
C            logan = .true.
C            goto 9999
         endif
      endif
C      write(*,*) 'k0cjdt, d_acj, d_bcj ', k0cjdt, d_acj, d_bcj
C
 9999 continue
      return
      end
C
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales