Télécharger funsdt.eso

Retour à la liste

Numérotation des lignes :

  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