Numérotation des lignes :

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)CC     INPUTCC     nordpo           = order of polynomial for cp and cv (see alsoC                        Tmaxcv)CC     Tmaxcv           = maximum temperature for cv polynomial expansionC                        cv(T) = cv(Tmaxcv) if T > TmaxcvCC     acv_r            = to compute cv and ether for r (and rs);C                        vector such thatC                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)CC     Rgas_r           = gas constant for r and (rs)CC     acv_b            = to compute cv and ether for rss;C                        vector such thatC                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)CC     Rgas_b           = gas constant for rssCC     q                = released chemical energyCC     T_ag             = guess temperature for a (= rs, state ahead theC                        reactive shock and behind the non-reactive oneC                        in the ZND model)CC     T_b              = temperature for rssCCC     OUTPUTCC     p_b, u_b, rho_b  = pressure, temperature, velocity and density forC                        rss (burnt state behind the SDT)CC     d                = detonation speedCC     k0sdt            = fundamental speed of a WDF going from a to bCC     p_a, T_a, u_a,C     rho_a            = pressure, temperature, velocity and density forC                        rs (unburnt state ahead the SDT)Cc      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))) thenCC******* Constant stateC         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_rC         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_bC         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)CC        State a (ahead the reactive shock and behind the non-reactiveC         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.0D0C            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_aC         write(*,*) 'ksdt, k0sdt ', ksdt, k0sdt      endifC 9999 continue      return      end

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