Télécharger flusti.eso

Retour à la liste

Numérotation des lignes :

  1. C FLUSTI SOURCE BECC 11/05/26 21:15:25 6981
  2. SUBROUTINE FLUSTI(
  3. & PC_L, GAM_L,
  4. & PC_R, GAM_R,
  5. & Z,
  6. & D_L, D_B, D_A, D_R,
  7. & RHO_L, P_L, U_L, UT1_L, UT2_L,
  8. & RHO_B, P_B, U_B,
  9. & RHO_A, P_A, U_A,
  10. & RHO_R, P_R, U_R, UT1_R, UT2_R,
  11. & UINT,
  12. & FLURHO, FLURU, FLURT1, FLURT2, FLURET,
  13. & LOGDEB, LOGAN)
  14. *
  15. *************************************************************************
  16. *
  17. * project : CAST3M, EUROPLEXUS...
  18. *
  19. * name : flusti
  20. *
  21. * description : euler equations for a mixture of stiffened gases
  22. * flux in the non-reactive case.
  23. *
  24. * given z = x/t, we compute the flux
  25. * f(x/t) - ((x/t) * u)
  26. * we compute density, pressure and velocity by
  27. * averaging the intermediate states using the
  28. * weve speeds.
  29. * we compute the other variables using these
  30. * quantities and the eos
  31. *
  32. * language : fortran 77
  33. *
  34. * author : a. beccantini den/dm2s/sfme/ltmf
  35. *
  36. *************************************************************************
  37. *
  38. * called by :
  39. *
  40. *
  41. *************************************************************************
  42. *
  43. ***** input
  44. *
  45. * pc_l, gam_l = properties of the gas in the left
  46. *
  47. * pc_r, gam_r = properties of the gas in the right
  48. *
  49. * z = x / t
  50. *
  51. * d_l, d_b, d_a, d_r =
  52. * = wave speeds
  53. *
  54. * rho_, p_, u_
  55. * = density, pressure, velocity
  56. *
  57. * uint = interface speed
  58. *
  59. * logdeb = debugging ?
  60. *
  61. ***** output
  62. *
  63. * flurho, flurhou, flurt1, flurt2, fluret =
  64. * fluxes for rho, rho u, rho ut1, rhout_2, rho e_t
  65. *
  66. * logan = anomaly
  67. *
  68. *************************************************************************
  69. *
  70. * 26/11/2009 created
  71. * 25/05/2011 evolution in CAST3M
  72. *
  73. *************************************************************************
  74. *
  75. * n.b.: all variables are declared
  76. *
  77. C IMPLICIT NONE
  78. IMPLICIT INTEGER(I-N)
  79. REAL*8
  80. & PC_L, GAM_L,
  81. & PC_R, GAM_R,
  82. & Z,
  83. & D_L, D_B, D_A, D_R,
  84. & RHO_L, P_L, U_L, UT1_L, UT2_L,
  85. & RHO_B, P_B, U_B,
  86. & RHO_A, P_A, U_A,
  87. & RHO_R, P_R, U_R, UT1_R, UT2_R,
  88. & UINT,
  89. & FLURHO, FLURU, FLURT1, FLURT2, FLURET
  90. *
  91. REAL*8 FP_L, FP_B, FP_A, FP_R
  92. & , GAM_Z, PC_Z
  93. & , RHO_Z, P_Z, U_Z, E_Z, UT1_Z, UT2_Z
  94. & , ECIN_Z
  95. *
  96. * debugging ?
  97. *
  98. LOGICAL LOGDEB, LOGAN
  99. *
  100. ******** states
  101. *
  102. * l
  103. * ***
  104. * |* r
  105. * | * b *********
  106. * | ******** *|
  107. * | | | * a * |
  108. * | | | ************* |
  109. * | | | | | |
  110. * | | | | | |
  111. * d_l d_b un_b un_a d_a d_r
  112. *
  113. *c
  114. *c test
  115. *c
  116. * z = 1.13d0
  117. * d_l = -1.0d0
  118. * d_b = 1.10d0
  119. * d_a = 1.11d0
  120. * d_r = 1.12d0
  121. *
  122. * speed of the contact discontinuty u_ab = u_a = u_b
  123. * if there is no vacuum between a and b
  124. *
  125. *
  126. ***** weight functions
  127. *
  128. FP_L = 0.0D0
  129. FP_B = 0.0D0
  130. FP_A = 0.0D0
  131. FP_R = 0.0D0
  132. *
  133. * this can be simplified from a computational point of view
  134. * using max, min, abs...
  135. * but in this way it is easy to read it...
  136. *
  137. IF (Z .LT. D_L) THEN
  138. FP_L = 1.0D0
  139. ELSEIF(Z .LT. D_B) THEN
  140. * rarefaction if d_l < d_b
  141. FP_L = (D_B - Z) / (D_B - D_L)
  142. FP_B = 1.0D0 - FP_L
  143. ELSEIF(Z .LT. U_B)THEN
  144. * contact discontinuty, left
  145. FP_B = 1.0D0
  146. ELSEIF(Z .LT. U_A)THEN
  147. * vacuum (u_a < u_b)
  148. FP_B = (U_A - Z) / (U_A - U_B)
  149. FP_A = 1.0D0 - FP_B
  150. ELSEIF(Z .LT. D_A)THEN
  151. * rs, right
  152. FP_A = 1.0D0
  153. ELSEIF(Z .LT. D_R)THEN
  154. * rgnl rarefaction if d_r < d_a
  155. FP_A = (D_R - Z) / (D_R - D_A)
  156. FP_R = 1.0D0 - FP_A
  157. ELSE
  158. * rgnl, right
  159. FP_R = 1.0D0
  160. ENDIF
  161. *
  162. ***** shock-shock + (e.f) intermediate solution
  163. * or
  164. * er-parameterisation, with linearisation over
  165. * the rarefaction for p, rho, w
  166. *
  167. * nb expression of gam_z and pc_z are exact except
  168. * in the case of vacuum (u_b < z < u_a).
  169. * indeed
  170. * gam_z = gam_l if z < un_b
  171. * gam_z = gam_r if z > un_a
  172. *
  173. GAM_Z = ((FP_L + FP_B) * GAM_L) +
  174. & ((FP_A + FP_R) * GAM_R)
  175. PC_Z = ((FP_L + FP_B) * PC_L) +
  176. & ((FP_A + FP_R) * PC_R)
  177. UT1_Z = ((FP_L + FP_B) * UT1_L) +
  178. & ((FP_A + FP_R) * UT1_R)
  179. UT2_Z = ((FP_L + FP_B) * UT2_L) +
  180. & ((FP_A + FP_R) * UT2_R)
  181. RHO_Z = (FP_L * RHO_L) + (FP_B * RHO_B) +
  182. & (FP_A * RHO_A) + (FP_R * RHO_R)
  183. P_Z = (FP_L * P_L) + (FP_B * P_B) +
  184. & (FP_A * P_A) + (FP_R * P_R)
  185. U_Z = (FP_L * U_L) + (FP_B * U_B) +
  186. & (FP_A * U_A) + (FP_R * U_R)
  187. *
  188. IF (LOGDEB) THEN
  189. WRITE(*,*) 'FP_L, FP_B, FP_A, FP_R ', FP_L, FP_B, FP_A, FP_R
  190. WRITE(*,*) 'GAM_Z ', GAM_Z
  191. WRITE(*,*) 'PC_Z ', PC_Z
  192. WRITE(*,*) 'UT1_Z ', UT1_Z
  193. WRITE(*,*) 'UT2_Z ', UT2_Z
  194. WRITE(*,*) 'RHO_Z ', RHO_Z
  195. WRITE(*,*) 'P_Z ', P_Z
  196. WRITE(*,*) 'U_Z ', U_Z
  197. ENDIF
  198. *
  199. ***** computation of the internal energy
  200. *
  201. CALL ESTIFF(RHO_Z, P_Z, GAM_Z, PC_Z, E_Z, LOGDEB, LOGAN)
  202. IF (LOGAN) THEN
  203. WRITE(*,*) 'SUBROUTINE FLUSTI.F'
  204. WRITE(*,*) 'ANOMALY DETECTED'
  205. GOTO 9999
  206. ENDIF
  207. ECIN_Z = 0.5D0 * ((U_Z * U_Z) + (UT1_Z * UT1_Z) +
  208. & (UT2_Z * UT2_Z))
  209. *
  210. ***** interfacial flux
  211. * according to nkonga, comput methods appl. mech engnr 190, 2000
  212. * \dep{u}{x} + \dep{f(u)}{x} = 0
  213. * z = speed of the surface
  214. * f_z = f - z u
  215. *
  216. UINT = U_Z - Z
  217. FLURHO = RHO_Z * (U_Z - Z)
  218. * write(*,*) 'flurho', flurho
  219. FLURU = (FLURHO * U_Z) + P_Z
  220. FLURET = (FLURHO * (E_Z + ECIN_Z))
  221. & + (P_Z * U_Z)
  222. FLURT1 = FLURHO * UT1_Z
  223. FLURT2 = FLURHO * UT2_Z
  224. *
  225. * write(*,*) (et_z + ef_z +
  226. * & ecin_z), flurho
  227. * write(*,*) p_z, u_z
  228. 9999 RETURN
  229. END
  230.  
  231.  

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