Télécharger fundt.eso

Retour à la liste

Numérotation des lignes :

  1. C FUNDT SOURCE BECC 09/12/07 21:15:30 6579
  2. subroutine fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
  3. & acv_b, Rgas_b, q, T_ag,
  4. & P_r, T_r, u_r,
  5. & T_b, P_b, rho_b, u_b, d, k0sdt,
  6. & T_a, P_a, rho_a, u_a,
  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_b = temperature for rss
  38. C
  39. C
  40. C OUTPUT
  41. C
  42. C p_b, u_b, rho_b = pressure, temperature, velocity and density for
  43. C rss (burnt state behind the SDT)
  44. C
  45. C d = detonation speed
  46. C
  47. C k0sdt = fundamental speed of a WDF going from a to b
  48. C
  49. C p_a, T_a, u_a,
  50. C rho_a = pressure, temperature, velocity and density for
  51. C rs (unburnt state ahead the SDT)
  52. C
  53. c implicit none
  54. integer nordpo
  55. real*8 eps
  56. parameter (eps = 1.0D-12)
  57. real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
  58. & , acv_b(1:(nordpo+1)), Rgas_b, q, T_ag
  59. & , rho_r, P_r, T_r, u_r
  60. & , T_b, P_b, rho_b, u_b, d
  61. & , T_a, P_a, rho_a, u_a
  62. real*8 e_r, cv_r, h_r, cv_b, e_b, h_b, b, c, csi, c_b, b2, m, ksdt
  63. & , cv_ag, e_ag, gam_ag, k0sdt, qref
  64. logical logan, logdeb
  65. parameter (logdeb = .false.)
  66. C
  67. call prith1(nordpo, acv_r, Tmaxcv, T_r, e_r, cv_r)
  68. h_r = e_r + (Rgas_r * T_r)
  69. call prith1(nordpo, acv_b, Tmaxcv, T_b, e_b, cv_b)
  70. h_b = e_b + (Rgas_b * T_b)
  71. rho_r = P_r / (Rgas_r * T_r)
  72. C
  73. qref = Rgas_r * T_r
  74. if ((q .lt. (eps * qref)) .and.
  75. & (abs(T_b - T_r) .lt. (eps * T_r)) .and.
  76. & (abs(Rgas_r - Rgas_b) .lt. (eps * Rgas_r))) then
  77. C
  78. C******* Constant state
  79. C
  80. rho_b = rho_r
  81. T_b = T_r
  82. P_b = P_r
  83. u_b = u_r
  84. k0sdt = (Rgas_r + cv_r) / cv_r
  85. k0sdt = k0sdt * Rgas_r * T_r
  86. k0sdt = sqrt(k0sdt)
  87. d = k0sdt + u_r
  88. C
  89. T_a = T_r
  90. P_a = P_R
  91. rho_a = rho_r
  92. u_a = u_r
  93. else
  94. b = h_b - (h_r + q)
  95. b = b + (0.5D0 * ((Rgas_r * T_r) - (Rgas_b * T_b)))
  96. c = Rgas_r * T_r
  97. c_b = Rgas_b * T_b
  98. b2 = b * b
  99. csi = b2 + (c * c_b)
  100. if (csi .lt. 0.0D0) then
  101. write(*,*) 'fundt.f'
  102. write(*,*) 'delta = ', csi
  103. write(*,*) 'We take delta = 0'
  104. csi = 0.0D0
  105. endif
  106. csi = sqrt(csi) + b
  107. if (csi .lt. 0.0D0) then
  108. write(*,*) 'fundt.f'
  109. write(*,*) 'negative density'
  110. logan = .true.
  111. goto 9999
  112. endif
  113. csi = csi / c_b
  114. rho_b = rho_r * csi
  115. P_b = Rgas_b * rho_b * T_b
  116. C
  117. m = (P_b - P_r)
  118. m = m / ((1.0D0 / rho_r) - (1.0D0 / rho_b))
  119. if (m .lt. 0.0D0) then
  120. write(*,*) 'm < 0', m
  121. m = 0.0D0
  122. endif
  123. m = sqrt( m )
  124. ksdt = m / rho_r
  125. d = u_r + ksdt
  126. u_b = d - (m / rho_b)
  127. C
  128. C State a (ahead the reactive shock and behind the non-reactive
  129. C shock in the ZND model)
  130. C
  131. call prith1(nordpo, acv_r, Tmaxcv, T_ag, e_ag, cv_ag)
  132. c = 2.0D0 * (h_r - (e_ag - (cv_ag * T_ag)))
  133. c = c / (ksdt * ksdt)
  134. c = c + 1.0D0
  135. gam_ag = (cv_ag + Rgas_r) / cv_ag
  136. c = ((gam_ag - 1.0D0) / (gam_ag + 1.0D0)) * c
  137. b = Rgas_r * T_r / (ksdt * ksdt)
  138. b = b + 1.0D0
  139. b = (b * gam_ag) / (gam_ag + 1.0D0)
  140. C
  141. csi = b * b - c
  142. if (csi .lt. 0.0D0) then
  143. if (logdeb) then
  144. write(*,*) 'fundt.f'
  145. write(*,*) 'delta < 0 . ', csi
  146. write(*,*) 'Anomaly detected'
  147. endif
  148. csi = 0.0D0
  149. C logan = .true.
  150. C goto 9999
  151. endif
  152. csi = sqrt (csi)
  153. csi = b - csi
  154. if (csi .le. 0.0D0) then
  155. write(*,*) 'fundt.f'
  156. write(*,*) 'Anomaly detected'
  157. write(*,*) 'csi = ', csi
  158. logan = .true.
  159. goto 9999
  160. endif
  161. rho_a = rho_r / csi
  162. P_a = P_r + ((m * m) * ((1.0D0 / rho_r) - (1.0D0 / rho_a)))
  163. T_a = P_a / (Rgas_r * rho_a)
  164. u_a = d - m / rho_a
  165. k0sdt = d - u_a
  166. C write(*,*) 'ksdt, k0sdt ', ksdt, k0sdt
  167. endif
  168. C
  169. 9999 continue
  170. return
  171. end
  172.  
  173.  

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