Télécharger prithe.eso

Retour à la liste

Numérotation des lignes :

  1. C PRITHE SOURCE BECC 10/09/10 21:15:01 6472
  2. subroutine prithe(ndim,nesp,nLHS,nordpo,aprop,reacoe,Runiv,Tmaxcv
  3. $ ,w,gamma,Rgas,cvtot,acvtot,hftot,rho,u,p,T,y,ether,csi,logan)
  4. C
  5. C This subroutine creates some thermodynamical variables and
  6. C properties starting from the primitive variables.
  7. C
  8. C
  9. C
  10. C INPUT
  11. C *****
  12. C
  13. C ndim = dimension (1,2,3)
  14. C
  15. C nesp = number of species in the Euler equations
  16. C
  17. C nLHS = number of species in the LHS of the chemical reaction
  18. C nLHS < nesp
  19. C
  20. C nordpo = order of polynomial for cp and cv (see also Tmaxcv)
  21. C
  22. C aprop = gas properties: ,reacoe,Runiv,
  23. C aprop(0:nordpo+2,iesp) :
  24. C (coefcv(i),i=0,nordpo,1), Wi, hf0i)
  25. C
  26. C reacoe(1:nesp) = coefficient of each species in the global chemical
  27. C reaction
  28. C
  29. C Runiv = universal constant of the gas
  30. C
  31. C Tmaxcv = maximum temperature for cv polynomial expansion
  32. C cv(T) = cv(Tmaxcv) if T > Tmaxcv
  33. C
  34. C w = primitive variables
  35. C w(1) = rho
  36. C w(2) = ux
  37. C w(3) = uy
  38. C w(4) = uz
  39. C w(2+ndim) = p
  40. C w(3+ndim) = csi
  41. C w(3+ndim+1) = yf(1)
  42. C ...
  43. C w(3+ndim+nLHS) = yf(nLHS)
  44. C w(3+ndim+nLHS+1) = yi(nLHS+1)
  45. C ...
  46. C w(3+ndim+nesp-1) = yi(nesp-1)
  47. C w(3+ndim+nesp) = dy1t = yi(1) - yf(1)
  48. C w(4+ndim+nesp) = k0
  49. C
  50. C
  51. C
  52. C OUTPUT
  53. C ******
  54. C
  55. C gamma = specific heat ratio
  56. C
  57. C Rgas = gas constant
  58. C
  59. C cvtot = cv
  60. C
  61. C acvtot = coefficients such that for the mixture
  62. C cv = \sum_{i=1,nordpo+1} acvtot(i) T^{i-1}
  63. C ether = \sum_{i=1,nordpo+1} acvtot(i) T^{i} / (i)
  64. C
  65. C hftot = \sum_i Y_i hf0i = formation energy of the mixture at 0K
  66. C
  67. C rho = density
  68. C
  69. C u(1:ndim) = velocity
  70. C
  71. C p = pressure
  72. C
  73. C T = temperature
  74. C
  75. C y = mass fractions (1:nesp), \sum y(i) = 1
  76. C
  77. C ether = sensible energy at 0K
  78. C
  79. C csi = progress variable
  80. C
  81. C
  82. C ERROR VARIABLE
  83. C **************
  84. C
  85. C logan = if it becomes .false. an anomaly has been detected.
  86. C
  87. C************************************************************************
  88. C
  89. C HISTORY AND MODIFICATIONS
  90. C
  91. c 07/12/2009 Created.
  92. C
  93. C 30/07/2010 We check that species mass fractions are not lower
  94. C than yneg > 0 (and not zero).
  95. C
  96. C************************************************************************
  97. C
  98. C implicit none
  99. integer ndim, nesp, nLHS, nordpo
  100. & , iesp, icoef, idim
  101. real*8 aprop(0:(nordpo+2),nesp), reacoe(nesp), Runiv, Tmaxcv
  102. & , w(1:(4+ndim+nesp))
  103. & , gamma, Rgas, cvtot, acvtot(1:(nordpo+1)), hftot
  104. & , rho, u(1:ndim), p, T, y(1:nesp), ether
  105. & , csi, funt, dcv, coef, coef1, T1, yneg
  106. parameter (yneg = -1.0D-4)
  107. logical logan
  108. C
  109. C Initialisation of some variable
  110. C
  111. y(nesp)=1.0d0
  112. hftot = 0.0d0
  113. do icoef = 1 ,(nordpo+1), 1
  114. acvtot(icoef) = 0.0d0
  115. enddo
  116. Rgas = 0.0d0
  117. C
  118. C Density, velociy, pressure
  119. C
  120. rho = w(1)
  121. do idim = 1, ndim, 1
  122. u(idim) = w(1+idim)
  123. enddo
  124. p = w(2+ndim)
  125. C
  126. C Csi and mass fractions
  127. C
  128. csi = w(3+ndim)
  129. if ((csi . lt. 0.0d0) .or. (csi .gt. 1.0d0)) then
  130. write(*,*) 'prithe.f'
  131. write(*,*) 'csi = ', csi
  132. write(*,*) 'Negative csi.'
  133. logan = .true.
  134. goto 9999
  135. endif
  136. C coef = dy1t / (w_1 * a_1)
  137. coef = w(3+ndim+nesp) / (aprop(nordpo+1,1) * reacoe(1))
  138. do iesp = 1, nLHS, 1
  139. coef1 = aprop(nordpo+1,iesp) * reacoe(iesp)
  140. if (coef1 .le. 0.0D0) then
  141. write(*,*) 'prithe.f'
  142. write(*,*) 'Incompatibility in reacoe variable for ', iesp
  143. write(*,*) 'reacoe should be > 0. reacoe =', reacoe(iesp)
  144. logan=.true.
  145. goto 9999
  146. endif
  147. C y_j = y_{j,f} + ((1 - csi) * (y_{j,i} - y_{j,f}))
  148. y(iesp) = w(3+ndim+iesp) + ((1.0d0 - csi) * coef1 * coef)
  149. if(y(iesp) .lt. yneg)then
  150. write(*,*) 'prithe.f'
  151. write(*,*) 'Negative mass fraction, species ', iesp
  152. write(*,*) y(iesp)
  153. logan=.true.
  154. goto 9999
  155. endif
  156. y(nesp) = y(nesp) - y(iesp)
  157. enddo
  158. C
  159. do iesp = (nLHS + 1), (nesp - 1), 1
  160. coef1 = aprop(nordpo+1,iesp) * reacoe(iesp)
  161. if (coef1 .gt. 0.0D0) then
  162. write(*,*) 'prithe.f'
  163. write(*,*) 'Incompatibility in reacoe variable for ', iesp
  164. write(*,*) 'reacoe should be <= 0. reacoe =', reacoe(iesp)
  165. logan=.true.
  166. goto 9999
  167. endif
  168. C y_j = y_{j,i} + csi * (y_{j,f} - y_{j,i})
  169. y(iesp) = w(3+ndim+iesp) - (csi * coef1 * coef)
  170. if(y(iesp) .lt. yneg)then
  171. write(*,*) 'prithe.f'
  172. write(*,*) 'Negative mass fraction, species ', iesp
  173. write(*,*) y(iesp)
  174. logan=.true.
  175. goto 9999
  176. endif
  177. y(nesp) = y(nesp) - y(iesp)
  178. enddo
  179. C
  180. if(y(nesp) .lt. yneg)then
  181. write(*,*) 'prithe.f'
  182. write(*,*) 'Negative mass fraction, species ', nesp
  183. write(*,*) y(nesp)
  184. logan=.true.
  185. goto 9999
  186. endif
  187. C
  188. C Rgas, hftot, Temperature and acvtot
  189. C In order to avoid one more loop, this could have be done before
  190. C
  191. do iesp=1,nesp,1
  192. C Rgas = y(iesp) / W_iesp
  193. Rgas = Rgas + (y(iesp)/aprop(nordpo+1,iesp))
  194. C hftot = hftot + (y_i * hf0i)
  195. hftot = hftot + (y(iesp) * aprop(nordpo+2,iesp))
  196. do icoef = 0, nordpo, 1
  197. acvtot(icoef+1) = acvtot(icoef+1) +
  198. & (y(iesp) * aprop(icoef,iesp))
  199. enddo
  200. enddo
  201. Rgas = Rgas * Runiv
  202. T = p / (Rgas * rho)
  203. C
  204. C gamma, ether
  205. C
  206. C T1 = min (T, Tmaxcv)
  207. T1 = 0.5d0 * ((T + Tmaxcv) - abs (T - Tmaxcv))
  208. C
  209. funt = 1.0D0
  210. cvtot = acvtot(1)
  211. ether = cvtot * T1
  212. do icoef = 2, nordpo+1, 1
  213. funt = funt * T1
  214. dcv = acvtot(icoef) * funt
  215. cvtot = cvtot + dcv
  216. ether = ether + (dcv * T1 / icoef)
  217. enddo
  218. gamma = (Rgas + cvtot) / cvtot
  219. C
  220. C if ( T > Tmaxcv ) => T1 = (T - Tmaxcv)
  221. C else T1 = 0
  222. C i.e.
  223. C T1 = max ((T - Tmaxcv),0)
  224. C
  225. T1 = 0.5d0 * ((T - Tmaxcv) + abs (T - Tmaxcv))
  226. ether = ether + (cvtot * T1)
  227. C
  228. 9999 continue
  229. return
  230. end
  231. C
  232.  
  233.  
  234.  

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