1 : C CONJP3 SOURCE CHAT 05/01/12 22:17:40 5004 2 : SUBROUTINE CONJP3(JL,JR,WVEC_L,WVEC_R,NVECT,TVECT, 3 : & ga) 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : CONJP3 9 : C 10 : C DESCRIPTION : Voir KONJP2 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 interface 20 : c next to wall. 21 : c For full descriptions please see file 'conjp2.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 jl -- jakobian matrix 4 by 4 - derivatives of the numerical 47 : c flux function with respect to the primitive variables 48 : c from the left cell; 49 : c 50 : c jr -- 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 alpha,beta 59 : real*8 ga, gm1 60 : real*8 n_x,n_y 61 : real*8 un_l, un_r 62 : real*8 ml,mr,Mplus,Mmin 63 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l 64 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r 65 : real*8 Pplus,Pmin 66 : real*8 temp_l,temp_r,brac_l,brac_r 67 : real*8 aleft, arigh, am 68 : real*8 damr_l,damr_r,damu_l,damu_r 69 : real*8 damv_l,damv_r,damp_l,damp_r 70 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r 71 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_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 dPpr_l,dPpr_r,dPpu_l,dPpu_r 75 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r 76 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r 77 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r 78 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r 79 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r 80 : real*8 tcoef,bcoef 81 : integer i 82 : parameter(alpha = 0.1875D0, beta = 0.125D0) 83 : c----------------------- 84 : gm1=ga-1.0d0 85 : c----------------------- 86 : n_x=nvect(1) 87 : n_y=nvect(2) 88 : c---------------------------- 89 : c----------------------- 90 : rold_l=wvec_l(1) 91 : uold_l=wvec_l(2) 92 : vold_l=wvec_l(3) 93 : pold_l=wvec_l(4) 94 : c----------------------- 95 : rold_r=wvec_r(1) 96 : uold_r=wvec_r(2) 97 : vold_r=wvec_r(3) 98 : pold_r=wvec_r(4) 99 : c----------------------- 100 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0 101 : eold_l=eold_l+pold_l/(gm1*rold_l) 102 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0 103 : eold_r=eold_r+pold_r/(gm1*rold_r) 104 : c------------------------------------------------------------------- 105 : c Computation of the speed of sound and its derivatives 106 : c--------------------------------------------------------------------- 107 : aleft=dsqrt(ga*pold_l/rold_l) 108 : arigh=dsqrt(ga*pold_r/rold_r) 109 : am=0.5d0*(aleft+arigh) 110 : c-------------------------------------------------------------------- 111 : damr_r=-arigh/(4.0d0*rold_r) 112 : damu_r=0.0d0 113 : damv_r=0.0d0 114 : damp_r=ga/(4.0d0*arigh*rold_r) 115 : c----------------------- 116 : damr_l=-aleft/(4.0d0*rold_l) 117 : damu_l=0.0d0 118 : damv_l=0.0d0 119 : damp_l=ga/(4.0d0*aleft*rold_l) 120 : c---------------------------------------------------------------------- 121 : c computing numerical Mach number and its derivatives 122 : c---------------------------------------------------------------------- 123 : un_l=uold_l*n_x+vold_l*n_y 124 : un_r=uold_r*n_x+vold_r*n_y 125 : c---------------------------------------- 126 : ml=un_l/am 127 : mr=un_r/am 128 : c------------------------------------- 128 : c Mplus and Mmin are calligraphic 129 : c lettes M+ and M- from the paper 130 : c------------------------------------- 131 : if(dabs(ml) .ge. 1.0d0) then 132 : Mplus=(ml+dabs(ml))/2.0d0 133 : else 134 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0 135 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0) 136 : endif 137 : c----------------- 138 : if(dabs(mr) .ge. 1.0d0) then 139 : Mmin=(mr-dabs(mr))/2.0d0 140 : else 141 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0 142 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0) 143 : endif 144 : c--------------------------------------- 145 : c Derivatives of ml and mr with respect 146 : c to both sets of primitive variables: 147 : c left and right 148 : c-------------------------------------- 149 : temp_l=-un_l/(am*am) 150 : temp_r=-un_r/(am*am) 151 : c-------- 152 : dmlr_l=temp_l*damr_l 153 : dmlr_r=temp_l*damr_r 154 : dmrr_l=temp_r*damr_l 155 : dmrr_r=temp_r*damr_r 156 : c-------- 157 : dmlu_l=n_x/am+temp_l*damu_l 158 : dmlu_r=temp_l*damu_r 159 : dmru_l=temp_r*damu_l 160 : dmru_r=n_x/am+temp_r*damu_r 161 : c-------- 162 : dmlv_l=n_y/am+temp_l*damv_l 163 : dmlv_r=temp_l*damv_r 164 : dmrv_l=temp_r*damv_l 165 : dmrv_r=n_y/am+temp_r*damv_r 166 : c-------- 167 : dmlp_l=temp_l*damp_l 168 : dmlp_r=temp_l*damp_r 169 : dmrp_l=temp_r*damp_l 170 : dmrp_r=temp_r*damp_r 171 : c--------------------------------------------------------------- 171 : c Computing the calligraphic P+ and P- with their derivatives 172 : c--------------------------------------------------------------- 173 : if(ml .ge. 1.0d0) then 174 : Pplus = 1.0d0 175 : else 176 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 177 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0 178 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0) 179 : else 180 : Pplus = 0.0d0 181 : endif 182 : endif 183 : c--------------------------------------------------------------- 184 : if(mr .ge. 1.0d0) then 185 : Pmin = 0.0d0 186 : else 187 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 188 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0 189 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0) 190 : else 191 : Pmin = 1.0d0 192 : endif 193 : endif 194 : c--------------------------------------------------------------- 195 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0 196 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0) 197 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0) 198 : c-------------- 199 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0 200 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0) 201 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0) 202 : c--------------------------------------------------------------- 203 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 204 : dPpr_l=brac_l*dmlr_l 205 : dPpr_r=brac_l*dmlr_r 206 : c------------ 207 : dPpu_l=brac_l*dmlu_l 208 : dPpu_r=brac_l*dmlu_r 209 : c------------ 210 : dPpv_l=brac_l*dmlv_l 211 : dPpv_r=brac_l*dmlv_r 212 : c------------ 213 : dPpp_l=brac_l*dmlp_l 214 : dPpp_r=brac_l*dmlp_r 215 : c------------ 216 : else 217 : dPpr_l=0.0d0 218 : dPpr_r=0.0d0 219 : c----------- 220 : dPpu_l=0.0d0 221 : dPpu_r=0.0d0 222 : c----------- 223 : dPpv_l=0.0d0 224 : dPpv_r=0.0d0 225 : c----------- 226 : dPpp_l=0.0d0 227 : dPpp_r=0.0d0 228 : c----------- 229 : endif 230 : c--------------------------------------------------------------- 231 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 232 : dPmr_l=brac_r*dmrr_l 233 : dPmr_r=brac_r*dmrr_r 234 : c------------ 235 : dPmu_l=brac_r*dmru_l 236 : dPmu_r=brac_r*dmru_r 237 : c------------ 238 : dPmv_l=brac_r*dmrv_l 239 : dPmv_r=brac_r*dmrv_r 240 : c------------ 241 : dPmp_l=brac_r*dmrp_l 242 : dPmp_r=brac_r*dmrp_r 243 : c------------ 244 : else 245 : dPmr_l=0.0d0 246 : dPmr_r=0.0d0 247 : c----------- 248 : dPmu_l=0.0d0 249 : dPmu_r=0.0d0 250 : c----------- 251 : dPmv_l=0.0d0 252 : dPmv_r=0.0d0 253 : c----------- 254 : dPmp_l=0.0d0 255 : dPmp_r=0.0d0 256 : c----------- 257 : endif 258 : c--------------------------------------------------------------------- 259 : c computing pmid -- p_{1/2} and its derivatives 260 : c--------------------------------------------------------------------- 261 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r 262 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r 263 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r 264 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r 265 : c---------------------------- 266 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r 267 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r 268 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r 269 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r 270 : c--------------------------------------------------------------------- 271 : c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!! 272 : c--------------------------------------------------------------------- 273 : jl(1,1)=0.0D0 274 : jl(1,2)=0.0D0 275 : jl(1,3)=0.0D0 276 : jl(1,4)=0.0D0 277 : c------------------------------------ 278 : jr(1,1)=0.0D0 279 : jr(1,2)=0.0D0 280 : jr(1,3)=0.0D0 281 : jr(1,4)=0.0D0 282 : c------------------------------------ 283 : ccccc f222222222222 ------------- 284 : c------------------------------------ 285 : c--------------------------------------------------------- 286 : jl(2,1)=n_x*dpir_l 287 : c------------------- 288 : jl(2,2)=n_x*dpiu_l 289 : c------------------- 290 : jl(2,3)=n_x*dpiv_l 291 : c------------------- 292 : jl(2,4)=n_x*dpip_l 293 : c------------------------------------------------------------- 294 : jr(2,1)=n_x*dpir_r 295 : c------------------- 296 : jr(2,2)=n_x*dpiu_r 297 : c------------------- 298 : jr(2,3)=n_x*dpiv_r 299 : c------------------- 300 : jr(2,4)=n_x*dpip_r 301 : c------------------------------------------------------------- 302 : c------------ f33333333333333333333 --------------------- 303 : c------------------------------------------------------------- 304 : jl(3,1)=n_y*dpir_l 305 : c------------------- 306 : jl(3,2)=n_y*dpiu_l 307 : c------------------- 308 : jl(3,3)=n_y*dpiv_l 309 : c------------------- 310 : jl(3,4)=n_y*dpip_l 311 : c------------------------------------------------------------- 312 : jr(3,1)=n_y*dpir_r 313 : c------------------- 314 : jr(3,2)=n_y*dpiu_r 315 : c------------------- 316 : jr(3,3)=n_y*dpiv_r 317 : c------------------- 318 : jr(3,4)=n_y*dpip_r 319 : c------------------------------------------------------------- 320 : c ------ f44444444444444444444444444444 --------- 321 : c------------------------------------------------------------- 322 : jl(4,1)=0.0D0 323 : c--------------------- 324 : jl(4,2)=0.0D0 325 : c--------------------- 326 : jl(4,3)=0.0D0 327 : c--------------------- 328 : jl(4,4)=0.0D0 329 : c---------------------------------------------------------- 330 : c---------------------------------------------------------- 331 : jr(4,1)=0.0D0 332 : c--------------------- 333 : jr(4,2)=0.0D0 334 : c--------------------- 335 : jr(4,3)=0.0D0 336 : c--------------------- 337 : jr(4,4)=0.0D0 338 : c---------------------------------------------------------- 339 : c----------------------------- 340 : tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2) 341 : bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2) 342 : c----------------------------- 343 : do 11 i=1,4 344 : jl(i,1)=jl(i,1)+jr(i,1) 345 : jl(i,2)=jl(i,2)+jr(i,2)*(-tcoef/bcoef)+ 346 : & jr(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef 347 : jl(i,3)=jl(i,3)+jr(i,2)*(-2.0d0*nvect(2)*tvect(2)/bcoef)+ 348 : & jr(i,3)*tcoef/bcoef 349 : jl(i,4)=jl(i,4)+jr(i,4) 350 : 11 continue 351 : c---------------------------------------------------------------------- 352 : return 353 : end 354 : 355 : 356 : 357 : 358 : 359 : 360 : 361 : 362 : 363 : 364 : 365 : 366 : 367 : 368 : 369 : 370 : 371 : 372 : 373 :
© Cast3M 2003 - All rights reserved.
Disclaimer