Télécharger flufu1.eso

Retour à la liste

Numérotation des lignes :

  1. C FLUFU1 SOURCE BECC 09/12/07 21:15:21 6579
  2. subroutine flufu1(nordpo,
  3. & Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
  4. & Rgas_l, acv_l, acv_z,
  5. & z,
  6. & d_l, d_d, d_c, d_b, d_a, d_r,
  7. & rho_l, T_l, u_l, hf_l,
  8. & rho_d, T_d, u_d,
  9. & rho_c, T_c, u_c,
  10. & rho_b, T_b, u_b, hf_b,
  11. & rho_a, T_a, u_a,
  12. & rho_r, T_r, u_r, hf_r,
  13. & flurho, fluru, fluret)
  14. C
  15. c************************************************************************
  16. c
  17. c project :
  18. c
  19. c name : flufu1
  20. c
  21. c description : Reactive Euler equations for a mixture of reactive
  22. C thermally perfect gases.
  23. c
  24. c Given z = x/t, we compute the flux
  25. C F(x/t) - ((x/t) * U)
  26. C We compute density, temperature and velocity by
  27. C averaging the intermediate states using the
  28. C weve speeds.
  29. C We compute the other variables using these
  30. C quantities and the EOS
  31. c
  32. c language : fortran 77
  33. c
  34. c author : a. beccantini den/dm2s/sfme/ltmf
  35. c
  36. c************************************************************************
  37. c
  38. C called by : fbecre
  39. C
  40. c
  41. c************************************************************************
  42. c
  43. c**** input
  44. c
  45. c nordpo = polynomial order of the cv
  46. c
  47. C Tmaxcv = Tmax for cv expansion
  48. C
  49. C z = x / t
  50. C
  51. C Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
  52. C Rgas_l, acv_l = quantities to compute cv and cp
  53. C
  54. C acv_z = work vector
  55. C
  56. C d_l, d_d, un_d, un_c, d_c, d_b, d_a, d_r =
  57. C = wave speeds
  58. C
  59. C rho_, u_, T_, hf_
  60. C = density, velocity, temperature,
  61. C formation enthalpy
  62. C**** output
  63. C
  64. c flurho, flurhou, fluret = fluxes for rho, rho u, rho e_t
  65. C
  66. c
  67. c************************************************************************
  68. c
  69. C 07/12/2009 created
  70. C
  71. c************************************************************************
  72. c
  73. c n.b.: all variables are declared
  74. c
  75. c implicit none
  76. integer nordpo, iord
  77. real*8 Tmaxcv, Rgas_r, acv_r(0:nordpo)
  78. & , Rgas_b, acv_b(0:nordpo)
  79. & , Rgas_l, acv_l(0:nordpo)
  80. & , z, acv_z(0:nordpo)
  81. & , d_l, d_d, d_c, d_b, d_a, d_r
  82. & , rho_l, u_l, T_l, hf_l
  83. & , rho_d, T_d, u_d
  84. & , rho_c, T_c, u_c
  85. & , rho_b, T_b, u_b, hf_b
  86. & , rho_a, T_a, u_a
  87. & , rho_r, T_r, u_r, hf_r
  88. & , flurho, fluru, fluret
  89. real*8 u_cd, fp_l, fp_d, fp_c, fp_b, fp_a, fp_r
  90. & , Rgas_z, rho_z, T_z, u_z, p_z, e_z, cv_z, hf_z
  91. & , ecin_z
  92. logical logdeb
  93. parameter(logdeb = .false.)
  94. CC
  95. CC Test
  96. CC
  97. C z = 1.13D0
  98. C d_l = -1.0D0
  99. C d_d = 0.10D0
  100. C u_c = 0.2D0
  101. C u_d = 0.2D0
  102. C d_c = 1.09D0
  103. C d_b = 1.10D0
  104. C d_a = 1.11D0
  105. C d_r = 1.12D0
  106. C
  107. C Speed of the contact discontinuty u_cd = u_c = u_d
  108. C
  109. u_cd = 0.5D0 * (u_c + u_d)
  110. C
  111. C**** Weight functions
  112. C
  113. fp_l = 0.0d0
  114. fp_d = 0.0D0
  115. fp_c = 0.0d0
  116. fp_b = 0.0d0
  117. fp_a = 0.0d0
  118. fp_r = 0.0d0
  119. C
  120. C This can be simplified from a computational point of view
  121. C
  122. if (z .lt. d_l) then
  123. fp_l = 1.0d0
  124. elseif(z .lt. d_d) then
  125. C Rarefaction if d_l < d_d
  126. fp_l = (d_d - z) / (d_d - d_l)
  127. fp_d = 1.0d0 - fp_l
  128. elseif(z .lt. u_cd)then
  129. C Contact discontinuty, left
  130. fp_d = 1.0d0
  131. elseif(z .lt. d_c)then
  132. C Contact discontinuty, right
  133. fp_c = 1.0d0
  134. elseif(z .lt. d_b)then
  135. C Rarefaction if d_b < d_c
  136. fp_c = (d_b - z) / (d_b - d_c)
  137. fp_b = 1.0d0 - fp_c
  138. elseif(z .lt. d_a)then
  139. C RS, right
  140. fp_a = 1.0D0
  141. elseif(z .lt. d_r)then
  142. C Rgnl rarefaction if d_r < d_a
  143. fp_a = (d_r - z) / (d_r - d_a)
  144. fp_r = 1.0d0 - fp_a
  145. else
  146. C Rgnl, right
  147. fp_r = 1.0d0
  148. endif
  149. CC
  150. C
  151. C**** Shock-shock + (e.f) intermediate solution
  152. C
  153. Rgas_z = ((fp_l + fp_d) * Rgas_l) +
  154. & ((fp_c + fp_b) * Rgas_b) +
  155. & ((fp_a + fp_r) * Rgas_r)
  156. rho_z = (fp_l * rho_l) + (fp_d * rho_d) +
  157. & (fp_c * rho_c) + (fp_b * rho_b) +
  158. & (fp_a * rho_a) + (fp_r * rho_r)
  159. T_z = (fp_l * T_l) + (fp_d * T_d) +
  160. & (fp_c * T_c) + (fp_b * T_b) +
  161. & (fp_a * T_a) + (fp_r * T_r)
  162. p_z = Rgas_z * rho_z * T_z
  163. u_z = (fp_l * u_l) + (fp_d * u_d) +
  164. & (fp_c * u_c) + (fp_b * u_b) +
  165. & (fp_a * u_a) + (fp_r * u_r)
  166. if (logdeb) then
  167. write(*,*) 'd_l, d_d, d_c, d_b, d_a, d_r'
  168. write(*,'(6E12.4)') d_l, d_d, d_c, d_b, d_a, d_r
  169. write(*,*) 'fp_l, fp_d, fp_c, fp_b, fp_a, fp_r'
  170. write(*,'(6E12.4)') fp_l, fp_d, fp_c, fp_b, fp_a, fp_r
  171. write(*,*) 'State in x/t'
  172. write(*,*) 'Rgas_z, rho_z, T_z, p_z, u_z'
  173. write(*,*) Rgas_z, rho_z, T_z, p_z, u_z
  174. endif
  175. C
  176. C**** Computation of thermal energy
  177. C
  178. do iord = 0, nordpo , 1
  179. acv_z(iord) = ((fp_l + fp_d) * acv_l(iord)) +
  180. & ((fp_c + fp_b) * acv_b(iord)) +
  181. & ((fp_a + fp_r) * acv_r(iord))
  182. enddo
  183. C
  184. call prith1(nordpo, acv_z, Tmaxcv, T_z, e_z, cv_z)
  185. C
  186. hf_z = ((fp_l + fp_d) * hf_l) + ((fp_c + fp_b) * hf_b)
  187. & + ((fp_a + fp_r) * hf_r)
  188. ecin_z = 0.5d0 * u_z * u_z
  189. C
  190. C**** Interfacial flux
  191. C According to NKONGA, Comput Methods Appl. Mech Engnr 190, 2000
  192. C \dep{U}{x} + \dep{F(U)}{x} = 0
  193. C z = speed of the surface
  194. C f_z = F - z U
  195. C
  196. flurho = rho_z * (u_z - z)
  197. C write(*,*) 'flurho', flurho
  198. fluru = (flurho * u_z) + p_z
  199. fluret = (flurho * (e_z + hf_z +
  200. & ecin_z)) + (p_z * u_z)
  201. C write(*,*) (et_z + ef_z +
  202. C & ecin_z), flurho
  203. C write(*,*) p_z, u_z
  204. return
  205. end
  206.  
  207.  

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