Télécharger funre1.eso

Retour à la liste

Numérotation des lignes :

  1. C FUNRE1 SOURCE BECC 09/12/07 21:15:31 6579
  2. subroutine funre1(nordpo, tmaxcv, rgas_a, acv_a,
  3. & rgas_b, acv_b, q, k0, T_ag, T_bg,
  4. & r_a, p_a, u_a, T_a,
  5. & r_b, p_b, u_b, T_b, Dsw, logan)
  6. c************************************************************************
  7. c
  8. c projet :
  9. c
  10. c name : funre1
  11. c
  12. c description : computation of the weak state behind a reactive
  13. c shock
  14. c
  15. c language : fortran 77
  16. c
  17. c author : a. beccantini den/dm2s/sfme/ltmf
  18. c
  19. c************************************************************************
  20. c
  21. c callen by : racre2
  22. c
  23. c called program : none
  24. c
  25. c************************************************************************
  26. c
  27. c
  28. c**** input
  29. c
  30. c nordpo = polynomials degree
  31. c tmaxcv = threshold temperature for cv(t)
  32. c rgas_a = r of the gas in r and rs
  33. c acv_a = \sum_{i=1,nesp} y_i*acv_a{i,j}
  34. c rgas_b = r of the gas in rss
  35. c acv_b = \sum_{i=1,nesp} y_i*acv_b{i,j}
  36. c q = realeased heat per mass unit
  37. c k0 = fundamental flame speed
  38. C T_ag,T_bg = guess temperatures for a and b
  39. C (to compute properties)
  40. C r_a = density in rs
  41. C p_a = pressure in rs
  42. c u_a = velocity in rs
  43. C T_a = temperature in rs
  44. c**** output
  45. c r_b = approximate density in rss (behind the shock)
  46. c p_b = approximate pressure in rss
  47. c u_b = approximate speed in rss
  48. C T_b = approximate temperature in rss
  49. C (NB if T_bg = T_b and T_ag = T_a, rss thus
  50. C computed is exact)
  51. C Dsw = approximate shock wave speed... it could
  52. C differs from u_a + k0 if (k0lim < k0)
  53. C Nevertheless, if logsto = .true., if k0lim
  54. C < k0 we stop!
  55. c logan = an anomaly is detected
  56. c
  57. c************************************************************************
  58. c
  59. c created the 07/12/2009
  60. c
  61. c implicit none
  62. integer nordpo
  63. real*8 tmaxcv, rgas_a, acv_a(1:(nordpo+1))
  64. & ,rgas_b, acv_b(1:(nordpo+1))
  65. real*8 q, k0
  66. real*8 T_ag, T_bg
  67. real*8 r_a, p_a, u_a, T_a
  68. real*8 r_b, p_b, u_b, T_b, Dsw
  69. real*8 gam_a, e_a, cv_a, gam_b, e_b, cv_b
  70. real*8 ttq, k02ad, h_a
  71. real*8 a, b, c, delta, delta1, k0cj2, k0new
  72. real*8 csi
  73. logical logan
  74. CC
  75. CC**** We write the input
  76. CC
  77. C write(*,*) 'We enter funre1'
  78. C
  79. C**** Computation of gam_a and gam_b
  80. C
  81. call prith1(nordpo, acv_a, Tmaxcv, T_ag, e_a, cv_a)
  82. gam_a = (cv_a + Rgas_a) / cv_a
  83. call prith1(nordpo, acv_b, Tmaxcv, T_bg, e_b, cv_b)
  84. gam_b = (cv_b + Rgas_b) / cv_b
  85. ttq = q + ((cv_b * T_bg) - e_b) -
  86. & ((cv_a * T_ag) - e_a)
  87. C h_a and h_b computed as we deal with cp gases
  88. h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_a * T_a
  89. C
  90. r_a = P_a / (Rgas_a * T_a)
  91. C
  92. C**** Computation of r_b
  93. C
  94. k02ad = k0 * k0 / (Rgas_a * T_a)
  95. if (k02ad .lt. 1.0D0)then
  96. C if ( .true. ) then
  97. C if ( .false. ) then
  98. C write(*,*) 'funre1.f. Parte 1'
  99. c = k02ad
  100. b = k02ad + 1.0D0
  101. b = b * gam_b / (gam_b + 1.0D0)
  102. a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a)
  103. a = a + k02ad
  104. a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0)
  105. delta = (b * b) - (a * c)
  106. if (delta .lt. 0.0D0) then
  107. C if ( .true. ) then
  108. C write(*,*) 'funre1.f. Parte 1.2'
  109. C We are in CJDF regime
  110. a = 1.0D0 / (gam_b * gam_b)
  111. b = 1.0D0 - a
  112. b = b * (h_a + ttq)
  113. b = b - (Rgas_a * T_a)
  114. c = (Rgas_a * T_a)**2
  115. delta1 = (b * b) - (a * c)
  116. if (delta1 .lt. 0.0D0) then
  117. write(*,*) 'subroutine funre1.f'
  118. write(*,*) 'We cannot compute k0cj'
  119. write(*,*) 'ANOMALY DETECTED'
  120. logan = .true.
  121. goto 9999
  122. endif
  123. k0cj2 = (b - sqrt(delta1)) / a
  124. if (k0cj2 .lt. 0.0D0) then
  125. write(*,*) 'subroutine funre1.f'
  126. write(*,*) 'k0cj^2 < 0'
  127. write(*,*) 'ANOMALY DETECTED'
  128. logan = .true.
  129. goto 9999
  130. endif
  131. k02ad = k0cj2 / (Rgas_a * T_a)
  132. b = k02ad + 1.0D0
  133. b = b * gam_b / (gam_b + 1.0D0)
  134. a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a)
  135. a = a + k02ad
  136. a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0)
  137. csi = b / a
  138. k0new = k0cj2 ** 0.5D0
  139. else
  140. k0new = k0
  141. csi = b + sqrt(delta)
  142. csi = csi / a
  143. endif
  144. r_b = r_a * csi
  145. else
  146. C write(*,*) 'funre1.f. Parte 2'
  147. c = 2.0D0 * (h_a + ttq)
  148. c = c / (k0 * k0)
  149. c = c + 1.0D0
  150. c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
  151. b = Rgas_a * T_a / (k0 * k0)
  152. b = b + 1.0D0
  153. b = b * gam_b / (gam_b + 1.0D0)
  154. delta = (b * b) - c
  155. if (delta .lt. 0.0D0) then
  156. C if ( .true. ) then
  157. C write(*,*) 'funre1.f. Parte 2.2'
  158. a = 1.0D0 / (gam_b * gam_b)
  159. b = 1.0D0 - a
  160. b = b * (h_a + ttq)
  161. b = b - (Rgas_a * T_a)
  162. c = (Rgas_a * T_a)
  163. c = c * c
  164. delta1 = (b * b) - (a * c)
  165. if (delta1 .lt. 0.0D0) then
  166. write(*,*) 'subroutine funre1.f'
  167. write(*,*) 'We cannot compute k0cj'
  168. write(*,*) 'ANOMALY DETECTED'
  169. logan = .true.
  170. goto 9999
  171. endif
  172. k0cj2 = (b - sqrt(delta1)) / a
  173. if (k0cj2 .lt. 0.0D0) then
  174. write(*,*) 'subroutine funre1.f'
  175. write(*,*) 'k0cj^2 < 0'
  176. write(*,*) 'ANOMALY DETECTED'
  177. logan = .true.
  178. goto 9999
  179. endif
  180. c = 2.0D0 * (h_a + ttq)
  181. c = c / k0cj2
  182. c = c + 1.0D0
  183. c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
  184. b = Rgas_a * T_a / k0cj2
  185. b = b + 1.0D0
  186. b = b * gam_b / (gam_b + 1.0D0)
  187. delta1 = b * b - c
  188. C if (abs (delta1) > 1.0D-6) then
  189. if (abs (delta1) .gt. 1.0D-6) then
  190. write(*,*) 'subroutine funre1.f'
  191. write(*,*) 'CJDF'
  192. write(*,*) 'delta != 0'
  193. write(*,*) 'ANOMALY DETECTED'
  194. logan = .true.
  195. goto 9999
  196. endif
  197. csi = b
  198. k0new = k0cj2 ** 0.5D0
  199. else
  200. csi = b - sqrt (delta)
  201. k0new = k0
  202. endif
  203. r_b = r_a / csi
  204. endif
  205. csi = r_a / r_b
  206. Dsw = k0new + u_a
  207. u_b = Dsw - (k0new * csi)
  208. p_b = p_a + (r_a * k0new * k0new * (1.0d0 - csi))
  209. if (p_b .lt. 0.0D0) then
  210. write(*,*) 'subroutine funre1.f'
  211. write(*,*) 'p_b < 0'
  212. write(*,*) 'ANOMALY DETECTED'
  213. logan = .true.
  214. goto 9999
  215. endif
  216. T_b = P_b / (Rgas_b * r_b)
  217. C write(*,*) 'r_b, p_b, u_b, T_b, Dsw, k0new'
  218. C write(*,'(6E12.4)') r_b, p_b, u_b, T_b, Dsw
  219. C write(*,*) r_b, p_b, u_b, T_b, Dsw, k0new
  220. C stop
  221. CC write(*,*) 'We exit funre1'
  222. CC
  223. 9999 continue
  224. return
  225. end
  226.  
  227.  

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