Sources Esope | Cast3M

fundt
C FUNDT     SOURCE    BECC      09/12/07    21:15:30     6579
      subroutine fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
     &     acv_b, Rgas_b, q, T_ag,
     &     P_r, T_r, u_r,
     &     T_b, P_b, rho_b, u_b, d, k0sdt,
     &     T_a, P_a, rho_a, u_a,
     &     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 (and rs);
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 and (rs)
C
C     acv_b            = to compute cv and ether for rss;
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_b           = gas constant for rss
C
C     q                = released chemical energy
C
C     T_ag             = guess temperature for a (= rs, state ahead the
C                        reactive shock and behind the non-reactive one
C                        in the ZND model)
C
C     T_b              = temperature for rss
C
C
C     OUTPUT
C
C     p_b, u_b, rho_b  = pressure, temperature, velocity and density for
C                        rss (burnt state behind the SDT)
C
C     d                = detonation speed
C
C     k0sdt            = fundamental speed of a WDF going from a to b
C
C     p_a, T_a, u_a,
C     rho_a            = pressure, temperature, velocity and density for
C                        rs (unburnt state ahead the SDT)
C
c      implicit none
      integer nordpo
      real*8 eps
      parameter (eps = 1.0D-12)
      real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
     &     , acv_b(1:(nordpo+1)), Rgas_b, q, T_ag
     &     , rho_r, P_r, T_r, u_r
     &     , T_b, P_b, rho_b, u_b, d
     &     , T_a, P_a, rho_a, u_a
      real*8 e_r, cv_r, h_r, cv_b, e_b, h_b, b, c, csi, c_b, b2, m, ksdt
     &     , cv_ag, e_ag, gam_ag, k0sdt, qref
      logical logan, logdeb
      parameter (logdeb = .false.)
C
      call prith1(nordpo, acv_r, Tmaxcv, T_r, e_r, cv_r)
      h_r = e_r + (Rgas_r * T_r)
      call prith1(nordpo, acv_b, Tmaxcv, T_b, e_b, cv_b)
      h_b = e_b + (Rgas_b * T_b)
      rho_r = P_r / (Rgas_r * T_r)
C
      qref = Rgas_r * T_r
      if ((q .lt. (eps * qref)) .and.
     &     (abs(T_b - T_r) .lt. (eps * T_r)) .and.
     &     (abs(Rgas_r - Rgas_b) .lt. (eps * Rgas_r))) then
C
C******* Constant state
C
         rho_b = rho_r
         T_b = T_r
         P_b = P_r
         u_b = u_r
         k0sdt = (Rgas_r + cv_r) / cv_r
         k0sdt = k0sdt * Rgas_r * T_r
         k0sdt = sqrt(k0sdt)
         d = k0sdt + u_r
C
         T_a = T_r
         P_a = P_R
         rho_a = rho_r
         u_a = u_r
      else
         b = h_b - (h_r + q)
         b = b + (0.5D0 * ((Rgas_r * T_r)  - (Rgas_b * T_b)))
         c = Rgas_r * T_r
         c_b = Rgas_b * T_b
         b2 = b * b
         csi = b2 + (c * c_b)
         if (csi .lt. 0.0D0) then
            write(*,*) 'fundt.f'
            write(*,*) 'delta = ', csi
            write(*,*) 'We take delta = 0'
            csi = 0.0D0
         endif
         csi = sqrt(csi) + b
         if (csi .lt. 0.0D0) then
            write(*,*) 'fundt.f'
            write(*,*) 'negative density'
            logan = .true.
            goto 9999
         endif
         csi = csi / c_b
         rho_b = rho_r * csi
         P_b = Rgas_b * rho_b * T_b
C
         m = (P_b - P_r)
         m = m / ((1.0D0 / rho_r) - (1.0D0 / rho_b))
         if (m .lt. 0.0D0) then
            write(*,*) 'm &lt; 0', m
            m = 0.0D0
         endif
         m = sqrt( m )
         ksdt = m / rho_r
         d = u_r + ksdt
         u_b = d - (m / rho_b)
C
C        State a (ahead the reactive shock and behind the non-reactive
C         shock  in the ZND model)
C
         call prith1(nordpo, acv_r, Tmaxcv, T_ag, e_ag, cv_ag)
         c = 2.0D0 * (h_r - (e_ag - (cv_ag * T_ag)))
         c = c / (ksdt * ksdt)
         c = c + 1.0D0
         gam_ag = (cv_ag + Rgas_r) / cv_ag
         c = ((gam_ag - 1.0D0) / (gam_ag + 1.0D0)) * c
         b = Rgas_r * T_r / (ksdt * ksdt)
         b = b + 1.0D0
         b = (b * gam_ag) / (gam_ag + 1.0D0)
C
         csi = b * b - c
         if (csi .lt. 0.0D0) then
            if (logdeb) then
               write(*,*) 'fundt.f'
               write(*,*) 'delta &lt; 0 . ', csi
               write(*,*) 'Anomaly detected'
            endif
            csi = 0.0D0
C            logan = .true.
C            goto 9999
         endif
         csi = sqrt (csi)
         csi = b - csi
         if (csi .le. 0.0D0) then
            write(*,*) 'fundt.f'
            write(*,*) 'Anomaly detected'
            write(*,*) 'csi = ', csi
            logan = .true.
            goto 9999
         endif
         rho_a = rho_r / csi
         P_a = P_r + ((m * m) * ((1.0D0 / rho_r) - (1.0D0 / rho_a)))
         T_a = P_a / (Rgas_r * rho_a)
         u_a = d - m / rho_a
         k0sdt = d - u_a
C         write(*,*) 'ksdt, k0sdt ', ksdt, k0sdt
      endif
C
 9999 continue
      return
      end