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 54c 55c 56c.... constant viscosity 57c 58 if (matflg(2,1) .eq. 0) then 59c 60 if (iLSet .ne. 0)then !two fluid properties used in this model 61 Sclr = zero 62 isc=abs(iRANS)+6 63 do n = 1, nshl 64 Sclr = Sclr + shp(:,n) * ycl(:,n,isc) 65 enddo 66 test_it = 0.5*(one + Sclr/epsilon_ls + 67 & (sin(pi*Sclr/epsilon_ls))/pi ) 68 69 prop_blend = max( min(test_it(:), one ), zero ) 70 rmu = datmat(1,2,2) + (datmat(1,2,1)-datmat(1,2,2)) 71 & *prop_blend 72 elseif(irampViscOutlet.eq.1)then ! increase viscosity near outlet 73c.............ramp rmu near outlet (for a NGC geometry) 74 xx=zero 75 do n=1,nenl 76 xx(:)=xx(:) + shp(:,n) * xl(:,n,1) 77 enddo 78 fmax=10.0 79 where(xx(:).le. 0.42) !healfway btwn AIP and exit 80 rmu(:)=datmat(1,2,1) 81 elsewhere(xx(:).ge. 0.75) !2/3 of the way to the exit 82 rmu(:)=fmax*datmat(1,2,1) 83 elsewhere 84 rmu(:)= datmat(1,2,1)*( 85 & (55.65294821-55.65294821*fmax)*xx(:)*xx(:)*xx(:) 86 & +(-97.67092412+97.67092412*fmax)*xx(:)*xx(:) 87 & +(52.59203606-52.59203606*fmax)*xx(:) 88 & -7.982719760+8.982719760*fmax) 89 endwhere 90 else ! constant viscosity 91 rmu = datmat(1,2,1) 92 endif 93c 94 else 95c 96c.... generalized Sutherland viscosity 97c 98 rmu = datmat(1,2,1) * (T/datmat(2,2,1))*sqrt(T/datmat(2,2,1)) 99 & * ( datmat(2,2,1) + datmat(3,2,1) ) / (T + datmat(3,2,1)) 100c 101 endif 102c 103c.... calculate the second viscosity coefficient 104c 105 if (matflg(3,1) .eq. 0) then 106 rlm = -pt66 * rmu 107 else 108 rlm = (datmat(1,3,1) - pt66) * rmu 109 endif 110c 111c.... calculate the remaining quantities 112c 113 con = rmu * cp / pr 114c 115c-------------Eddy Viscosity Calculation----------------- 116c 117c.... dynamic model 118c 119 if (iLES .gt. 0. and. iRANS.eq.0) then ! simple LES 120 xmut = xmudmi(:,intp) 121 else if (iRANS .eq. 0 .and. iLES.eq.0 ) then !DNS 122 xmut = zero 123 else if (iRANS .lt. 0) then ! calculate RANS viscosity 124c 125c.... RANS 126c 127 do e = 1, npro 128 wallmask = 0 129 if(itwmod.eq.-2) then ! effective viscosity 130c mark the wall nodes for this element, if there are any 131 do n = 1, nshl 132c 133c note that we are using ycl here so that means that these 134c terms are not perturbed for MFG difference and therefore 135c NOT in the LHS. As they only give the evisc near the wall 136c I doubt this is a problem. 137c 138 u1=ycl(e,n,2) 139 u2=ycl(e,n,3) 140 u3=ycl(e,n,4) 141 if((u1.eq.zero).and.(u2.eq.zero).and.(u3.eq.zero)) 142 & then 143 wallmask(n)=1 144 endif 145 enddo 146 endif 147c 148 if( any(wallmask.eq.1) ) then 149c if there are wall nodes for this elt in an effective-viscosity wall 150c modeled case,then eddy viscosity has been stored at the wall nodes 151c in place of the spalart-allmaras variable; the eddy viscosity for 152c the whole element is taken to be the avg of wall values 153 evisc = zero 154 nwnode=0 155 do n = 1, nshl 156 if(wallmask(n).eq.1) then 157 evisc = evisc + ycl(e,n,6) 158 nwnode = nwnode + 1 159 endif 160 enddo 161 evisc = evisc/nwnode 162 xmut(e)= abs(evisc) 163c this is what we would use instead of the above if we were allowing 164c the eddy viscosity to vary through the element based on non-wall nodes 165c$$$ evisc = zero 166c$$$ Turb = zero 167c$$$ do n = 1, nshl 168c$$$ if(wallmask(n).eq.1) then 169c$$$ evisc = evisc + shape(e,n) * ycl(e,n,6) 170c$$$ else 171c$$$ Turb = Turb + shape(e,n) * ycl(e,n,6) 172c$$$ endif 173c$$$ enddo 174c$$$ xki = abs(Turb)/rmu(e) 175c$$$ xki3 = xki * xki * xki 176c$$$ fv1 = xki3 / (xki3 + saCv1P3) 177c$$$ rmu(e) = rmu(e) + fv1*abs(Turb) 178c$$$ rmu(e) = rmu(e) + abs(evisc) 179 else 180c else one of the following is the case: 181c using effective-viscosity, but no wall nodes on this elt 182c using slip-velocity 183c using no model; walls are resolved 184c in all of these cases, eddy viscosity is calculated normally 185 savar = zero 186 do n = 1, nshl 187 savar = savar + shp(e,n) * ycl(e,n,6) 188 enddo 189 xki = abs(savar)/rmu(e) 190 xki3 = xki * xki * xki 191 fv1 = xki3 / (xki3 + saCv1P3) 192 xmut(e) = fv1*abs(savar) 193 endif 194 enddo ! end loop over elts 195 196 if (iLES.gt.0) then ! this is DES so we have to blend in 197 ! xmudmi based on max edge length of 198 ! element 199 call EviscDES (xl,xmut,xmudmi) 200 endif 201 endif ! check for LES or RANS 202 203 rlm = rlm - pt66*xmuT 204 rmu = rmu + xmuT 205 rlm2mu = rlm + two * rmu 206 con = con + xmuT*cp/pr 207c 208c.... return 209c 210 return 211 end 212c 213c 214c 215 subroutine getDiffSclr (T, cp, rmu, rlm, 216 & rlm2mu, con, rho, Sclr) 217c 218c---------------------------------------------------------------------- 219c 220c This routine calculates the fluid material properties. 221c 222c input: 223c T (npro) : temperature 224c cp (npro) : specific heat at constant pressure 225c 226c output: 227c rmu (npro) : Mu 228c rlm (npro) : Lambda 229c rlm2mu (npro) : Lambda + 2 Mu 230c con (npro) : Conductivity 231c 232c Note: material type flags 233c matflg(2): 234c eq. 0, constant viscosity 235c eq. 1, generalized Sutherland viscosity 236c matflg(3): 237c eq. 0, Stokes approximation 238c eq. 1, shear proportional bulk viscosity 239c 240c 241c Farzin Shakib, Winter 1987. 242c Zdenek Johan, Winter 1991. (Fortran 90) 243c---------------------------------------------------------------------- 244c 245 include "common.h" 246c 247 dimension T(npro), cp(npro), 248 & rmu(npro), rlm(npro), 249 & rlm2mu(npro), con(npro), 250 & rho(npro), Sclr(npro) 251 252 253 254c 255c 256c.... constant viscosity 257c 258 if (matflg(2,1) .eq. 0) then 259c 260 rmu = datmat(1,2,1) 261c 262 else 263c 264c.... generalized Sutherland viscosity 265c 266 rmu = datmat(1,2,1) * (T/datmat(2,2,1))*sqrt(T/datmat(2,2,1)) 267 & * ( datmat(2,2,1) + datmat(3,2,1) ) / (T + datmat(3,2,1)) 268c 269 endif 270c 271*************************check**************************** 272c if (iRANS(1).lt.zero) then 273c rmu = saSigmaInv*rho*((rmu/rho)+Sclr) 274c endif 275c This augmentation of viscosity is performed in e3viscsclr 276c The Spalart -Allmaras model will need molecular viscosity 277c in subsequent calculations. 278c.... calculate the second viscosity coefficient 279c 280 if (matflg(3,1) .eq. 0) then 281c 282 rlm = -pt66 * rmu 283c 284 else 285c 286 rlm = (datmat(1,3,1) - pt66) * rmu 287c 288 endif 289c 290c.... calculate the remaining quantities 291c 292 293 294 295 rlm2mu = rlm + two * rmu 296 con = rmu * cp / pr 297 298 299 300c 301c.... return 302c 303 return 304 end 305 306 subroutine EviscDES(xl,xmut,xmudmi) 307 308 include "common.h" 309 real*8 xmut(npro),xl(npro,nenl,nsd),xmudmi(npro,ngauss) 310 311 312 do i=1,npro 313 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 314 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 315 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 316 emax=max(dx,max(dy,dz)) 317 if(emax.lt.eles) then ! pure les 318 xmut(i)=xmudmi(i,intp) 319 else if(emax.lt.two*eles) then ! blend 320 xi=(emax-eles)/(eles) 321 xmut(i)=xi*xmut(i)+(one-xi)*xmudmi(1,intp) 322 endif ! leave at RANS value as edge is twice pure les 323 enddo 324 325 return 326 end 327