Télécharger funcj.eso

Retour à la liste

Numérotation des lignes :

  1. C FUNCJ SOURCE BECC 09/12/07 21:15:28 6579
  2. subroutine funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
  3. & k0, q,
  4. & T_ag, T_bg,
  5. & p_r, T_r, u_r,
  6. & logdef,
  7. & r_acj, T_acj, p_acj, u_acj, d_acj,
  8. & r_bcj, T_bcj, p_bcj, u_bcj, d_bcj,
  9. & k0cjdt, logan)
  10. C
  11. C INPUT
  12. C
  13. C nordpo = order of polynomial for cp and cv (see also
  14. C Tmaxcv)
  15. C
  16. C Tmaxcv = maximum temperature for cv polynomial expansion
  17. C cv(T) = cv(Tmaxcv) if T > Tmaxcv
  18. C
  19. C acv_r = to compute cv and ether for r;
  20. C vector such that
  21. C cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
  22. C ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
  23. C
  24. C Rgas_r = gas constant for r
  25. C
  26. C acv_b = to compute cv and ether for b (rss) ;
  27. C
  28. C Rgas_b = gas constant for (rss)
  29. C
  30. C k0 = fundamental flame speed
  31. C
  32. C q = released heat rs->rss
  33. C
  34. C T_ag, T_bg = guess temperature for a (rs) and b (rss)
  35. C used to evaluate gamma and linearisation of
  36. C enthalpies
  37. C
  38. C p_r, T_r, u_r = pressure, temperature and velocities for r
  39. C
  40. C
  41. C
  42. C OUTPUT
  43. C
  44. C logdef = logical. If .true., CJDF is admissible.
  45. C Otherwise, it is not and the output is CJDT
  46. C state
  47. C
  48. C d_bcj = reactive shock speed
  49. C
  50. C d_acj = non-reactive shock speed in CJDF case
  51. C reactive shock speed in CJDT case
  52. C
  53. C r_acj, T_acj,
  54. C p_acj, u_acj = unburnt state behind the CJDF or CJDT
  55. C
  56. C r_bcj, T_bcj,
  57. C p_bcj, u_bcj = unburnt state behind the CJDF or CJDT
  58. C
  59. C logan = if .true., an anomaly has been detected
  60. C
  61. c implicit none
  62. integer nordpo
  63. real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
  64. & , acv_b(1:(nordpo+1)), Rgas_b, k0, q
  65. & , T_ag, T_bg, p_r, T_r, u_r
  66. & , r_acj, T_acj, p_acj, u_acj, d_acj
  67. & , r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
  68. real*8 zero
  69. parameter(zero = 0.0D0)
  70. real*8 et_a, cv_a, et_b, cv_b, gam_a, gam_b
  71. & , qti, qtiti, a, b, c
  72. & , ht_acj, m
  73. & , cv_r, et_r, k0cjdt, kcjdt
  74. & , r_r
  75. logical logan, logdef, logdeb
  76. parameter(logdeb = .false.)
  77. if (logdeb) then
  78. write(*,*) 'Subroutine funcj.f'
  79. endif
  80. C
  81. logan = .false.
  82. C
  83. C Properties (gam_a, gam_b)
  84. C
  85. call prith1(nordpo, acv_r, Tmaxcv, T_ag, et_a, cv_a)
  86. gam_a = (cv_a + Rgas_r) / cv_a
  87. call prith1(nordpo, acv_b, Tmaxcv, T_bg, et_b, cv_b)
  88. gam_b = (cv_b + Rgas_b) / cv_b
  89. C
  90. C**** Correction of q as it is zero
  91. C
  92. q = max(q, (zero * Rgas_r * T_r))
  93. qti = q + ((cv_b * T_bg) - et_b)
  94. C qtiti used in case of CJDF
  95. qtiti = qti - ((cv_a * T_ag) - et_a)
  96. qtiti = max(qtiti, 0.0D0)
  97. if (logdeb )then
  98. write(*,*) 'gam_a, gam_b', gam_a, gam_b
  99. write(*,*) 'q, qti, qtiti'
  100. write(*,*) q, qti, qtiti
  101. endif
  102. C
  103. C**** Computation of the CJDT
  104. C It should be done here, in order to know K0,CJDT
  105. C
  106. call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)
  107. a = 1.0D0 / gam_b
  108. a = a * a
  109. c = Rgas_r * T_r
  110. c = c * c
  111. b = 1.0D0 - a
  112. b = b * ((et_r + (Rgas_r * T_r)) + qti)
  113. b = b - (Rgas_r * T_r)
  114. kcjdt = (b*b) - (a*c)
  115. if ((kcjdt .lt. 0.0D0))then
  116. kcjdt = 0.0D0
  117. if (logdeb) then
  118. write(*,*) 'subroutine funcj'
  119. write(*,*) 'Computation of kcjdt'
  120. write(*,*) 'Delta < 0', kcjdt
  121. write(*,*) 'Delta = ', kcjdt
  122. write(*,*) 'b*b = ', b*b
  123. write(*,*) 'sqrt(|Delta|)/b = ', (sqrt(abs(kcjdt))/b)
  124. endif
  125. C logan = .true.
  126. C goto 9999
  127. endif
  128. kcjdt = b + sqrt(kcjdt)
  129. kcjdt = kcjdt / a
  130. if (kcjdt .lt. 0.0D0)then
  131. write(*,*) 'subroutine funcj'
  132. write(*,*) 'Computation of kcjdt'
  133. write(*,*) 'kcjdt*kcjdt < 0'
  134. logan = .true.
  135. goto 9999
  136. endif
  137. kcjdt = sqrt(kcjdt)
  138. d_bcj = kcjdt + u_r
  139. C
  140. C**** Computation of T_bcj (case of cjdt detonation)
  141. C
  142. r_r = P_r / (Rgas_r * T_r)
  143. r_bcj = (gam_b - 1.0D0) / (gam_b + 1.0D0)
  144. r_bcj = r_bcj *
  145. & ((2.0D0 * ((et_r + (Rgas_r * T_r)) + qti) / (kcjdt*kcjdt))
  146. & + 1.0D0)
  147. if (r_bcj .lt. 0.0D0) then
  148. write(*,*) 'subroutine funcj'
  149. write(*,*) 'Computation of r_bcj'
  150. write(*,*) 'r_bcj < 0'
  151. logan = .true.
  152. goto 9999
  153. endif
  154. r_bcj = sqrt(r_bcj)
  155. r_bcj = r_r / r_bcj
  156. m = r_r * kcjdt
  157. p_bcj = p_r -
  158. & ((m * m) * ((1.0D0 / r_bcj) - (1.0D0 / r_r)))
  159. if (p_bcj .lt. 0.0D0) then
  160. write(*,*) 'subroutine funcj'
  161. write(*,*) 'Computation of p_bcj'
  162. write(*,*) 'p_bcj < 0'
  163. logan = .true.
  164. goto 9999
  165. endif
  166. T_bcj = p_bcj / (Rgas_b * r_bcj)
  167. u_bcj = d_bcj - (m / r_bcj)
  168. C
  169. C******* Computation of VN state
  170. C
  171. b = gam_a / (gam_a + 1.0D0)
  172. b = b * (1.0D0 + (P_r / (r_r * kcjdt * kcjdt)))
  173. c = 2.0D0 * ((et_r + (Rgas_r * T_r)) - (et_a - (cv_a * T_ag)))
  174. c = c / (kcjdt * kcjdt)
  175. c = c + 1.0D0
  176. c = c * (gam_a - 1.0D0) / (gam_a + 1.0D0)
  177. r_acj = b*b - c
  178. if (r_acj .lt. 0.0D0) then
  179. r_acj = 0.0D0
  180. if (logdeb) then
  181. write(*,*) 'subroutine funcj'
  182. write(*,*) 'Computation of r_acj'
  183. write(*,*) 'delta < 0', r_acj
  184. endif
  185. C logan = .true.
  186. C goto 9999
  187. endif
  188. r_acj = b - sqrt( r_acj )
  189. r_acj = r_R / r_acj
  190. m = kcjdt * r_R
  191. P_acj = (1.0D0 / r_acj) - (1.0D0 / r_r)
  192. P_acj = P_R - (m**2 * P_acj)
  193. if (P_acj .lt. 0.0D0) then
  194. write(*,*) 'subroutine funcj'
  195. write(*,*) 'Computation of P_acj'
  196. write(*,*) 'P_acj < 0'
  197. logan = .true.
  198. goto 9999
  199. endif
  200. T_acj = P_acj / (Rgas_R * r_acj)
  201. u_acj = D_BCJ - (m / r_acj)
  202. k0cjdt = d_bcj - u_acj
  203. CC
  204. CC******* Test. We write the values thus determined
  205. CC
  206. C write(*,*) 'T_acj, p_acj, r_acj, u_acj '
  207. C write(*,*) T_acj, p_acj, r_acj, u_acj
  208. C write(*,*) 'T_bcj, p_bcj, r_bcj, u_bcj, d_bcj'
  209. C write(*,*) T_bcj, p_bcj, r_bcj, u_bcj, d_bcj
  210. C stop
  211. C
  212. if (k0 .ge. k0cjdt) then
  213. logdef = .false.
  214. d_acj = d_bcj
  215. else
  216. C
  217. C CJDF
  218. C
  219. logdef = .true.
  220. a = 1.0D0 / gam_b
  221. a = a * a
  222. b = 1.0D0 - a
  223. b = b * (gam_a / (gam_a - 1.0D0))
  224. b = b - 1.0D0
  225. b = b * (k0 * k0)
  226. c = 2.0D0 * (k0 * k0)
  227. c = c * (1.0D0 - a)
  228. c = c * qtiti
  229. c = c - ((k0**4) * a)
  230. T_acj = b * b + c
  231. if (T_acj .lt. 0.0D0) then
  232. T_acj = 0.0D0
  233. if (logdeb) then
  234. write(*,*) 'subroutine funcj'
  235. write(*,*) 'computation of T_acj'
  236. write(*,*) 'Delta < 0', T_acj
  237. endif
  238. C logan = .true.
  239. C goto 9999
  240. endif
  241. T_acj = sqrt(t_acj)
  242. T_acj = T_acj + b
  243. T_acj = T_acj / Rgas_r
  244. C
  245. C CJDF state in case of right-travelling shock
  246. C
  247. call funsho(nordpo, Tmaxcv, acv_r, Rgas_r,
  248. & P_r, T_r, u_r, .true.,
  249. & T_acj, p_acj, r_acj, u_acj, d_acj)
  250. C
  251. r_bcj = (gam_b - 1.0D0) / (gam_b + 1.0D0)
  252. ht_acj = gam_a / (gam_a - 1.0D0) * Rgas_r * T_acj
  253. r_bcj = r_bcj *
  254. & ((2.0D0 * (ht_acj + qtiti) / (k0*k0))
  255. & + 1.0D0)
  256. if (r_bcj .lt. 0.0D0) then
  257. write(*,*) 'subroutine funcj'
  258. write(*,*) 'computation of r_bcj'
  259. write(*,*) 'Delta < 0'
  260. logan = .true.
  261. goto 9999
  262. endif
  263. r_bcj = sqrt(r_bcj)
  264. r_bcj = r_acj / r_bcj
  265. m = r_acj * k0
  266. p_bcj = p_acj -
  267. & ((m * m) * ((1.0D0 / r_bcj) - (1.0D0 / r_acj)))
  268. if (p_bcj .lt. 0.0D0) then
  269. write(*,*) 'subroutine funcj'
  270. write(*,*) 'Computation of P_bcj'
  271. write(*,*) 'p_bcj < 0'
  272. logan = .true.
  273. goto 9999
  274. endif
  275. T_bcj = p_bcj / (Rgas_b * r_bcj)
  276. d_bcj = k0 + u_acj
  277. u_bcj = d_bcj - (m/r_bcj)
  278. CC
  279. CC**** Test. We write the values thus determined
  280. CC
  281. C write(*,*) 'T_acj, p_acj, r_acj, u_acj, d_acj '
  282. C write(*,*) T_acj, p_acj, r_acj, u_acj, d_acj
  283. C write(*,*) 'T_bcj, p_bcj, r_bcj, u_bcj, d_bcj'
  284. C write(*,*) T_bcj, p_bcj, r_bcj, u_bcj, d_bcj
  285. C stop
  286. C
  287. if (d_bcj .gt. d_acj) then
  288. if (logdeb) then
  289. write(*,*) 'subroutine funcj'
  290. write(*,*) 'Deflagration regime'
  291. write(*,*) 'd_bcj, d_acj ', d_bcj, d_acj
  292. write(*,*)
  293. $ 'Reactive shock faster than the non-reactive one'
  294. endif
  295. C logan = .true.
  296. C goto 9999
  297. endif
  298. endif
  299. C write(*,*) 'k0cjdt, d_acj, d_bcj ', k0cjdt, d_acj, d_bcj
  300. C
  301. 9999 continue
  302. return
  303. end
  304. C
  305.  
  306.  

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