1 : C CONJWL SOURCE CHAT 05/01/12 22:18:02 5004 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=dsqrt(ga*pold_l/rold_l) 110 : arigh=dsqrt(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------------------------------------- 131 : c Mplus and Mmin are calligraphic 132 : c lettes M+ and M- from the paper 133 : c------------------------------------- 134 : if(dabs(ml) .ge. 1.0d0) then 135 : Mplus=(ml+dabs(ml))/2.0d0 136 : else 137 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0 138 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0) 139 : endif 140 : c----------------- 141 : if(dabs(mr) .ge. 1.0d0) then 142 : Mmin=(mr-dabs(mr))/2.0d0 143 : else 144 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0 145 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0) 146 : endif 147 : c--------------------------------------- 148 : c Derivatives of ml and mr with respect 149 : c to both sets of primitive variables: 150 : c left and right 151 : c-------------------------------------- 152 : temp_l=-un_l/(am*am) 153 : temp_r=-un_r/(am*am) 154 : c-------- 155 : dmlr_l=temp_l*damr_l 156 : dmlr_r=temp_l*damr_r 157 : dmrr_l=temp_r*damr_l 158 : dmrr_r=temp_r*damr_r 159 : c-------- 160 : dmlu_l=n_x/am+temp_l*damu_l 161 : dmlu_r=temp_l*damu_r 162 : dmru_l=temp_r*damu_l 163 : dmru_r=n_x/am+temp_r*damu_r 164 : c-------- 165 : dmlv_l=n_y/am+temp_l*damv_l 166 : dmlv_r=temp_l*damv_r 167 : dmrv_l=temp_r*damv_l 168 : dmrv_r=n_y/am+temp_r*damv_r 169 : c-------- 170 : dmlp_l=temp_l*damp_l 171 : dmlp_r=temp_l*damp_r 172 : dmrp_l=temp_r*damp_l 173 : dmrp_r=temp_r*damp_r 174 : c--------------------------------------------------------------- 174 : c Computing the calligraphic P+ and P- with their derivatives 175 : c--------------------------------------------------------------- 176 : if(ml .ge. 1.0d0) then 177 : Pplus = 1.0d0 178 : else 179 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 180 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0 181 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0) 182 : else 183 : Pplus = 0.0d0 184 : endif 185 : endif 186 : c--------------------------------------------------------------- 187 : if(mr .ge. 1.0d0) then 188 : Pmin = 0.0d0 189 : else 190 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 191 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0 192 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0) 193 : else 194 : Pmin = 1.0d0 195 : endif 196 : endif 197 : c--------------------------------------------------------------- 198 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0 199 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0) 200 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0) 201 : c-------------- 202 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0 203 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0) 204 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0) 205 : c--------------------------------------------------------------- 206 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 207 : dPpr_l=brac_l*dmlr_l 208 : dPpr_r=brac_l*dmlr_r 209 : c------------ 210 : dPpu_l=brac_l*dmlu_l 211 : dPpu_r=brac_l*dmlu_r 212 : c------------ 213 : dPpv_l=brac_l*dmlv_l 214 : dPpv_r=brac_l*dmlv_r 215 : c------------ 216 : dPpp_l=brac_l*dmlp_l 217 : dPpp_r=brac_l*dmlp_r 218 : c------------ 219 : else 220 : dPpr_l=0.0d0 221 : dPpr_r=0.0d0 222 : c----------- 223 : dPpu_l=0.0d0 224 : dPpu_r=0.0d0 225 : c----------- 226 : dPpv_l=0.0d0 227 : dPpv_r=0.0d0 228 : c----------- 229 : dPpp_l=0.0d0 230 : dPpp_r=0.0d0 231 : c----------- 232 : endif 233 : c--------------------------------------------------------------- 234 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 235 : dPmr_l=brac_r*dmrr_l 236 : dPmr_r=brac_r*dmrr_r 237 : c------------ 238 : dPmu_l=brac_r*dmru_l 239 : dPmu_r=brac_r*dmru_r 240 : c------------ 241 : dPmv_l=brac_r*dmrv_l 242 : dPmv_r=brac_r*dmrv_r 243 : c------------ 244 : dPmp_l=brac_r*dmrp_l 245 : dPmp_r=brac_r*dmrp_r 246 : c------------ 247 : else 248 : dPmr_l=0.0d0 249 : dPmr_r=0.0d0 250 : c----------- 251 : dPmu_l=0.0d0 252 : dPmu_r=0.0d0 253 : c----------- 254 : dPmv_l=0.0d0 255 : dPmv_r=0.0d0 256 : c----------- 257 : dPmp_l=0.0d0 258 : dPmp_r=0.0d0 259 : c----------- 260 : endif 261 : c--------------------------------------------------------------------- 262 : c computing pmid -- p_{1/2} and its derivatives 263 : c--------------------------------------------------------------------- 264 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r 265 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r 266 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r 267 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r 268 : c---------------------------- 269 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r 270 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r 271 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r 272 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r 273 : c--------------------------------------------------------------------- 274 : c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!! 275 : c--------------------------------------------------------------------- 276 : jl(1,1)=0.0D0 277 : jl(1,2)=0.0D0 278 : jl(1,3)=0.0D0 279 : jl(1,4)=0.0D0 280 : c------------------------------------ 281 : jr(1,1)=0.0D0 282 : jr(1,2)=0.0D0 283 : jr(1,3)=0.0D0 284 : jr(1,4)=0.0D0 285 : c------------------------------------ 286 : ccccc f222222222222 ------------- 287 : c------------------------------------ 288 : c--------------------------------------------------------- 289 : jl(2,1)=n_x*dpir_l 290 : c------------------- 291 : jl(2,2)=n_x*dpiu_l 292 : c------------------- 293 : jl(2,3)=n_x*dpiv_l 294 : c------------------- 295 : jl(2,4)=n_x*dpip_l 296 : c------------------------------------------------------------- 297 : jr(2,1)=n_x*dpir_r 298 : c------------------- 299 : jr(2,2)=n_x*dpiu_r 300 : c------------------- 301 : jr(2,3)=n_x*dpiv_r 302 : c------------------- 303 : jr(2,4)=n_x*dpip_r 304 : c------------------------------------------------------------- 305 : c------------ f33333333333333333333 --------------------- 306 : c------------------------------------------------------------- 307 : jl(3,1)=n_y*dpir_l 308 : c------------------- 309 : jl(3,2)=n_y*dpiu_l 310 : c------------------- 311 : jl(3,3)=n_y*dpiv_l 312 : c------------------- 313 : jl(3,4)=n_y*dpip_l 314 : c------------------------------------------------------------- 315 : jr(3,1)=n_y*dpir_r 316 : c------------------- 317 : jr(3,2)=n_y*dpiu_r 318 : c------------------- 319 : jr(3,3)=n_y*dpiv_r 320 : c------------------- 321 : jr(3,4)=n_y*dpip_r 322 : c------------------------------------------------------------- 323 : c ------ f44444444444444444444444444444 --------- 324 : c------------------------------------------------------------- 325 : jl(4,1)=0.0D0 326 : c--------------------- 327 : jl(4,2)=0.0D0 328 : c--------------------- 329 : jl(4,3)=0.0D0 330 : c--------------------- 331 : jl(4,4)=0.0D0 332 : c---------------------------------------------------------- 333 : c---------------------------------------------------------- 334 : jr(4,1)=0.0D0 335 : c--------------------- 336 : jr(4,2)=0.0D0 337 : c--------------------- 338 : jr(4,3)=0.0D0 339 : c--------------------- 340 : jr(4,4)=0.0D0 341 : c---------------------------------------------------------- 342 : c matrixes wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww 343 : c---------------------------------------------------------- 344 : wl(1,1)=1.0d0 345 : wl(1,2)=0.0d0 346 : wl(1,3)=0.0d0 347 : wl(1,4)=0.0d0 348 : c------------------------------ 349 : wl(2,1)=-uold_l/rold_l 350 : wl(2,2)=1.0d0/rold_l 351 : wl(2,3)=0.0d0 352 : wl(2,4)=0.0d0 353 : c------------------------------ 354 : wl(3,1)=-vold_l/rold_l 355 : wl(3,2)=0.0d0 356 : wl(3,3)=1.0d0/rold_l 357 : wl(3,4)=0.0d0 358 : c------------------------------ 359 : wl(4,1)=gm1*(uold_l*uold_l+vold_l*vold_l)/2.0d0 360 : wl(4,2)=-uold_l*gm1 361 : wl(4,3)=-vold_l*gm1 362 : wl(4,4)=gm1 363 : c------------------------------ 364 : c------------------------------ 365 : wr(1,1)=1.0d0 366 : wr(1,2)=0.0d0 367 : wr(1,3)=0.0d0 368 : wr(1,4)=0.0d0 369 : c------------------------------ 370 : wr(2,1)=-uold_r/rold_r 371 : wr(2,2)=1.0d0/rold_r 372 : wr(2,3)=0.0d0 373 : wr(2,4)=0.0d0 374 : c------------------------------ 375 : wr(3,1)=-vold_r/rold_r 376 : wr(3,2)=0.0d0 377 : wr(3,3)=1.0d0/rold_r 378 : wr(3,4)=0.0d0 379 : c------------------------------ 380 : wr(4,1)=gm1*(uold_r*uold_r+vold_r*vold_r)/2.0d0 381 : wr(4,2)=-uold_r*gm1 382 : wr(4,3)=-vold_r*gm1 383 : wr(4,4)=gm1 384 : c---------------------------------------------- 385 : c---------------------------------------------- 386 : do 1 i=1,4 387 : do 2 j=1,4 388 : jtl(i,j)=0.0d0 389 : jtr(i,j)=0.0d0 390 : do 3 k=1,4 391 : jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j) 392 : jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j) 393 : 3 continue 394 : 2 continue 395 : 1 continue 396 : c----------------------------- 397 : tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2) 398 : bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2) 399 : c----------------------------- 400 : do 11 i=1,4 401 : jtl(i,1)=jtl(i,1)+jtr(i,1) 402 : jtl(i,2)=jtl(i,2)+jtr(i,2)*(-tcoef/bcoef)+ 403 : & jtr(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef 404 : jtl(i,3)=jtl(i,3)+jtr(i,2)*(-2.0d0*nvect(2)*tvect(2)/bcoef)+ 405 : & jtr(i,3)*tcoef/bcoef 406 : jtl(i,4)=jtl(i,4)+jtr(i,4) 407 : 11 continue 408 : c---------------------------------------------------------------------- 409 : return 410 : end 411 : 412 : 413 : 414 : 415 : 416 : 417 : 418 : 419 : 420 : 421 : 422 : 423 : 424 : 425 : 426 : 427 : 428 : 429 :
© Cast3M 2003 - All rights reserved.
Disclaimer