Télécharger racnre.eso

Retour à la liste

Numérotation des lignes :

  1. C RACNRE SOURCE BECC 09/12/07 21:15:59 6579
  2. subroutine racnre(nordpo, Tmaxcv,
  3. & Rgas_l, acv_l, Rgas_r, acv_r,
  4. & r_l, p_l, u_l, T_l, gam_l,
  5. & r_r, p_r, u_r, T_r, gam_r,
  6. & d_l, d_r,
  7. & r_b, p_b, u_b, T_b, d_b,
  8. & r_a, p_a, u_a, T_a, d_a,
  9. & logan, lognc)
  10. C
  11. c************************************************************************
  12. c
  13. c project :
  14. c
  15. c name : racnre
  16. c
  17. c description : computation of the shock-shock field-by-field
  18. c decomposition of the non-reactive Riemann problem
  19. C for thermally perfect gases.
  20. C
  21. C We have to distinguish between 4 different states
  22. C
  23. C l = left state
  24. C |
  25. C | Left GNL wave
  26. C |
  27. C b = left star state;
  28. C connected to l by a GNL wave
  29. C |
  30. C | Contact discontinuity
  31. C |
  32. C a = right star state;
  33. C connected to r by a GNL wave
  34. C |
  35. C | Right GNL wave
  36. C |
  37. C r = right state;
  38. C unburnt state
  39. c
  40. c language : fortran 77
  41. c
  42. c author : a. beccantini den/dm2s/sfme/ltmf
  43. c
  44. c************************************************************************
  45. c
  46. c called by : ...
  47. c
  48. c called program : ...
  49. c
  50. c************************************************************************
  51. c
  52. c
  53. c**** input
  54. c
  55. c nordpo = polynomials degree of cv(t)
  56. c tmaxcv = threshold temperature for cv(t)
  57. c rgas_l = gas constant (J/kg/K)in l (and in ls)
  58. c acv_l = cv_l = \sum_{i=0,nordpo} (acv_l{i} T^i)
  59. c rgas_r = gas constant r (and in rs)
  60. c acv_r = cv_r = \sum_{i=0,nordpo} (acv_r{i} T^i)
  61. c
  62. C r_l = density in l
  63. C p_l = pressure in
  64. c u_l = velocity in l
  65. C T_l = temperature in l (redundant with r_l and
  66. c p_l)
  67. c gam_l = specific heat ratio in l (redundant
  68. c with r_l and p_l)
  69. c
  70. C r_r = density in r
  71. C p_r = pressure in r
  72. c u_r = velocity in r
  73. C T_r = temperature in r (redundant with r_r and
  74. c p_r)
  75. c gam_r = specific heat ratio in r (redundant
  76. c with r_r and p_r)
  77. c
  78. c**** output
  79. C
  80. C d_r, d_l = speed of states l and r
  81. c
  82. c
  83. C r_b = density in b
  84. C p_b = pressure in b
  85. c u_b = velocity in b
  86. c T_b = temperature in b
  87. c d_b = shock speed of the left shock l->b
  88. c
  89. C r_a = density in a
  90. C p_a = pressure in a
  91. c u_a = velocity in a
  92. c T_a = temperature in a
  93. c d_a = shock speed of the
  94. C right shock a->r
  95. c
  96. c logan = an anomaly is detected
  97. c lognc = the approach has some problems
  98. C
  99. c
  100. c************************************************************************
  101. c
  102. c 07/12/2009 Created.
  103. c
  104. c implicit none
  105. integer nordpo
  106. integer nmax
  107. parameter (nmax = 20)
  108. integer iter, icon
  109. real*8 eps
  110. parameter(eps = 1.0D-6)
  111. real*8 Tmaxcv, Rgas_l, Rgas_r
  112. & , acv_l(0:nordpo), acv_r(0:nordpo)
  113. & , r_l, p_l, u_l, T_l, gam_l
  114. & , r_r, p_r, u_r, T_r, gam_r
  115. & , d_l, d_r
  116. & , r_b, p_b, u_b, T_b, d_b
  117. & , r_a, p_a, u_a, T_a, d_a
  118. C
  119. real*8 Tref, pref, Dtmin, coefder
  120. real*8 T_ag, T_bg
  121. real*8 P_ag, r_ag, u_ag, d_ag
  122. real*8 T_agp, P_agp, r_agp, u_agp, d_agp
  123. real*8 dpdua
  124. real*8 P_bg, r_bg, u_bg, d_bg
  125. real*8 T_bgp, P_bgp, r_bgp, u_bgp, d_bgp
  126. real*8 dpdub
  127. real*8 a, b
  128. real*8 erro
  129. logical lognc, logan, logdeb
  130. parameter(logdeb = .false.)
  131. CC
  132. CC**** We write the inputs
  133. CC
  134. if( logdeb ) then
  135. write(*,*) ' '
  136. write(*,*) 'RP'
  137. write(*,*) 'Left state'
  138. write(*,*) 'R, cv'
  139. write(*,'(6E12.4)') rgas_l, (acv_l(icon), icon=0,nordpo)
  140. write(*,*) 'rho, p, u, T, gamma'
  141. write(*,'(6E12.4)') r_l, p_l, u_l, T_l, gam_l
  142. write(*,*) 'Right state'
  143. write(*,*) 'R, cv'
  144. write(*,'(6E12.4)') rgas_r, (acv_r(icon), icon=0,nordpo)
  145. write(*,*) 'rho, p, u, T, gamma'
  146. write(*,'(6E12.4)') r_r, p_r, u_r, T_r, gam_r
  147. endif
  148. C stop
  149. C
  150. C**** Reference scales for a compressible solver.
  151. C For a low Mach number solver, the cut off speed is better than
  152. C the sound speed
  153. C
  154. Tref = 0.5d0 * ((T_l + T_r) + abs (T_l - T_r))
  155. Dtmin = 1.0D-3 * Tref
  156. coefder = 1.001D0
  157. pref = 0.5d0 * ((p_l + p_r) + abs (p_l - p_r))
  158. C
  159. C Initialisation of temperatures in a, b, c, d
  160. C
  161. T_a = T_r
  162. T_b = T_l
  163. C
  164. do iter = 1, nmax
  165. C
  166. C Guess temperatures at this iteration
  167. C
  168. T_ag = T_a
  169. T_bg = T_b
  170. C
  171. C******* Right state
  172. C
  173. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  174. & P_r, T_r, u_r, .true.,
  175. & T_ag, P_ag, r_ag, u_ag, d_ag)
  176. C
  177. T_agp = max((T_ag + Dtmin), (T_ag * coefder))
  178. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  179. & P_r, T_r, u_r, .true.,
  180. & T_agp, P_agp, r_agp, u_agp, d_agp)
  181. C
  182. dpdua = (P_agp - P_ag) / (u_agp - u_ag)
  183. if (dpdua .lt. 0.0d0) then
  184. write(*,*) 'dpdua > 0 , ', dpdua
  185. logan = .true.
  186. write(*,*) 'subroutine racnre.f'
  187. write(*,*) 'an error has been detected'
  188. goto 9999
  189. endif
  190. C
  191. C******* Left state
  192. C
  193. call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
  194. & P_l, T_l, u_l, .false.,
  195. & T_bg, P_bg, r_bg, u_bg, d_bg)
  196. C
  197. T_bgp = max((T_bg + Dtmin), (T_bg * coefder))
  198. call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
  199. & P_l, T_l, u_l, .false.,
  200. & T_bgp, P_bgp, r_bgp, u_bgp, d_bgp)
  201. C
  202. dpdub = (P_bgp - P_bg) / (u_bgp - u_bg)
  203. if (dpdub .gt. 0.0d0) then
  204. write(*,*) 'dpdub > 0 , ', dpdub
  205. logan = .true.
  206. write(*,*) 'subroutine racnre.f'
  207. write(*,*) 'an error has been detected'
  208. goto 9999
  209. endif
  210. C
  211. C******* Intersection
  212. C
  213. a = dpdua - dpdub
  214. b = ((dpdua * u_ag) - (dpdub * u_bg)) + (p_bg - p_ag)
  215. u_b = b / a
  216. p_b = p_bg + (dpdub * (u_b - u_bg))
  217. p_a = p_ag + (dpdua * (u_b - u_ag))
  218. if ( abs (p_a - p_b) .gt. (eps * pref)) then
  219. write(*,*) 'pa - pb too big, ', (p_a - p_b)
  220. logan = .true.
  221. write(*,*) 'subroutine racnre.f'
  222. write(*,*) 'an error has been detected'
  223. goto 9999
  224. endif
  225. C p_a, p_b should be positive
  226. p_a = max (p_a, (eps * p_ag))
  227. p_b = max (p_b, (eps * p_bg))
  228. C
  229. C******* Computation of temperature T_b, T_a
  230. C
  231. T_b = T_bg + ((p_b - p_bg) * (T_bgp - T_bg) / (p_bgp - p_bg))
  232. T_b = max (T_b, (eps * T_bg))
  233. C
  234. T_a = T_ag + ((p_a - p_ag) * (T_agp - T_ag) / (p_agp - p_ag))
  235. T_a = max (T_a, (eps * T_ag))
  236. C
  237. C******* TEST
  238. C
  239. erro = abs (T_a - T_ag) + abs (T_b - T_bg)
  240. erro = erro / Tref
  241. if (logdeb) then
  242. write(*,*) 'Intersection'
  243. write(*,*) 'T_ag, T_bg = ', T_ag, T_bg
  244. write(*,*) 'dpdua, dpdub = ', dpdua, dpdub
  245. write(*,*) 'p_a, p_b = ', p_a, p_b
  246. write(*,*) 'T_a, T_b = ', T_a, T_b
  247. WRITE(*,*) 'erro =', erro
  248. endif
  249. if (erro .le. 1.0D-14) then
  250. lognc = .false.
  251. goto 9998
  252. endif
  253. enddo
  254. if (logdeb ) then
  255. write(*,*) 'Warning. Convergence non-reached'
  256. write(*,*) 'erro = ', erro
  257. endif
  258. lognc = .true.
  259. 9998 continue
  260. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  261. & P_r, T_r, u_r, .true.,
  262. & T_a, P_a, r_a, u_a, d_a)
  263. call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
  264. & P_l, T_l, u_l, .false.,
  265. & T_b, P_b, r_b, u_b, d_b)
  266. d_l = gam_l * P_l / r_l
  267. d_l = u_l - sqrt(d_l)
  268. d_l = min(d_l, d_b)
  269. d_r = gam_r * P_r / r_r
  270. d_r = u_r + sqrt(d_r)
  271. d_r = max(d_r, d_a)
  272. C
  273. if (logdeb) then
  274. write(*,*) 'Final result'
  275. write(*,*) 'iter = ', iter
  276. write(*,*) 'r. rho, p, u, T, d'
  277. write(*,*) r_r, p_r, u_r, T_r, d_r
  278. write(*,*) 'a. rho, p, u, T, d'
  279. write(*,*) r_a, p_a, u_a, T_a, d_a
  280. write(*,*) 'b. rho, p, u, T, d'
  281. write(*,*) r_b, p_b, u_b, T_b, d_b
  282. write(*,*) 'l. rho, p, u, T, d'
  283. write(*,*) r_l, p_l, u_l, T_l, d_l
  284. C write(*,*) 'acv'
  285. C do iter = 0, nordpo, 1
  286. C write(*,*) acv_b(iter), acv_l(iter)
  287. C enddo
  288. endif
  289. 9999 continue
  290. return
  291. end
  292.  
  293.  

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