Télécharger conjwl.eso

Retour à la liste

Numérotation des lignes :

  1. C CONJWL SOURCE CB215821 16/04/21 21:16:00 8920
  2. SUBROUTINE CONJWL(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVECT,
  3. & ga)
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : CONJWL
  9. C
  10. C DESCRIPTION : Voir KONJA2
  11. C
  12. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. c----------------------------------------------------------------------
  18. c GENERAL DESCRIPTION:
  19. c This subroutine provides jacobian jtl at the interfface
  20. c next to wall.
  21. c For full descriptions please see file 'conjak.eso'.
  22. c----------------------------------------------------------------------
  23. c INPUT:
  24. c
  25. c alpha -- parameter of the AUSM+ scheme in the Pressure function;
  26. c ( -3/4 <= alpha <= 3/16 )
  27. c
  28. c beta -- parameter of the AUSM+ scheme in the Mach function;
  29. c ( -1/16 <= beta <= 1/2 )
  30. c
  31. c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the
  32. c left cell;
  33. c
  34. c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the
  35. c right cell;
  36. c
  37. c nvect -- normal vector to the interface (2 components in 2D);
  38. c
  39. c tvect -- tangential vector to the interface;
  40. c
  41. c ga -- ratio of the specific heats (assumed constant)
  42. c----------------------------------------------------------------------
  43. c
  44. c OUTPUT:
  45. c
  46. c jtl -- jakobian matrix 4 by 4 - derivatives of the numerical
  47. c flux function with respect to the conservative variables
  48. c from the left cell;
  49. c
  50. c jtr -- jakobian matrix 4 by 4 - derivatives of the numerical
  51. c flux function with respect to the conservative variables
  52. c from the right cell.
  53. c----------------------------------------------------------------------
  54. IMPLICIT INTEGER(I-N)
  55. real*8 wvec_l(4),wvec_r(4)
  56. real*8 nvect(2),tvect(2)
  57. real*8 jl(4,4),jr(4,4)
  58. real*8 wl(4,4),wr(4,4)
  59. real*8 jtl(4,4),jtr(4,4)
  60. real*8 alpha,beta
  61. real*8 ga, gm1
  62. real*8 n_x,n_y
  63. real*8 un_l, un_r
  64. real*8 ml,mr,Mplus,Mmin
  65. real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  66. real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  67. real*8 Pplus,Pmin
  68. real*8 temp_l,temp_r,brac_l,brac_r
  69. real*8 aleft, arigh, am
  70. real*8 damr_l,damr_r,damu_l,damu_r
  71. real*8 damv_l,damv_r,damp_l,damp_r
  72. real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  73. real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  74. real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  75. real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  76. real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  77. real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_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 dpir_l,dpir_r,dpiu_l,dpiu_r
  81. real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  82. real*8 tcoef,bcoef
  83. integer i,j,k
  84. parameter(alpha = 0.1875D0, beta = 0.125D0)
  85. c-----------------------
  86. gm1=ga-1.0d0
  87. c-----------------------
  88. n_x=nvect(1)
  89. n_y=nvect(2)
  90. c----------------------------
  91. c-----------------------
  92. rold_l=wvec_l(1)
  93. uold_l=wvec_l(2)
  94. vold_l=wvec_l(3)
  95. pold_l=wvec_l(4)
  96. c-----------------------
  97. rold_r=wvec_r(1)
  98. uold_r=wvec_r(2)
  99. vold_r=wvec_r(3)
  100. pold_r=wvec_r(4)
  101. c-----------------------
  102. eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  103. eold_l=eold_l+pold_l/(gm1*rold_l)
  104. eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  105. eold_r=eold_r+pold_r/(gm1*rold_r)
  106. c-------------------------------------------------------------------
  107. c Computation of the speed of sound and its derivatives
  108. c---------------------------------------------------------------------
  109. aleft=SQRT(ga*pold_l/rold_l)
  110. arigh=SQRT(ga*pold_r/rold_r)
  111. am=0.5d0*(aleft+arigh)
  112. c--------------------------------------------------------------------
  113. damr_r=-arigh/(4.0d0*rold_r)
  114. damu_r=0.0d0
  115. damv_r=0.0d0
  116. damp_r=ga/(4.0d0*arigh*rold_r)
  117. c-----------------------
  118. damr_l=-aleft/(4.0d0*rold_l)
  119. damu_l=0.0d0
  120. damv_l=0.0d0
  121. damp_l=ga/(4.0d0*aleft*rold_l)
  122. c----------------------------------------------------------------------
  123. c computing numerical Mach number and its derivatives
  124. c----------------------------------------------------------------------
  125. un_l=uold_l*n_x+vold_l*n_y
  126. un_r=uold_r*n_x+vold_r*n_y
  127. c----------------------------------------
  128. ml=un_l/am
  129. mr=un_r/am
  130. c print*,'ml,mr',ml,mr
  131. c-------------------------------------
  132. c Mplus and Mmin are calligraphic
  133. c lettes M+ and M- from the paper
  134. c-------------------------------------
  135. if(ABS(ml) .ge. 1.0d0) then
  136. Mplus=(ml+ABS(ml))/2.0d0
  137. else
  138. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  139. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  140. endif
  141. c-----------------
  142. if(ABS(mr) .ge. 1.0d0) then
  143. Mmin=(mr-ABS(mr))/2.0d0
  144. else
  145. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  146. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  147. endif
  148. c---------------------------------------
  149. c Derivatives of ml and mr with respect
  150. c to both sets of primitive variables:
  151. c left and right
  152. c--------------------------------------
  153. temp_l=-un_l/(am*am)
  154. temp_r=-un_r/(am*am)
  155. c--------
  156. dmlr_l=temp_l*damr_l
  157. dmlr_r=temp_l*damr_r
  158. dmrr_l=temp_r*damr_l
  159. dmrr_r=temp_r*damr_r
  160. c--------
  161. dmlu_l=n_x/am+temp_l*damu_l
  162. dmlu_r=temp_l*damu_r
  163. dmru_l=temp_r*damu_l
  164. dmru_r=n_x/am+temp_r*damu_r
  165. c--------
  166. dmlv_l=n_y/am+temp_l*damv_l
  167. dmlv_r=temp_l*damv_r
  168. dmrv_l=temp_r*damv_l
  169. dmrv_r=n_y/am+temp_r*damv_r
  170. c--------
  171. dmlp_l=temp_l*damp_l
  172. dmlp_r=temp_l*damp_r
  173. dmrp_l=temp_r*damp_l
  174. dmrp_r=temp_r*damp_r
  175. c---------------------------------------------------------------
  176. c Computing the calligraphic P+ and P- with their derivatives
  177. c---------------------------------------------------------------
  178. if(ml .ge. 1.0d0) then
  179. Pplus = 1.0d0
  180. else
  181. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  182. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  183. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  184. else
  185. Pplus = 0.0d0
  186. endif
  187. endif
  188. c---------------------------------------------------------------
  189. if(mr .ge. 1.0d0) then
  190. Pmin = 0.0d0
  191. else
  192. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  193. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  194. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  195. else
  196. Pmin = 1.0d0
  197. endif
  198. endif
  199. c---------------------------------------------------------------
  200. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  201. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  202. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  203. c--------------
  204. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  205. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  206. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  207. c---------------------------------------------------------------
  208. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  209. dPpr_l=brac_l*dmlr_l
  210. dPpr_r=brac_l*dmlr_r
  211. c------------
  212. dPpu_l=brac_l*dmlu_l
  213. dPpu_r=brac_l*dmlu_r
  214. c------------
  215. dPpv_l=brac_l*dmlv_l
  216. dPpv_r=brac_l*dmlv_r
  217. c------------
  218. dPpp_l=brac_l*dmlp_l
  219. dPpp_r=brac_l*dmlp_r
  220. c------------
  221. else
  222. dPpr_l=0.0d0
  223. dPpr_r=0.0d0
  224. c-----------
  225. dPpu_l=0.0d0
  226. dPpu_r=0.0d0
  227. c-----------
  228. dPpv_l=0.0d0
  229. dPpv_r=0.0d0
  230. c-----------
  231. dPpp_l=0.0d0
  232. dPpp_r=0.0d0
  233. c-----------
  234. endif
  235. c---------------------------------------------------------------
  236. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  237. dPmr_l=brac_r*dmrr_l
  238. dPmr_r=brac_r*dmrr_r
  239. c------------
  240. dPmu_l=brac_r*dmru_l
  241. dPmu_r=brac_r*dmru_r
  242. c------------
  243. dPmv_l=brac_r*dmrv_l
  244. dPmv_r=brac_r*dmrv_r
  245. c------------
  246. dPmp_l=brac_r*dmrp_l
  247. dPmp_r=brac_r*dmrp_r
  248. c------------
  249. else
  250. dPmr_l=0.0d0
  251. dPmr_r=0.0d0
  252. c-----------
  253. dPmu_l=0.0d0
  254. dPmu_r=0.0d0
  255. c-----------
  256. dPmv_l=0.0d0
  257. dPmv_r=0.0d0
  258. c-----------
  259. dPmp_l=0.0d0
  260. dPmp_r=0.0d0
  261. c-----------
  262. endif
  263. c---------------------------------------------------------------------
  264. c computing pmid -- p_{1/2} and its derivatives
  265. c---------------------------------------------------------------------
  266. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  267. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  268. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  269. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  270. c----------------------------
  271. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  272. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  273. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  274. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  275. c---------------------------------------------------------------------
  276. c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  277. c---------------------------------------------------------------------
  278. jl(1,1)=0.0D0
  279. jl(1,2)=0.0D0
  280. jl(1,3)=0.0D0
  281. jl(1,4)=0.0D0
  282. c------------------------------------
  283. jr(1,1)=0.0D0
  284. jr(1,2)=0.0D0
  285. jr(1,3)=0.0D0
  286. jr(1,4)=0.0D0
  287. c------------------------------------
  288. ccccc f222222222222 -------------
  289. c------------------------------------
  290. c---------------------------------------------------------
  291. jl(2,1)=n_x*dpir_l
  292. c-------------------
  293. jl(2,2)=n_x*dpiu_l
  294. c-------------------
  295. jl(2,3)=n_x*dpiv_l
  296. c-------------------
  297. jl(2,4)=n_x*dpip_l
  298. c-------------------------------------------------------------
  299. jr(2,1)=n_x*dpir_r
  300. c-------------------
  301. jr(2,2)=n_x*dpiu_r
  302. c-------------------
  303. jr(2,3)=n_x*dpiv_r
  304. c-------------------
  305. jr(2,4)=n_x*dpip_r
  306. c-------------------------------------------------------------
  307. c------------ f33333333333333333333 ---------------------
  308. c-------------------------------------------------------------
  309. jl(3,1)=n_y*dpir_l
  310. c-------------------
  311. jl(3,2)=n_y*dpiu_l
  312. c-------------------
  313. jl(3,3)=n_y*dpiv_l
  314. c-------------------
  315. jl(3,4)=n_y*dpip_l
  316. c-------------------------------------------------------------
  317. jr(3,1)=n_y*dpir_r
  318. c-------------------
  319. jr(3,2)=n_y*dpiu_r
  320. c-------------------
  321. jr(3,3)=n_y*dpiv_r
  322. c-------------------
  323. jr(3,4)=n_y*dpip_r
  324. c-------------------------------------------------------------
  325. c ------ f44444444444444444444444444444 ---------
  326. c-------------------------------------------------------------
  327. jl(4,1)=0.0D0
  328. c---------------------
  329. jl(4,2)=0.0D0
  330. c---------------------
  331. jl(4,3)=0.0D0
  332. c---------------------
  333. jl(4,4)=0.0D0
  334. c----------------------------------------------------------
  335. c----------------------------------------------------------
  336. jr(4,1)=0.0D0
  337. c---------------------
  338. jr(4,2)=0.0D0
  339. c---------------------
  340. jr(4,3)=0.0D0
  341. c---------------------
  342. jr(4,4)=0.0D0
  343. c----------------------------------------------------------
  344. c matrixes wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
  345. c----------------------------------------------------------
  346. wl(1,1)=1.0d0
  347. wl(1,2)=0.0d0
  348. wl(1,3)=0.0d0
  349. wl(1,4)=0.0d0
  350. c------------------------------
  351. wl(2,1)=-uold_l/rold_l
  352. wl(2,2)=1.0d0/rold_l
  353. wl(2,3)=0.0d0
  354. wl(2,4)=0.0d0
  355. c------------------------------
  356. wl(3,1)=-vold_l/rold_l
  357. wl(3,2)=0.0d0
  358. wl(3,3)=1.0d0/rold_l
  359. wl(3,4)=0.0d0
  360. c------------------------------
  361. wl(4,1)=gm1*(uold_l*uold_l+vold_l*vold_l)/2.0d0
  362. wl(4,2)=-uold_l*gm1
  363. wl(4,3)=-vold_l*gm1
  364. wl(4,4)=gm1
  365. c------------------------------
  366. c------------------------------
  367. wr(1,1)=1.0d0
  368. wr(1,2)=0.0d0
  369. wr(1,3)=0.0d0
  370. wr(1,4)=0.0d0
  371. c------------------------------
  372. wr(2,1)=-uold_r/rold_r
  373. wr(2,2)=1.0d0/rold_r
  374. wr(2,3)=0.0d0
  375. wr(2,4)=0.0d0
  376. c------------------------------
  377. wr(3,1)=-vold_r/rold_r
  378. wr(3,2)=0.0d0
  379. wr(3,3)=1.0d0/rold_r
  380. wr(3,4)=0.0d0
  381. c------------------------------
  382. wr(4,1)=gm1*(uold_r*uold_r+vold_r*vold_r)/2.0d0
  383. wr(4,2)=-uold_r*gm1
  384. wr(4,3)=-vold_r*gm1
  385. wr(4,4)=gm1
  386. c----------------------------------------------
  387. c----------------------------------------------
  388. do 1 i=1,4
  389. do 2 j=1,4
  390. jtl(i,j)=0.0d0
  391. jtr(i,j)=0.0d0
  392. do 3 k=1,4
  393. jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
  394. jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
  395. 3 continue
  396. 2 continue
  397. 1 continue
  398. c-----------------------------
  399. tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2)
  400. bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2)
  401. c-----------------------------
  402. do 11 i=1,4
  403. jtl(i,1)=jtl(i,1)+jtr(i,1)
  404. jtl(i,2)=jtl(i,2)+jtr(i,2)*(-tcoef/bcoef)+
  405. & jtr(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef
  406. jtl(i,3)=jtl(i,3)+jtr(i,2)*(-2.0d0*nvect(2)*tvect(2)/bcoef)+
  407. & jtr(i,3)*tcoef/bcoef
  408. jtl(i,4)=jtl(i,4)+jtr(i,4)
  409. 11 continue
  410. c----------------------------------------------------------------------
  411. return
  412. end
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  

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