Numérotation des lignes :

C FUNSDT    SOURCE    BECC      09/12/07    21:15:34     6579      subroutine funsdt(nordpo, Tmaxcv, acv_r, Rgas_r,     &     acv_b, Rgas_b, q, T_ag, T_bg,     &     P_r, T_r, u_r,     &     T_a, P_a, r_a, u_a, d, k0sdt,     &     T_b, P_b, r_b, u_b,     &     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_bg             = guess temperature for b (= rss, state behind theC                        reactive shock in the ZND model)CC     T_a              = temperature for rsCCC     OUTPUTCC     p_b, T_b, u_b,C     r_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, u_a,C     r_a            = pressure, velocity and density forC                        rs (unburnt state ahead the SDT)Cc      implicit none      integer nordpo      real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r     &     , acv_b(1:(nordpo+1)), Rgas_b, q, T_ag, T_bg     &     , r_r, P_r, T_r, u_r     &     , T_a, P_a, r_a, u_a, d     &     , T_b, P_b, r_b, u_b      real*8  k0sdt, e_a, cv_a, h_a, gam_a     &     , gam_b, e_b, cv_b     &     , c, b, csi, m      logical loganCC**** Unburnt stateC      r_r = P_r / (Rgas_r * T_r)      call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,     &     P_r, T_r, u_r, .true.,     &     T_a, P_a, r_a, u_a, d)CC**** Now we search for the burnt state which travels with theC     same speedC      k0sdt = d - u_a      call prith1(nordpo, acv_r, Tmaxcv, T_ag, e_a, cv_a)      gam_a = (cv_a + Rgas_r) / cv_a      call prith1(nordpo, acv_b, Tmaxcv, T_bg, e_b, cv_b)      gam_b = (cv_b + Rgas_b) / cv_bC      h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_r * T_a      c = h_a + q  - (e_b - (cv_b * T_bg)) +     &     (e_a - (cv_a * T_ag))      c = 2.0D0 * c      c = c / (k0sdt * k0sdt)      c = c + 1.0D0      c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)      if (c .lt. 0.0D0) then         write(*,*) 'funsdt.f'         write(*,*) 'c &lt; 0'         write(*,*) 'Anomaly detected'         logan = .true.         goto 9999      endifC      b = Rgas_r * T_a      b = b / (k0sdt * k0sdt)      b = b + 1.0D0      b = b * gam_b / (gam_b + 1.0D0)      csi = b * b - c      if (csi .lt. 0.0D0) thenC        Non itersection between the rayleigh lineC        and the Hugoniot adiabatic         if (csi .lt. 1.0D-3) thenC        csi can be considered as zero            csi = 0.0D0         else            write(*,*) 'csi = ', csi            write(*,*) 'csi &lt; 0'            write(*,*) 'Anomaly detected'            logan = .true.            goto 9999         endif      else         csi = sqrt(csi)      endifC     Computation of the SDT solutionC     For the WDT solution, csi = b + csi      csi = b - csi      r_b = r_a / csi      m = k0sdt * r_a      p_b = p_a + ((m*m) * ((1.0D0 / r_a) - (1.0D0 / r_b)))      t_b = P_b / (Rgas_b * r_b)      u_b = d - (m / r_b)C 9999 continue      return      end

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