1 subroutine getDiff (T, cp, rho, ycl, 2 & rmu, rlm, rlm2mu, con, shp, 3 & xmudmi, xl) 4 5c---------------------------------------------------------------------- 6c 7c This routine calculates the fluid material properties. 8c 9c input: 10c T (npro) : temperature 11c cp (npro) : specific heat at constant pressure 12c ************************************************************** 13c rho (npro) : density 14c ycl (npro,nshl,ndof): Y variables 15c shp (npro,nshl) : element shape-functions 16c ************************************************************* 17c output: 18c rmu (npro) : Mu 19c rlm (npro) : Lambda 20c rlm2mu (npro) : Lambda + 2 Mu 21c con (npro) : Conductivity 22c 23c Note: material type flags 24c matflg(2): 25c eq. 0, constant viscosity 26c eq. 1, generalized Sutherland viscosity 27c matflg(3): 28c eq. 0, Stokes approximation 29c eq. 1, shear proportional bulk viscosity 30c 31c 32c Farzin Shakib, Winter 1987. 33c Zdenek Johan, Winter 1991. (Fortran 90) 34c---------------------------------------------------------------------- 35c 36 use turbSA 37 use pointer_data 38 include "common.h" 39c 40 dimension T(npro), cp(npro), 41 & rho(npro), Sclr(npro), 42 & rmu(npro), rlm(npro), 43 & rlm2mu(npro), con(npro), 44 & ycl(npro,nshl,ndof), shp(npro,nshl), 45 & xmudmi(npro,ngauss), xl(npro,nenl,nsd), 46 & xx(npro) 47c 48 dimension xmut(npro) 49 real*8 prop_blend(npro),test_it(npro) 50 51 integer n, e 52 integer wallmask(nshl) 53 real*8 xki, xki3, fv1, evisc, lvisc 54 xx=zero 55 do n=1,nenl 56 xx(:)=xx(:) + shp(:,n) * xl(:,n,1) 57 enddo 58c 59c 60c.... constant viscosity 61c 62 if (matflg(2,1) .eq. 0) then 63c 64 if (iLSet .ne. 0)then !two fluid properties used in this model 65 Sclr = zero 66 isc=abs(iRANS)+6 67 do n = 1, nshl 68 Sclr = Sclr + shp(:,n) * ycl(:,n,isc) 69 enddo 70 test_it = 0.5*(one + Sclr/epsilon_ls + 71 & (sin(pi*Sclr/epsilon_ls))/pi ) 72 73 prop_blend = max( min(test_it(:), one ), zero ) 74 rmu = datmat(1,2,2) + (datmat(1,2,1)-datmat(1,2,2)) 75 & *prop_blend 76 elseif(irampViscOutlet.eq.1)then ! increase viscosity near outlet 77c.............ramp rmu near outlet (for a Duct geometry) 78 fmax=10.0 79! fmax=2000.0 80 81c if (myrank .eq. master)then 82c write(*,*) 'viscosity', datmat(1,2,1) 83c endif 84 85c... geometry6 86 if(iDuctgeometryType .eq. 6)then 87 88c if (myrank .eq. master) write(*,*) 'getdiff, geometry6' 89 90 where(xx(:).le. 0.42) !halfway btwn AIP and exit 91 rmu(:)= datmat(1,2,1) 92 elsewhere(xx(:).ge. 0.75) !2/3 of the way to the exit 93 rmu(:)=fmax*datmat(1,2,1) 94 elsewhere 95 rmu(:)= datmat(1,2,1)*( 96 & (55.65294821-55.65294821*fmax)*xx(:)*xx(:)*xx(:) 97 & +(-97.67092412+97.67092412*fmax)*xx(:)*xx(:) 98 & +(52.59203606-52.59203606*fmax)*xx(:) 99 & -7.982719760+8.982719760*fmax) 100 101 endwhere 102 103c do i = 1,npro 104c if(xx(i) .lt. 0.75 .and. xx(i) .gt. 0.42) 105c & write(*,*) xx(i), rmu(i) 106c enddo 107 108c... geometry8 109 elseif (iDuctgeometryType .eq. 8)then 110 xstart=1.5 !1.6*4.5*0.0254+0.85*0.5 111 xmidwy=2.0 !1.6*4.5*0.0254+0.85*1.0 112 where(xx(:).le.xstart) 113 rmu(:)=datmat(1,2,1) 114 elsewhere(xx(:).ge.xmidwy) 115 rmu(:)=fmax*datmat(1,2,1) 116 elsewhere 117 rmu(:)=datmat(1,2,1) 118 & *(1+(fmax-1)*(0.5+ 119 & tanh(5*(xx(:)-(xmidwy+xstart)/2)/(xmidwy-xstart))/2)) 120 endwhere 121 122 123 124 endif 125c..................................................... 126 else ! constant viscosity 127 rmu = datmat(1,2,1) 128 endif 129! 130! boundary layer thickening via molecular viscosity 131! 132 scaleCntrl=1.0 133 Lvisc=0.2 134 xbltb=-0.2159-two*Lvisc 135 xblte=-0.2159-Lvisc 136 where((xx(:).gt.xbltb) .and. (xx(:).lt.xblte)) 137 rmu(:)=scaleCntrl*datmat(1,2,1) 138 endwhere 139 140! xvisc1 = -0.3048 141! xvisc2 = -0.2159 142! where(xx(:).lt.xvisc1) 143! rmu(:)=scaleCntrl*datmat(1,2,1) 144! elsewhere(xx(:).gt.xvisc1 .and. xx(:).lt.xvisc2) 145! rmu(:)=( scaleCntrl - (scaleCntrl - 1)* 146! & (xx(:) - xvisc1)/(xvisc2 - xvisc1))*datmat(1,2,1) 147! endwhere 148 149 !if(myrank.eq.master) then 150 ! write(*,*) 'adjusting viscosity in region by ', scaleCntrl 151 !endif 152 else 153c 154c.... generalized Sutherland viscosity 155c 156 rmu = datmat(1,2,1) * (T/datmat(2,2,1))*sqrt(T/datmat(2,2,1)) 157 & * ( datmat(2,2,1) + datmat(3,2,1) ) / (T + datmat(3,2,1)) 158c 159 endif 160c 161c.... calculate the second viscosity coefficient 162c 163 if (matflg(3,1) .eq. 0) then 164 rlm = -pt66 * rmu 165 else 166 rlm = (datmat(1,3,1) - pt66) * rmu 167 endif 168c 169c.... calculate the remaining quantities 170c 171 con = rmu * cp / pr 172c 173c-------------Eddy Viscosity Calculation----------------- 174c 175c.... dynamic model 176c 177 if (iLES .gt. 0. and. iRANS.eq.0) then ! simple LES 178 xmut = xmudmi(:,intp) 179 else if (iRANS .eq. 0 .and. iLES.eq.0 ) then !DNS 180 xmut = zero 181 else if (iRANS .lt. 0) then ! calculate RANS viscosity 182c 183c.... RANS 184c 185 do e = 1, npro 186 wallmask = 0 187 if(itwmod.eq.-2) then ! effective viscosity 188c mark the wall nodes for this element, if there are any 189 do n = 1, nshl 190c 191c note that we are using ycl here so that means that these 192c terms are not perturbed for MFG difference and therefore 193c NOT in the LHS. As they only give the evisc near the wall 194c I doubt this is a problem. 195c 196 u1=ycl(e,n,2) 197 u2=ycl(e,n,3) 198 u3=ycl(e,n,4) 199 if((u1.eq.zero).and.(u2.eq.zero).and.(u3.eq.zero)) 200 & then 201 wallmask(n)=1 202 endif 203 enddo 204 endif 205c 206 if( any(wallmask.eq.1) ) then 207c if there are wall nodes for this elt in an effective-viscosity wall 208c modeled case,then eddy viscosity has been stored at the wall nodes 209c in place of the spalart-allmaras variable; the eddy viscosity for 210c the whole element is taken to be the avg of wall values 211 evisc = zero 212 nwnode=0 213 do n = 1, nshl 214 if(wallmask(n).eq.1) then 215 evisc = evisc + ycl(e,n,6) 216 nwnode = nwnode + 1 217 endif 218 enddo 219 evisc = evisc/nwnode 220 xmut(e)= abs(evisc) 221c this is what we would use instead of the above if we were allowing 222c the eddy viscosity to vary through the element based on non-wall nodes 223c$$$ evisc = zero 224c$$$ Turb = zero 225c$$$ do n = 1, nshl 226c$$$ if(wallmask(n).eq.1) then 227c$$$ evisc = evisc + shape(e,n) * ycl(e,n,6) 228c$$$ else 229c$$$ Turb = Turb + shape(e,n) * ycl(e,n,6) 230c$$$ endif 231c$$$ enddo 232c$$$ xki = abs(Turb)/rmu(e) 233c$$$ xki3 = xki * xki * xki 234c$$$ fv1 = xki3 / (xki3 + saCv1P3) 235c$$$ rmu(e) = rmu(e) + fv1*abs(Turb) 236c$$$ rmu(e) = rmu(e) + abs(evisc) 237 else 238c else one of the following is the case: 239c using effective-viscosity, but no wall nodes on this elt 240c using slip-velocity 241c using no model; walls are resolved 242c in all of these cases, eddy viscosity is calculated normally 243 savar = zero 244 do n = 1, nshl 245 savar = savar + shp(e,n) * ycl(e,n,6) 246 enddo 247 xki = abs(savar)/rmu(e) 248 xki3 = xki * xki * xki 249 fv1 = xki3 / (xki3 + saCv1P3) 250 xmut(e) = fv1*abs(savar) 251 endif 252 enddo ! end loop over elts 253 254 if (iLES.gt.0) then ! this is DES so we have to blend in 255 ! xmudmi based on max edge length of 256 ! element 257 call EviscDES (xl,xmut,xmudmi) 258 endif 259 endif ! check for LES or RANS 260 261 rlm = rlm - pt66*xmuT 262 rmu = rmu + xmuT 263 rlm2mu = rlm + two * rmu 264 con = con + xmuT*cp/pr 265c 266c.... return 267c 268 return 269 end 270c 271c 272c 273 subroutine getDiffSclr (T, cp, rmu, rlm, 274 & rlm2mu, con, rho, Sclr) 275c 276c---------------------------------------------------------------------- 277c 278c This routine calculates the fluid material properties. 279c 280c input: 281c T (npro) : temperature 282c cp (npro) : specific heat at constant pressure 283c 284c output: 285c rmu (npro) : Mu 286c rlm (npro) : Lambda 287c rlm2mu (npro) : Lambda + 2 Mu 288c con (npro) : Conductivity 289c 290c Note: material type flags 291c matflg(2): 292c eq. 0, constant viscosity 293c eq. 1, generalized Sutherland viscosity 294c matflg(3): 295c eq. 0, Stokes approximation 296c eq. 1, shear proportional bulk viscosity 297c 298c 299c Farzin Shakib, Winter 1987. 300c Zdenek Johan, Winter 1991. (Fortran 90) 301c---------------------------------------------------------------------- 302c 303 include "common.h" 304c 305 dimension T(npro), cp(npro), 306 & rmu(npro), rlm(npro), 307 & rlm2mu(npro), con(npro), 308 & rho(npro), Sclr(npro) 309 310 311 312c 313c 314c.... constant viscosity 315c 316 if (matflg(2,1) .eq. 0) then 317c 318 rmu = datmat(1,2,1) 319c 320 else 321c 322c.... generalized Sutherland viscosity 323c 324 rmu = datmat(1,2,1) * (T/datmat(2,2,1))*sqrt(T/datmat(2,2,1)) 325 & * ( datmat(2,2,1) + datmat(3,2,1) ) / (T + datmat(3,2,1)) 326c 327 endif 328c 329*************************check**************************** 330c if (iRANS(1).lt.zero) then 331c rmu = saSigmaInv*rho*((rmu/rho)+Sclr) 332c endif 333c This augmentation of viscosity is performed in e3viscsclr 334c The Spalart -Allmaras model will need molecular viscosity 335c in subsequent calculations. 336c.... calculate the second viscosity coefficient 337c 338 if (matflg(3,1) .eq. 0) then 339c 340 rlm = -pt66 * rmu 341c 342 else 343c 344 rlm = (datmat(1,3,1) - pt66) * rmu 345c 346 endif 347c 348c.... calculate the remaining quantities 349c 350 351 352 353 rlm2mu = rlm + two * rmu 354 con = rmu * cp / pr 355 356 357 358c 359c.... return 360c 361 return 362 end 363 364 subroutine EviscDES(xl,xmut,xmudmi) 365 366 include "common.h" 367 real*8 xmut(npro),xl(npro,nenl,nsd),xmudmi(npro,ngauss) 368 369 370 do i=1,npro 371 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 372 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 373 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 374 emax=max(dx,max(dy,dz)) 375 if(emax.lt.eles) then ! pure les 376 xmut(i)=xmudmi(i,intp) 377 else if(emax.lt.two*eles) then ! blend 378 xi=(emax-eles)/(eles) 379 xmut(i)=xi*xmut(i)+(one-xi)*xmudmi(1,intp) 380 endif ! leave at RANS value as edge is twice pure les 381 enddo 382 383 return 384 end 385