Télécharger ausmf.eso

Retour à la liste

Numérotation des lignes :

  1. C AUSMF SOURCE CB215821 16/04/21 21:15:10 8920
  2. c----------------------------------------------------------------------
  3. c GENERAL DESCRIPTION:
  4. c This subroutine provides the jacobian which is the derivatives
  5. c of the numerical flux function of the AUSM+ scheme
  6. c defined at the cell interface
  7. c with respect to the primitive variables of the left
  8. c cell (relative to the cell interface).
  9. c
  10. c EQUATIONS: 2D Euler equations of gas dynamics
  11. c
  12. c
  13. c REFERENCE: JCP, 129, 364-382 (1996)
  14. c " A Sequel to AUSM: AUSM+ ".
  15. c----------------------------------------------------------------------
  16. c INPUT:
  17. c
  18. c
  19. c (rog,uxg,uyg,pg) -- vector of the primitive variables at the
  20. c left cell;
  21. c
  22. c (rod,uxd,uyd,pd) -- vector of the primitive variables at the
  23. c right cell;
  24. c
  25. c (cnx,cny) -- normal vector to the interface (2 components in 2D);
  26. c
  27. c (ctx,cty) -- tangential vector to the interface;
  28. c
  29. c ga -- ratio of the specific heats (assumed constant)
  30. c----------------------------------------------------------------------
  31. c
  32. c OUTPUT:
  33. c
  34. c dfro -- vector of length 4 which represents derivatives of
  35. c flux-function for conservative variable (rho)
  36. c with respect to the primitive variables
  37. c
  38. c dfrun -- vector of length 4 which represents derivatives of
  39. c flux-function for conservative variable (rho u_n)
  40. c with respect to the primitive variables
  41. c
  42. c dfrut -- vector of length 4 which represents derivatives of
  43. c flux-function for conservative variable (rho u_t)
  44. c with respect to the primitive variables
  45. c
  46. c dfret -- vector of length 4 which represents derivatives of
  47. c flux-function for conservative variable (rho e_t)
  48. c with respect to the primitive variables
  49. c----------------------------------------------------------------------
  50. SUBROUTINE AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
  51. & ga,cnx,cny,ctx,cty,
  52. & dfro,dfrun,dfrut,dfret)
  53. IMPLICIT INTEGER(I-N)
  54. real*8 rog,uxg,uyg,pg
  55. real*8 rod,uxd,uyd,pd
  56. real*8 cnx,cny,ctx,cty
  57. real*8 jl(4,4),f(4)
  58. real*8 dfro(4),dfrun(4)
  59. real*8 dfrut(4),dfret(4)
  60. real*8 alpha,beta
  61. real*8 ga, gm1,temph
  62. real*8 n_x,n_y,t_x,t_y
  63. real*8 un_l, un_r, ut_l, ut_r
  64. real*8 ml,mr,Mplus,Mmin,mmid
  65. real*8 mpl_m, mmin_m,am
  66. real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  67. real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  68. real*8 Pplus,Pmin,pmid
  69. real*8 hr_l,hr_r
  70. real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
  71. real*8 aleft, arigh
  72. real*8 damr_l,damu_l
  73. real*8 damv_l,damp_l
  74. real*8 dmlr_l,dmlu_l
  75. real*8 dmlv_l,dmlp_l
  76. real*8 dmrr_l,dmru_l
  77. real*8 dmrv_l,dmrp_l
  78. real*8 dMpr_l,dMpu_l
  79. real*8 dMpv_l,dMpp_l
  80. real*8 dMmr_l,dMmu_l
  81. real*8 dMmv_l,dMmp_l
  82. real*8 dmir_l,dmiu_l
  83. real*8 dmiv_l,dmip_l
  84. real*8 d3mr_l,d3mu_l
  85. real*8 d3mv_l,d3mp_l
  86. real*8 d2mr_l,d2mu_l
  87. real*8 d2mv_l,d2mp_l
  88. real*8 dPpr_l,dPpu_l
  89. real*8 dPpv_l,dPpp_l
  90. real*8 dPmr_l,dPmu_l
  91. real*8 dPmv_l,dPmp_l
  92. real*8 dpir_l,dpiu_l
  93. real*8 dpiv_l,dpip_l
  94. integer i
  95. c
  96. c alpha -- parameter of the AUSM+ scheme in the Pressure function;
  97. c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  98. c
  99. c beta -- parameter of the AUSM+ scheme in the Mach function;
  100. c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
  101. c
  102. parameter(alpha = 0.1875D0, beta = 0.125D0)
  103. c-----------------------
  104. gm1=ga-1.0d0
  105. c-----------------------
  106. n_x=cnx
  107. n_y=cny
  108. t_x=ctx
  109. t_y=cty
  110. c----------------------------
  111. c----------------------------
  112. rold_l=rog
  113. uold_l=uxg
  114. vold_l=uyg
  115. pold_l=pg
  116. c-----------------------
  117. rold_r=rod
  118. uold_r=uxd
  119. vold_r=uyd
  120. pold_r=pd
  121. c------------------------------------------------------------------
  122. c Computation of the specific total energy on the left and right.
  123. c------------------------------------------------------------------
  124. eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  125. eold_l=eold_l+pold_l/(gm1*rold_l)
  126. eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  127. eold_r=eold_r+pold_r/(gm1*rold_r)
  128. c-------------------------------------------------------------------
  129. c Computation of the speed of sound and its derivatives;
  130. c numerical speed of sound at the interface is taken as an average
  131. c of the speeds of sounds of the neighbouring cells
  132. c-------------------------------------------------------------------
  133. aleft=SQRT(ga*pold_l/rold_l)
  134. arigh=SQRT(ga*pold_r/rold_r)
  135. am=0.5d0*(aleft+arigh)
  136. c-------------------------------------------------------------------
  137. damr_l=-aleft/(4.0d0*rold_l)
  138. damu_l=0.0d0
  139. damv_l=0.0d0
  140. damp_l=ga/(4.0d0*aleft*rold_l)
  141. c-------------------------------------------------------------------
  142. c Computing numerical Mach number and its derivatives,
  143. c see p.370, under (A1).
  144. c-------------------------------------------------------------------
  145. un_l=uold_l*n_x+vold_l*n_y
  146. un_r=uold_r*n_x+vold_r*n_y
  147. ut_l=uold_l*t_x+vold_l*t_y
  148. ut_r=uold_r*t_x+vold_r*t_y
  149. c-------------------------------------------------------------------
  150. ml=un_l/am
  151. mr=un_r/am
  152. c-------------------------------------------------------------------
  153. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  154. c see (19a) and (19b), p.367.
  155. c-------------------------------------------------------------------
  156. if(ABS(ml) .ge. 1.0d0) then
  157. Mplus=(ml+ABS(ml))/2.0d0
  158. else
  159. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  160. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  161. endif
  162. c-------------------------------------------------------------------
  163. if(ABS(mr) .ge. 1.0d0) then
  164. Mmin=(mr-ABS(mr))/2.0d0
  165. else
  166. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  167. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  168. endif
  169. c-------------------------------------------------------------------
  170. c Derivatives of ml and mr with respect to both sets of primitive
  171. c variables: left and right.
  172. c-------------------------------------------------------------------
  173. temp_l=-un_l/(am*am)
  174. temp_r=-un_r/(am*am)
  175. c--------
  176. dmlr_l=temp_l*damr_l
  177. dmrr_l=temp_r*damr_l
  178. c--------
  179. dmlu_l=n_x/am+temp_l*damu_l
  180. dmru_l=temp_r*damu_l
  181. c--------
  182. dmlv_l=n_y/am+temp_l*damv_l
  183. dmrv_l=temp_r*damv_l
  184. c--------
  185. dmlp_l=temp_l*damp_l
  186. dmrp_l=temp_r*damp_l
  187. c-----------------------------------------------------------
  188. c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
  189. c-----------------------------------------------------------
  190. mmid=Mplus+Mmin
  191. c-----------------------------------------------------------
  192. c Computing the derivatives of M+ and M-
  193. c-----------------------------------------------------------
  194. if(ml .ge. 1.0d0) then
  195. dMpr_l=dmlr_l
  196. dMpu_l=dmlu_l
  197. dMpv_l=dmlv_l
  198. dMpp_l=dmlp_l
  199. c--------------------
  200. else
  201. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  202. temph=(ml+1.0d0)/2.0d0
  203. dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  204. dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  205. dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  206. dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  207. c--------------------
  208. else
  209. dMpr_l=0.0d0
  210. dMpu_l=0.0d0
  211. dMpv_l=0.0d0
  212. dMpp_l=0.0d0
  213. c---------------------
  214. endif
  215. endif
  216. c-----------------------------------------------------------
  217. if(mr .ge. 1.0d0) then
  218. dMmr_l=0.0d0
  219. dMmu_l=0.0d0
  220. dMmv_l=0.0d0
  221. dMmp_l=0.0d0
  222. c---------------------
  223. else
  224. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  225. temph=(-mr+1.0d0)/2.0d0
  226. dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  227. dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  228. dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  229. dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  230. c--------------------
  231. else
  232. dMmr_l=dmrr_l
  233. dMmu_l=dmru_l
  234. dMmv_l=dmrv_l
  235. dMmp_l=dmrp_l
  236. c--------------------
  237. endif
  238. endif
  239. c-----------------------------------------------------------------
  240. c computing the derivatives of m_{1/2} (notation as in the paper)
  241. c-----------------------------------------------------------------
  242. dmir_l=dMpr_l+dMmr_l
  243. c-------------
  244. dmiu_l=dMpu_l+dMmu_l
  245. c-------------
  246. dmiv_l=dMpv_l+dMmv_l
  247. c-------------
  248. dmip_l=dMpp_l+dMmp_l
  249. c----------------------------------------------------------------
  250. c computing the main convective variables and their derivatives
  251. c mpl_m is m^{+}_{1/2} (paper's notation) and
  252. c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
  253. c----------------------------------------------------------------
  254. if(mmid .ge. 0.0d0) then
  255. mpl_m = mmid
  256. d2mr_l=dmir_l
  257. d2mu_l=dmiu_l
  258. d2mv_l=dmiv_l
  259. d2mp_l=dmip_l
  260. c------------
  261. else
  262. mpl_m = 0.0d0
  263. d2mr_l=0.0d0
  264. d2mu_l=0.0d0
  265. d2mv_l=0.0d0
  266. d2mp_l=0.0d0
  267. c------------
  268. endif
  269. c---------------------------------------------
  270. if(mmid .le. 0.0d0) then
  271. mmin_m = mmid
  272. d3mr_l=dmir_l
  273. d3mu_l=dmiu_l
  274. d3mv_l=dmiv_l
  275. d3mp_l=dmip_l
  276. c------------
  277. else
  278. mmin_m = 0.0d0
  279. d3mr_l=0.0d0
  280. d3mu_l=0.0d0
  281. d3mv_l=0.0d0
  282. d3mp_l=0.0d0
  283. c------------
  284. endif
  285. c---------------------------------------------------------------
  286. c Computing the calligraphic P+ and P- with their derivatives,
  287. c see (21a) & (21b) on p.368.
  288. c---------------------------------------------------------------
  289. if(ml .ge. 1.0d0) then
  290. Pplus = 1.0d0
  291. else
  292. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  293. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  294. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  295. else
  296. Pplus = 0.0d0
  297. endif
  298. endif
  299. c---------------------------------------------------------------
  300. if(mr .ge. 1.0d0) then
  301. Pmin = 0.0d0
  302. else
  303. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  304. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  305. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  306. else
  307. Pmin = 1.0d0
  308. endif
  309. endif
  310. c---------------------------------------------------------------
  311. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  312. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  313. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  314. c--------------
  315. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  316. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  317. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  318. c---------------------------------------------------------------
  319. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  320. dPpr_l=brac_l*dmlr_l
  321. dPpu_l=brac_l*dmlu_l
  322. dPpv_l=brac_l*dmlv_l
  323. dPpp_l=brac_l*dmlp_l
  324. else
  325. dPpr_l=0.0d0
  326. dPpu_l=0.0d0
  327. dPpv_l=0.0d0
  328. dPpp_l=0.0d0
  329. endif
  330. c---------------------------------------------------------------
  331. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  332. dPmr_l=brac_r*dmrr_l
  333. dPmu_l=brac_r*dmru_l
  334. dPmv_l=brac_r*dmrv_l
  335. dPmp_l=brac_r*dmrp_l
  336. else
  337. dPmr_l=0.0d0
  338. dPmu_l=0.0d0
  339. dPmv_l=0.0d0
  340. dPmp_l=0.0d0
  341. endif
  342. c-------------------------------------------------------------------
  343. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  344. c-------------------------------------------------------------------
  345. pmid=Pplus*pold_l + Pmin*pold_r
  346. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  347. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  348. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  349. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  350. c---------------------------------------------------------------------
  351. c computing JACOBIAN as a derivative of the numerical flux function
  352. c with respect to the primitive variables.
  353. c Notation: jl(i,j) --- is the derivative of the i-component of the
  354. c flux function with respect to the j-component of the
  355. c vector of primitive variables of the left state.
  356. c jr(i,j) --- is the derivative of the i-component of the
  357. c flux function with respect to the j-component of the
  358. c vector of primitive variables of the right state.
  359. c---------------------------------------------------------------------
  360. f(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
  361. c---------------------------------------------------------------------
  362. jl(1,1)=damr_l*f(1)/am+am*(d2mr_l*rold_l+mpl_m)
  363. jl(1,1)=jl(1,1)+am*d3mr_l*rold_r
  364. jl(1,2)=damu_l*f(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
  365. jl(1,3)=damv_l*f(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
  366. jl(1,4)=damp_l*f(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
  367. c------------------------------------
  368. c------------------------------------
  369. br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
  370. br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
  371. f(2)=n_x*(am*br1+pmid)+am*t_x*br2
  372. c---------------------------------------------------------
  373. br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
  374. br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
  375. jl(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
  376. jl(2,1)=jl(2,1)+t_x*damr_l*br2+am*t_x*br4
  377. c-------------------
  378. br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*n_x
  379. br3=br3+d3mu_l*rold_r*un_r
  380. br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*t_x
  381. br4=br4+d3mu_l*rold_r*ut_r
  382. jl(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
  383. jl(2,2)=jl(2,2)+t_x*damu_l*br2+am*t_x*br4
  384. >div rtqle="font: normal normal 1em/1.2em monospace; margin:0; padding:0; background:none; vertical-align:top;">c-------------------
  • br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*n_y
  • br3=br3+d3mv_l*rold_r*un_r
  • br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*t_y
  • br4=br4+d3mv_l*rold_r*ut_r
  • jl(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
  • jl(2,3)=jl(2,3)+t_x*damv_l*br2+am*t_x*br4
  • c-------------------
  • br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
  • br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
  • jl(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
  • jl(2,4)=jl(2,4)+t_x*damp_l*br2+am*t_x*br4
  • c-------------------------------------------------------------
  • c-------------------------------------------------------------
  • br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
  • br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
  • f(3)=n_y*(am*br1+pmid)+am*t_y*br2
  • br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
  • br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
  • jl(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
  • jl(3,1)=jl(3,1)+t_y*damr_l*br2+am*t_y*br4
  • c-------------------
  • br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*n_x
  • br3=br3+d3mu_l*rold_r*un_r
  • br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*t_x
  • br4=br4+d3mu_l*rold_r*ut_r
  • jl(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
  • jl(3,2)=jl(3,2)+t_y*damu_l*br2+am*t_y*br4
  • c-------------------
  • br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*n_y
  • br3=br3+d3mv_l*rold_r*un_r
  • br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*t_y
  • br4=br4+d3mv_l*rold_r*ut_r
  • jl(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
  • jl(3,3)=jl(3,3)+t_y*damv_l*br2+am*t_y*br4
  • c-------------------
  • br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
  • br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
  • jl(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
  • jl(3,4)=jl(3,4)+t_y*damp_l*br2+am*t_y*br4
  • c-------------------------------------------------------------
  • c ------ f44444444444444444444444444444 ---------
  • c-------------------------------------------------------------
  • hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+ga*pold_l/gm1
  • hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+ga*pold_r/gm1
  • f(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
  • c---------------------
  • br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
  • br1=br1+d3mr_l*hr_r
  • jl(4,1)=damr_l*f(4)/am+am*br1
  • c---------------------
  • br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
  • br1=br1+d3mu_l*hr_r
  • jl(4,2)=damu_l*f(4)/am+am*br1
  • c---------------------
  • br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
  • br1=br1+d3mv_l*hr_r
  • jl(4,3)=damv_l*f(4)/am+am*br1
  • c---------------------
  • br1=d2mp_l*hr_l+mpl_m*ga/gm1
  • br1=br1+d3mp_l*hr_r
  • jl(4,4)=damp_l*f(4)/am+am*br1
  • c----------------------------------------------------------
  • c----------------------------------------------------------
  • do 1 i=1,4
  • dfro(i)=jl(1,i)
  • dfrun(i)=jl(2,i)*n_x+jl(3,i)*n_y
  • dfrut(i)=jl(2,i)*t_x+jl(3,i)*t_y
  • dfret(i)=jl(4,i)
  • 1 continue
  • c----------------------------------------------------------------------
  • return
  • end
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales