Télécharger conjp5.eso

Retour à la liste

Numérotation des lignes :

  1. C CONJP5 SOURCE CB215821 16/04/21 21:15:59 8920
  2. c Voir konjp4.eso
  3. c----------------------------------------------------------------------
  4. c GENERAL DESCRIPTION:
  5. c This subroutine provides the jacobians which are the derivatives
  6. c of the numerical flux function defined at "wall" cell interface
  7. c with respect to the primitive variables of the left and right
  8. c cells (relative to the cell interface).
  9. c
  10. c EQUATIONS: 3D 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 alpha -- parameter of the AUSM+ scheme in the Pressure function;
  19. c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  20. c
  21. c beta -- parameter of the AUSM+ scheme in the Mach function;
  22. c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
  23. c
  24. c wvec_l -- vector of the primitive variables (rho,ux,uy,uz,p) at the
  25. c left cell;
  26. c
  27. c wvec_r -- vector of the primitive variables (rho,ux,uy,uz,p) at the
  28. c right cell;
  29. c
  30. c nvect -- normal vector to the interface (3 components in 3D);
  31. c
  32. c tvec1 -- tangential vector to the interface;
  33. c
  34. c tvec2 -- tangential vector to the interface;
  35. c
  36. c ga -- ratio of the specific heats (assumed constant)
  37. c----------------------------------------------------------------------
  38. c
  39. c OUTPUT:
  40. c
  41. c jl -- jakobian matrix 5 by 5 - derivatives of the numerical
  42. c flux function with respect to the primitive variables
  43. c from the left cell;
  44. c
  45. c jr -- jakobian matrix 5 by 5 - derivatives of the numerical
  46. c flux function with respect to the primitive variables
  47. c from the right cell.
  48. c----------------------------------------------------------------------
  49. SUBROUTINE CONJP5(JL,JR,WVEC_L,WVEC_R,NVECT,TVECT1,TVECT2,
  50. & ga)
  51. IMPLICIT INTEGER(I-N)
  52. real*8 wvec_l(5),wvec_r(5)
  53. real*8 nvect(3),tvect1(3),tvect2(3)
  54. real*8 jl(5,5),jr(5,5)
  55. real*8 alpha,beta
  56. real*8 ga, gm1
  57. real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
  58. real*8 un_l,un_r
  59. real*8 ml,mr,Mplus,Mmin
  60. real*8 am
  61. real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
  62. real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
  63. real*8 Pplus,Pmin
  64. real*8 temp_l,temp_r,brac_l,brac_r
  65. real*8 aleft, arigh
  66. real*8 damr_l,damr_r,damu_l,damu_r
  67. real*8 damv_l,damv_r,damp_l,damp_r
  68. real*8 damw_l,damw_r
  69. real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  70. real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  71. real*8 dmlw_l,dmlw_r
  72. real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  73. real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  74. real*8 dmrw_l,dmrw_r
  75. real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  76. real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  77. real*8 dPpw_l,dPpw_r
  78. real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  79. real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  80. real*8 dPmw_l,dPmw_r
  81. real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  82. real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  83. real*8 dpiw_l,dpiw_r
  84. real*8 c11,c12,c13,c21,c22,c23,c31,c32,c33,det
  85. real*8 zc11,zc12,zc13,zc21,zc22,zc23,zc31,zc32,zc33
  86. integer i
  87. parameter(alpha = 0.1875D0, beta = 0.125D0)
  88. c-----------------------
  89. gm1=ga-1.0d0
  90. c-----------------------
  91. n_x=nvect(1)
  92. n_y=nvect(2)
  93. n_z=nvect(3)
  94. c-----------------------
  95. t1_x=tvect1(1)
  96. t1_y=tvect1(2)
  97. t1_z=tvect1(3)
  98. c-----------------------
  99. t2_x=tvect2(1)
  100. t2_y=tvect2(2)
  101. t2_z=tvect2(3)
  102. c--------------------------------------
  103. rold_l=wvec_l(1)
  104. uold_l=wvec_l(2)
  105. vold_l=wvec_l(3)
  106. wold_l=wvec_l(4)
  107. pold_l=wvec_l(5)
  108. c-----------------------
  109. rold_r=wvec_r(1)
  110. uold_r=wvec_r(2)
  111. vold_r=wvec_r(3)
  112. wold_r=wvec_r(4)
  113. pold_r=wvec_r(5)
  114. c------------------------------------------------------------------
  115. c Computation of the specific total energy on the left and right.
  116. c------------------------------------------------------------------
  117. eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  118. eold_l=eold_l+pold_l/(gm1*rold_l)
  119. eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  120. eold_r=eold_r+pold_r/(gm1*rold_r)
  121. c-------------------------------------------------------------------
  122. c Computation of the speed of sound and its derivatives;
  123. c numerical speed of sound at the interface is taken as an average
  124. c of the speeds of sounds of the neighbouring cells
  125. c---------------------------------------------------------------------
  126. aleft=SQRT(ga*pold_l/rold_l)
  127. arigh=SQRT(ga*pold_r/rold_r)
  128. am=0.5d0*(aleft+arigh)
  129. c--------------------------------------------------------------------
  130. damr_r=-arigh/(4.0d0*rold_r)
  131. damu_r=0.0d0
  132. damv_r=0.0d0
  133. damw_r=0.0d0
  134. damp_r=ga/(4.0d0*arigh*rold_r)
  135. c-----------------------
  136. damr_l=-aleft/(4.0d0*rold_l)
  137. damu_l=0.0d0
  138. damv_l=0.0d0
  139. damw_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+wold_l*n_z
  146. un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
  147. c----------------------------------------
  148. ml=un_l/am
  149. mr=un_r/am
  150. c-------------------------------------------------------------------
  151. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  152. c see (19a) and (19b), p.367.
  153. c-------------------------------------------------------------------
  154. if(ABS(ml) .ge. 1.0d0) then
  155. Mplus=(ml+ABS(ml))/2.0d0
  156. else
  157. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  158. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  159. endif
  160. c-----------------
  161. if(ABS(mr) .ge. 1.0d0) then
  162. Mmin=(mr-ABS(mr))/2.0d0
  163. else
  164. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  165. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  166. endif
  167. c-------------------------------------------------------------------
  168. c Derivatives of ml and mr with respect to both sets of primitive
  169. c variables: left and right.
  170. c-------------------------------------------------------------------
  171. temp_l=-un_l/(am*am)
  172. temp_r=-un_r/(am*am)
  173. c--------
  174. dmlr_l=temp_l*damr_l
  175. dmlr_r=temp_l*damr_r
  176. dmrr_l=temp_r*damr_l
  177. dmrr_r=temp_r*damr_r
  178. c--------
  179. dmlu_l=n_x/am+temp_l*damu_l
  180. dmlu_r=temp_l*damu_r
  181. dmru_l=temp_r*damu_l
  182. dmru_r=n_x/am+temp_r*damu_r
  183. c--------
  184. dmlv_l=n_y/am+temp_l*damv_l
  185. dmlv_r=temp_l*damv_r
  186. dmrv_l=temp_r*damv_l
  187. dmrv_r=n_y/am+temp_r*damv_r
  188. c--------
  189. dmlw_l=n_z/am+temp_l*damw_l
  190. dmlw_r=temp_l*damw_r
  191. dmrw_l=temp_r*damw_l
  192. dmrw_r=n_z/am+temp_r*damw_r
  193. c--------
  194. dmlp_l=temp_l*damp_l
  195. dmlp_r=temp_l*damp_r
  196. dmrp_l=temp_r*damp_l
  197. dmrp_r=temp_r*damp_r
  198. c---------------------------------------------------------------
  199. c Computing the calligraphic P+ and P- with their derivatives
  200. c see (21a) & (21b) on p.368.
  201. c---------------------------------------------------------------
  202. if(ml .ge. 1.0d0) then
  203. Pplus = 1.0d0
  204. else
  205. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  206. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  207. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  208. else
  209. Pplus = 0.0d0
  210. endif
  211. endif
  212. c---------------------------------------------------------------
  213. if(mr .ge. 1.0d0) then
  214. Pmin = 0.0d0
  215. else
  216. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  217. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  218. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  219. else
  220. Pmin = 1.0d0
  221. endif
  222. endif
  223. c---------------------------------------------------------------
  224. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  225. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  226. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  227. c--------------
  228. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  229. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  230. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  231. c---------------------------------------------------------------
  232. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  233. dPpr_l=brac_l*dmlr_l
  234. dPpr_r=brac_l*dmlr_r
  235. c------------
  236. dPpu_l=brac_l*dmlu_l
  237. dPpu_r=brac_l*dmlu_r
  238. c------------
  239. dPpv_l=brac_l*dmlv_l
  240. dPpv_r=brac_l*dmlv_r
  241. c------------
  242. dPpw_l=brac_l*dmlw_l
  243. dPpw_r=brac_l*dmlw_r
  244. c------------
  245. dPpp_l=brac_l*dmlp_l
  246. dPpp_r=brac_l*dmlp_r
  247. c------------
  248. else
  249. dPpr_l=0.0d0
  250. dPpr_r=0.0d0
  251. c-----------
  252. dPpu_l=0.0d0
  253. dPpu_r=0.0d0
  254. c-----------
  255. dPpv_l=0.0d0
  256. dPpv_r=0.0d0
  257. c-----------
  258. dPpw_l=0.0d0
  259. dPpw_r=0.0d0
  260. c-----------
  261. dPpp_l=0.0d0
  262. dPpp_r=0.0d0
  263. c-----------
  264. endif
  265. c---------------------------------------------------------------
  266. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  267. dPmr_l=brac_r*dmrr_l
  268. dPmr_r=brac_r*dmrr_r
  269. c------------
  270. dPmu_l=brac_r*dmru_l
  271. dPmu_r=brac_r*dmru_r
  272. c------------
  273. dPmv_l=brac_r*dmrv_l
  274. dPmv_r=brac_r*dmrv_r
  275. c------------
  276. dPmw_l=brac_r*dmrw_l
  277. dPmw_r=brac_r*dmrw_r
  278. c------------
  279. dPmp_l=brac_r*dmrp_l
  280. dPmp_r=brac_r*dmrp_r
  281. c------------
  282. else
  283. dPmr_l=0.0d0
  284. dPmr_r=0.0d0
  285. c-----------
  286. dPmu_l=0.0d0
  287. dPmu_r=0.0d0
  288. c-----------
  289. dPmv_l=0.0d0
  290. dPmv_r=0.0d0
  291. c-----------
  292. dPmw_l=0.0d0
  293. dPmw_r=0.0d0
  294. c-----------
  295. dPmp_l=0.0d0
  296. dPmp_r=0.0d0
  297. c-----------
  298. endif
  299. c---------------------------------------------------------------------
  300. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  301. c---------------------------------------------------------------------
  302. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  303. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  304. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  305. dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
  306. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  307. c----------------------------
  308. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  309. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  310. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  311. dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
  312. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  313. c---------------------------------------------------------------------
  314. c computing JACOBIAN as a derivative of the numerical flux function
  315. c with respect to the primitive variables.
  316. c Notation: jl(i,j) --- is the derivative of the i-component of the
  317. c flux function with respect to the j-component of the
  318. c vector of primitive variables of the left state.
  319. c jr(i,j) --- is the derivative of the i-component of the
  320. c flux function with respect to the j-component of the
  321. c vector of primitive variables of the right state.
  322. c---------------------------------------------------------------------
  323. jl(1,1)=0.0D0
  324. jl(1,2)=0.0D0
  325. jl(1,3)=0.0D0
  326. jl(1,4)=0.0D0
  327. jl(1,5)=0.0D0
  328. c------------------------------------
  329. jr(1,1)=0.0D0
  330. jr(1,2)=0.0D0
  331. jr(1,3)=0.0D0
  332. jr(1,4)=0.0D0
  333. jr(1,5)=0.0D0
  334. c------------------------------------
  335. c------------------------------------
  336. jl(2,1)=dpir_l*n_x
  337. c------------------------
  338. jl(2,2)=dpiu_l*n_x
  339. c------------------------
  340. jl(2,3)=dpiv_l*n_x
  341. c------------------------
  342. jl(2,4)=dpiw_l*n_x
  343. c------------------------
  344. jl(2,5)=dpip_l*n_x
  345. c-------------------------------------------------
  346. c-------------------------------------------------
  347. jl(3,1)=dpir_l*n_y
  348. c------------------------
  349. jl(3,2)=dpiu_l*n_y
  350. c------------------------
  351. jl(3,3)=dpiv_l*n_y
  352. c------------------------
  353. jl(3,4)=dpiw_l*n_y
  354. c------------------------
  355. jl(3,5)=dpip_l*n_y
  356. c-------------------------------------------------
  357. c-------------------------------------------------
  358. jl(4,1)=dpir_l*n_z
  359. c------------------------
  360. jl(4,2)=dpiu_l*n_z
  361. c------------------------
  362. jl(4,3)=dpiv_l*n_z
  363. c------------------------
  364. jl(4,4)=dpiw_l*n_z
  365. c------------------------
  366. jl(4,5)=dpip_l*n_z
  367. c-------------------------------------------------------
  368. c derivatives with respect to the right
  369. c set of primitive variables
  370. c-------------------------------------------------------
  371. jr(2,1)=dpir_r*n_x
  372. c------------------------
  373. jr(2,2)=dpiu_r*n_x
  374. c------------------------
  375. jr(2,3)=dpiv_r*n_x
  376. c------------------------
  377. jr(2,4)=dpiw_r*n_x
  378. c------------------------
  379. jr(2,5)=dpip_r*n_x
  380. c-------------------------------------------------------
  381. jr(3,1)=dpir_r*n_y
  382. c------------------------
  383. jr(3,2)=dpiu_r*n_y
  384. c------------------------
  385. jr(3,3)=dpiv_r*n_y
  386. c------------------------
  387. jr(3,4)=dpiw_r*n_y
  388. c------------------------
  389. jr(3,5)=dpip_r*n_y
  390. c--------------------------------------------------------
  391. jr(4,1)=dpir_r*n_z
  392. c------------------------
  393. jr(4,2)=dpiu_r*n_z
  394. c------------------------
  395. jr(4,3)=dpiv_r*n_z
  396. c------------------------
  397. jr(4,4)=dpiw_r*n_z
  398. c------------------------
  399. jr(4,5)=dpip_r*n_z
  400. c-------------------------------------------------------
  401. c-------------------------------------------------------
  402. jl(5,1)=0.0D0
  403. jl(5,2)=0.0D0
  404. jl(5,3)=0.0D0
  405. jl(5,4)=0.0D0
  406. jl(5,5)=0.0D0
  407. c-------------------------------------------------
  408. c-------------------------------------------------
  409. jr(5,1)=0.0D0
  410. jr(5,2)=0.0D0
  411. jr(5,3)=0.0D0
  412. jr(5,4)=0.0D0
  413. jr(5,5)=0.0D0
  414. c-------------------------------------------------------------
  415. c-----------------------------------------------
  416. c Taking into account the dependancy of variables
  417. c-----------------------------------------------
  418. c11=t1_y*t2_z - t1_z*t2_y
  419. c12=n_y*t2_z - n_z*t2_y
  420. c13=n_y*t1_z - n_z*t1_y
  421. c-------------------------------------
  422. c21=t1_x*t2_z - t1_z*t2_x
  423. c22=n_x*t2_z - n_z*t2_x
  424. c23=n_x*t1_z - n_z*t1_x
  425. c-------------------------------------
  426. c31=t1_x*t2_y - t1_y*t2_x
  427. c32=n_x*t2_y - n_y*t2_x
  428. c33=n_x*t1_y - n_y*t1_x
  429. det=n_x*c11 - n_y*c21 + n_z*c31
  430. c----------------------------------------------------------------------
  431. ZC11=-NVECT(1)*C11-TVECT1(1)*C12+TVECT2(1)*C13
  432. ZC12=-NVECT(2)*C11-TVECT1(2)*C12+TVECT2(2)*C13
  433. ZC13=-NVECT(3)*C11-TVECT1(3)*C12+TVECT2(3)*C13
  434. C---------------------------------
  435. ZC21=NVECT(1)*C21+TVECT1(1)*C22-TVECT2(1)*C23
  436. ZC22=NVECT(2)*C21+TVECT1(2)*C22-TVECT2(2)*C23
  437. ZC23=NVECT(3)*C21+TVECT1(3)*C22-TVECT2(3)*C23
  438. C---------------------------------
  439. ZC31=-NVECT(1)*C31-TVECT1(1)*C32+TVECT2(1)*C33
  440. ZC32=-NVECT(2)*C31-TVECT1(2)*C32+TVECT2(2)*C33
  441. ZC33=-NVECT(3)*C31-TVECT1(3)*C32+TVECT2(3)*C33
  442. c---------------------------------------------------------------------
  443. do 11 i=1,5
  444. jl(i,1)=jl(i,1)+jr(i,1)
  445. jl(i,2)=jl(i,2)+jr(i,2)*zc11/det+jr(i,3)*zc21/det+
  446. & jr(i,4)*zc31/det
  447. jl(i,3)=jl(i,3)+jr(i,2)*zc12/det+jr(i,3)*zc22/det+
  448. & jr(i,4)*zc32/det
  449. jl(i,4)=jl(i,4)+jr(i,2)*zc13/det+jr(i,3)*zc23/det+
  450. & jr(i,4)*zc33/det
  451. jl(i,5)=jl(i,5)+jr(i,5)
  452. 11 continue
  453. c----------------------------------------------------------------------
  454. return
  455. end
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  

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