1 : C AUSMW SOURCE CHAT 05/01/12 21:29:23 5004 2 : c---------------------------------------------------------------------- 3 : c INPUT: 4 : c 5 : c rog : left density 6 : C 7 : c uxg : left x-component of speed 8 : c 9 : c uyg : left y-component of speed 10 : c 11 : c pg : left pressure 12 : c 13 : c rod : right density 14 : c 15 : c uxd : right x-component of speed 16 : c 17 : c uyd : right y-component of speed 18 : c 19 : c pd : right pressure 20 : c 21 : c ga : gamma of the gas 22 : c 23 : c cnx : x-component of the surface normal 24 : c 25 : c cny : y-component of the surface normal 26 : c 27 : c OUTPUT: 28 : c 29 : c dfrun : derivative of the numerical estimation of (rho u_n u_n + 30 : c p) (flux-function for conservative variable (rho u_n)) 31 : c with respect to the primitive variables 32 : c 33 : c---------------------------------------------------------------------- 34 : SUBROUTINE AUSMW(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD, 35 : & ga,cnx,cny,dfrun) 36 : IMPLICIT INTEGER(I-N) 37 : real*8 rog,uxg,uyg,pg 38 : real*8 rod,uxd,uyd,pd 39 : real*8 cnx,cny 40 : real*8 jl(4,4),jr(4,4),dfrun(4) 41 : real*8 alpha,beta 42 : real*8 ga, gm1 43 : real*8 n_x,n_y,t_x,t_y 44 : real*8 un_l, un_r 45 : real*8 ml,mr,Mplus,Mmin 46 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l 47 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r 48 : real*8 Pplus,Pmin 49 : real*8 temp_l,temp_r,brac_l,brac_r 50 : real*8 aleft, arigh, am 51 : real*8 damr_l,damr_r,damu_l,damu_r 52 : real*8 damv_l,damv_r,damp_l,damp_r 53 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r 54 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r 55 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r 56 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r 57 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r 58 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r 59 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r 60 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r 61 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r 62 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r 63 : real*8 tcoef,bcoef 64 : integer i 65 : parameter(alpha = 0.1875D0, beta = 0.125D0) 66 : c----------------------- 67 : gm1=ga-1.0d0 68 : c----------------------- 69 : n_x=cnx 70 : n_y=cny 71 : t_x=-cny 72 : t_y=cnx 73 : c---------------------------- 74 : c----------------------- 75 : rold_l=rog 76 : uold_l=uxg 77 : vold_l=uyg 78 : pold_l=pg 79 : c----------------------- 80 : rold_r=rod 81 : uold_r=uxd 82 : vold_r=uyd 83 : pold_r=pd 84 : c----------------------- 85 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0 86 : eold_l=eold_l+pold_l/(gm1*rold_l) 87 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0 88 : eold_r=eold_r+pold_r/(gm1*rold_r) 89 : c------------------------------------------------------------------- 90 : c Computation of the speed of sound and its derivatives 91 : c--------------------------------------------------------------------- 92 : aleft=dsqrt(ga*pold_l/rold_l) 93 : arigh=dsqrt(ga*pold_r/rold_r) 94 : am=0.5d0*(aleft+arigh) 95 : c-------------------------------------------------------------------- 96 : damr_r=-arigh/(4.0d0*rold_r) 97 : damu_r=0.0d0 98 : damv_r=0.0d0 99 : damp_r=ga/(4.0d0*arigh*rold_r) 100 : c----------------------- 101 : damr_l=-aleft/(4.0d0*rold_l) 102 : damu_l=0.0d0 103 : damv_l=0.0d0 104 : damp_l=ga/(4.0d0*aleft*rold_l) 105 : c---------------------------------------------------------------------- 106 : c computing numerical Mach number and its derivatives 107 : c---------------------------------------------------------------------- 108 : un_l=uold_l*n_x+vold_l*n_y 109 : un_r=uold_r*n_x+vold_r*n_y 110 : c---------------------------------------- 111 : ml=un_l/am 112 : mr=un_r/am 113 : c------------------------------------- 113 : c Mplus and Mmin are calligraphic 114 : c lettes M+ and M- from the paper 115 : c------------------------------------- 116 : if(dabs(ml) .ge. 1.0d0) then 117 : Mplus=(ml+dabs(ml))/2.0d0 118 : else 119 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0 120 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0) 121 : endif 122 : c----------------- 123 : if(dabs(mr) .ge. 1.0d0) then 124 : Mmin=(mr-dabs(mr))/2.0d0 125 : else 126 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0 127 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0) 128 : endif 129 : c--------------------------------------- 130 : c Derivatives of ml and mr with respect 131 : c to both sets of primitive variables: 132 : c left and right 133 : c-------------------------------------- 134 : temp_l=-un_l/(am*am) 135 : temp_r=-un_r/(am*am) 136 : c-------- 137 : dmlr_l=temp_l*damr_l 138 : dmlr_r=temp_l*damr_r 139 : dmrr_l=temp_r*damr_l 140 : dmrr_r=temp_r*damr_r 141 : c-------- 142 : dmlu_l=n_x/am+temp_l*damu_l 143 : dmlu_r=temp_l*damu_r 144 : dmru_l=temp_r*damu_l 145 : dmru_r=n_x/am+temp_r*damu_r 146 : c-------- 147 : dmlv_l=n_y/am+temp_l*damv_l 148 : dmlv_r=temp_l*damv_r 149 : dmrv_l=temp_r*damv_l 150 : dmrv_r=n_y/am+temp_r*damv_r 151 : c-------- 152 : dmlp_l=temp_l*damp_l 153 : dmlp_r=temp_l*damp_r 154 : dmrp_l=temp_r*damp_l 155 : dmrp_r=temp_r*damp_r 156 : c--------------------------------------------------------------- 156 : c Computing the calligraphic P+ and P- with their derivatives 157 : c--------------------------------------------------------------- 158 : if(ml .ge. 1.0d0) then 159 : Pplus = 1.0d0 160 : else 161 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 162 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0 163 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0) 164 : else 165 : Pplus = 0.0d0 166 : endif 167 : endif 168 : c--------------------------------------------------------------- 169 : if(mr .ge. 1.0d0) then 170 : Pmin = 0.0d0 171 : else 172 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 173 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0 174 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0) 175 : else 176 : Pmin = 1.0d0 177 : endif 178 : endif 179 : c--------------------------------------------------------------- 180 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0 181 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0) 182 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0) 183 : c-------------- 184 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0 185 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0) 186 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0) 187 : c--------------------------------------------------------------- 188 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 189 : dPpr_l=brac_l*dmlr_l 190 : dPpr_r=brac_l*dmlr_r 191 : c------------ 192 : dPpu_l=brac_l*dmlu_l 193 : dPpu_r=brac_l*dmlu_r 194 : c------------ 195 : dPpv_l=brac_l*dmlv_l 196 : dPpv_r=brac_l*dmlv_r 197 : c------------ 198 : dPpp_l=brac_l*dmlp_l 199 : dPpp_r=brac_l*dmlp_r 200 : c------------ 201 : else 202 : dPpr_l=0.0d0 203 : dPpr_r=0.0d0 204 : c----------- 205 : dPpu_l=0.0d0 206 : dPpu_r=0.0d0 207 : c----------- 208 : dPpv_l=0.0d0 209 : dPpv_r=0.0d0 210 : c----------- 211 : dPpp_l=0.0d0 212 : dPpp_r=0.0d0 213 : c----------- 214 : endif 215 : c--------------------------------------------------------------- 216 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 217 : dPmr_l=brac_r*dmrr_l 218 : dPmr_r=brac_r*dmrr_r 219 : c------------ 220 : dPmu_l=brac_r*dmru_l 221 : dPmu_r=brac_r*dmru_r 222 : c------------ 223 : dPmv_l=brac_r*dmrv_l 224 : dPmv_r=brac_r*dmrv_r 225 : c------------ 226 : dPmp_l=brac_r*dmrp_l 227 : dPmp_r=brac_r*dmrp_r 228 : c------------ 229 : else 230 : dPmr_l=0.0d0 231 : dPmr_r=0.0d0 232 : c----------- 233 : dPmu_l=0.0d0 234 : dPmu_r=0.0d0 235 : c----------- 236 : dPmv_l=0.0d0 237 : dPmv_r=0.0d0 238 : c----------- 239 : dPmp_l=0.0d0 240 : dPmp_r=0.0d0 241 : c----------- 242 : endif 243 : c--------------------------------------------------------------------- 244 : c computing pmid -- p_{1/2} and its derivatives 245 : c--------------------------------------------------------------------- 246 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r 247 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r 248 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r 249 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r 250 : c---------------------------- 251 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r 252 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r 253 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r 254 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r 255 : c--------------------------------------------------------------------- 256 : c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!! 257 : c--------------------------------------------------------------------- 258 : jl(1,1)=0.0D0 259 : jl(1,2)=0.0D0 260 : jl(1,3)=0.0D0 261 : jl(1,4)=0.0D0 262 : c------------------------------------ 263 : jr(1,1)=0.0D0 264 : jr(1,2)=0.0D0 265 : jr(1,3)=0.0D0 266 : jr(1,4)=0.0D0 267 : c------------------------------------ 268 : ccccc f222222222222 ------------- 269 : c------------------------------------ 270 : c--------------------------------------------------------- 271 : c-------------------------------------------- 272 : jl(2,1)=n_x*dpir_l 273 : c------------------- 274 : jl(2,2)=n_x*dpiu_l 275 : c------------------- 276 : jl(2,3)=n_x*dpiv_l 277 : c------------------- 278 : jl(2,4)=n_x*dpip_l 279 : c------------------------------------------------------------- 280 : jr(2,1)=n_x*dpir_r 281 : c------------------- 282 : jr(2,2)=n_x*dpiu_r 283 : c------------------- 284 : jr(2,3)=n_x*dpiv_r 285 : c------------------- 286 : jr(2,4)=n_x*dpip_r 287 : c------------------------------------------------------------- 288 : c------------ f33333333333333333333 --------------------- 289 : c------------------------------------------------------------- 290 : jl(3,1)=n_y*dpir_l 291 : c------------------- 292 : jl(3,2)=n_y*dpiu_l 293 : c------------------- 294 : jl(3,3)=n_y*dpiv_l 295 : c------------------- 296 : jl(3,4)=n_y*dpip_l 297 : c------------------------------------------------------------- 298 : jr(3,1)=n_y*dpir_r 299 : c------------------- 300 : jr(3,2)=n_y*dpiu_r 301 : c------------------- 302 : jr(3,3)=n_y*dpiv_r 303 : c------------------- 304 : jr(3,4)=n_y*dpip_r 305 : c------------------------------------------------------------- 306 : c ------ f44444444444444444444444444444 --------- 307 : c------------------------------------------------------------- 308 : jl(4,1)=0.0D0 309 : c--------------------- 310 : jl(4,2)=0.0D0 311 : c--------------------- 312 : jl(4,3)=0.0D0 313 : c--------------------- 314 : jl(4,4)=0.0D0 315 : c---------------------------------------------------------- 316 : c---------------------------------------------------------- 317 : jr(4,1)=0.0D0 318 : c--------------------- 319 : jr(4,2)=0.0D0 320 : c--------------------- 321 : jr(4,3)=0.0D0 322 : c--------------------- 323 : jr(4,4)=0.0D0 324 : c---------------------------------------------------------- 325 : tcoef=n_x*t_y+t_x*n_y 326 : bcoef=n_x*t_y-t_x*n_y 327 : c----------------------------- 328 : c----------------------------------- 329 : do 11 i=1,4 330 : jl(i,1)=jl(i,1)+jr(i,1) 331 : jl(i,2)=jl(i,2)+jr(i,2)*(-tcoef/bcoef)+ 332 : & jr(i,3)*2.0d0*n_x*t_x/bcoef 333 : jl(i,3)=jl(i,3)+jr(i,2)*(-2.0d0*n_y*t_y/bcoef)+ 334 : & jr(i,3)*tcoef/bcoef 335 : jl(i,4)=jl(i,4)+jr(i,4) 336 : 11 continue 337 : c-------------------------------------- 338 : do 12 i=1,4 339 : dfrun(i)=jl(2,i)*n_x+jl(3,i)*n_y 340 : 12 continue 341 : c---------------------------------------------------------------------- 342 : return 343 : end 344 : 345 : 346 : 347 : 348 : 349 : 350 : 351 : 352 : 353 : 354 : 355 : 356 : 357 : 358 : 359 : 360 : 361 : 362 :
© Cast3M 2003 - All rights reserved.
Disclaimer