Télécharger flufu2.eso

Retour à la liste

Numérotation des lignes :

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

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