Télécharger la plaquette Cast3M


fr en

Télécharger funsdt.eso

Retour à la liste

Numérotation des lignes : 

funsdt

  1. C FUNSDT    SOURCE    BECC      09/12/07    21:15:34     6579
  2.       subroutine funsdt(nordpo, Tmaxcv, acv_r, Rgas_r,
  3.      &     acv_b, Rgas_b, q, T_ag, T_bg,
  4.      &     P_r, T_r, u_r,
  5.      &     T_a, P_a, r_a, u_a, d, k0sdt,
  6.      &     T_b, P_b, r_b, u_b,
  7.      &     logan)
  8. C
  9. C     INPUT
  10. C
  11. C     nordpo           = order of polynomial for cp and cv (see also
  12. C                        Tmaxcv)
  13. C
  14. C     Tmaxcv           = maximum temperature for cv polynomial expansion
  15. C                        cv(T) = cv(Tmaxcv) if T > Tmaxcv
  16. C
  17. C     acv_r            = to compute cv and ether for r (and rs);
  18. C                        vector such that
  19. C                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
  20. C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
  21. C
  22. C     Rgas_r           = gas constant for r and (rs)
  23. C
  24. C     acv_b            = to compute cv and ether for rss;
  25. C                        vector such that
  26. C                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
  27. C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
  28. C
  29. C     Rgas_b           = gas constant for rss
  30. C
  31. C     q                = released chemical energy
  32. C
  33. C     T_ag             = guess temperature for a (= rs, state ahead the
  34. C                        reactive shock and behind the non-reactive one
  35. C                        in the ZND model)
  36. C
  37. C     T_bg             = guess temperature for b (= rss, state behind the
  38. C                        reactive shock in the ZND model)
  39. C
  40. C     T_a              = temperature for rs
  41. C
  42. C
  43. C     OUTPUT
  44. C
  45. C     p_b, T_b, u_b,
  46. C     r_b            = pressure, temperature, velocity and density for
  47. C                        rss (burnt state behind the SDT)
  48. C
  49. C     d                = detonation speed
  50. C
  51. C     k0sdt            = fundamental speed of a WDF going from a to b
  52. C
  53. C     p_a, u_a,
  54. C     r_a            = pressure, velocity and density for
  55. C                        rs (unburnt state ahead the SDT)
  56. C
  57. c      implicit none
  58.       integer nordpo
  59.       real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
  60.      &     , acv_b(1:(nordpo+1)), Rgas_b, q, T_ag, T_bg
  61.      &     , r_r, P_r, T_r, u_r
  62.      &     , T_a, P_a, r_a, u_a, d
  63.      &     , T_b, P_b, r_b, u_b
  64.       real*8  k0sdt, e_a, cv_a, h_a, gam_a
  65.      &     , gam_b, e_b, cv_b
  66.      &     , c, b, csi, m
  67.       logical logan
  68. C
  69. C**** Unburnt state
  70. C
  71.       r_r = P_r / (Rgas_r * T_r)
  72.       call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  73.      &     P_r, T_r, u_r, .true.,
  74.      &     T_a, P_a, r_a, u_a, d)
  75. C
  76. C**** Now we search for the burnt state which travels with the
  77. C     same speed
  78. C
  79.       k0sdt = d - u_a
  80.       call prith1(nordpo, acv_r, Tmaxcv, T_ag, e_a, cv_a)
  81.       gam_a = (cv_a + Rgas_r) / cv_a
  82.       call prith1(nordpo, acv_b, Tmaxcv, T_bg, e_b, cv_b)
  83.       gam_b = (cv_b + Rgas_b) / cv_b
  84. C
  85.       h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_r * T_a
  86.       c = h_a + q  - (e_b - (cv_b * T_bg)) +
  87.      &     (e_a - (cv_a * T_ag))
  88.       c = 2.0D0 * c
  89.       c = c / (k0sdt * k0sdt)
  90.       c = c + 1.0D0
  91.       c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
  92.       if (c .lt. 0.0D0) then
  93.          write(*,*) 'funsdt.f'
  94.          write(*,*) 'c < 0'
  95.          write(*,*) 'Anomaly detected'
  96.          logan = .true.
  97.          goto 9999
  98.       endif
  99. C
  100.       b = Rgas_r * T_a
  101.       b = b / (k0sdt * k0sdt)
  102.       b = b + 1.0D0
  103.       b = b * gam_b / (gam_b + 1.0D0)
  104.       csi = b * b - c
  105.       if (csi .lt. 0.0D0) then
  106. C        Non itersection between the rayleigh line
  107. C        and the Hugoniot adiabatic
  108.          if (csi .lt. 1.0D-3) then
  109. C        csi can be considered as zero
  110.             csi = 0.0D0
  111.          else
  112.             write(*,*) 'csi = ', csi
  113.             write(*,*) 'csi < 0'
  114.             write(*,*) 'Anomaly detected'
  115.             logan = .true.
  116.             goto 9999
  117.          endif
  118.       else
  119.          csi = sqrt(csi)
  120.       endif
  121. C     Computation of the SDT solution
  122. C     For the WDT solution, csi = b + csi
  123.       csi = b - csi
  124.       r_b = r_a / csi
  125.       m = k0sdt * r_a
  126.       p_b = p_a + ((m*m) * ((1.0D0 / r_a) - (1.0D0 / r_b)))
  127.       t_b = P_b / (Rgas_b * r_b)
  128.       u_b = d - (m / r_b)
  129. C
  130.  9999 continue
  131.       return
  132.       end
  133.  
  134.  

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