1 subroutine getDiff( dwl,yl, shape, xmudmi, xl, rmu, rho) 2c----------------------------------------------------------------------- 3c compute and add the contribution of the turbulent 4c eddy viscosity to the molecular viscosity. 5c----------------------------------------------------------------------- 6 use turbSA 7 include "common.h" 8 9 real*8 yl(npro,nshl,ndof), rmu(npro), xmudmi(npro,ngauss), 10 & shape(npro,nshl), rho(npro), 11 & dwl(npro,nshl), sclr(npro), 12 & xl(npro,nenl,nsd) 13 integer n, e 14 15 real*8 epsilon_ls, kay, epsilon 16 & h_param, prop_blend(npro),test_it(npro) 17c 18c 19c.... get the material properties (2 fluid models will need to determine 20c the "interpolated in phase space" properties....constant for now. 21c two options exist in the interpolation 1) smooth (recommended) 22c interpolation of nodal data, 2) discontinuous "sampling" phase at 23c quadrature points. 24c 25CAD 26CAD prop_blend is a smoothing function to avoid possible large density 27CAD gradients, e.g., water and air flows where density ratios can approach 28CAD 1000. 29CAD 30CAD epsilon_ls is an adjustment to control the width of the band over which 31CAD the properties are blended. 32 33 34 35 if (iLSet .eq. 0)then 36 37 rho = datmat(1,1,1) ! single fluid model, i.e., only 1 density 38 rmu = datmat(1,2,1) 39 40 else ! two fluid properties used in this model 41 42! Smooth the tranistion of properties for a "distance" of epsilon_ls 43! around the interface. Here "distance" is define as the value of the 44! levelset function. If the levelset function is properly defined, 45! this is the true distance normal from the front. Of course, the 46! distance is in a driection normal to the front. 47 48 Sclr = zero 49 isc=abs(iRANS)+6 50 do n = 1, nshl 51 Sclr = Sclr + shape(:,n) * yl(:,n,isc) 52 enddo 53 do i= 1, npro 54 if (sclr(i) .lt. - epsilon_ls)then 55 prop_blend(i) = zero 56 elseif (abs(sclr(i)) .le. epsilon_ls)then 57 prop_blend(i) = 0.5*(one + Sclr(i)/epsilon_ls + 58 & (sin(pi*Sclr(i)/epsilon_ls))/pi ) 59 elseif (sclr(i) .gt. epsilon_ls) then 60 prop_blend(i) = one 61 endif 62 enddo 63c 64 rho = datmat(1,1,2) + (datmat(1,1,1)-datmat(1,1,2))*prop_blend 65 rmu = datmat(1,2,2) + (datmat(1,2,1)-datmat(1,2,2))*prop_blend 66 67 endif 68 69CAD At this point we have a rho that is bounded by the two values for 70CAD density 1, datmat(1,1,1), the fluid, and density 2, datmat(1,1,2) 71CAD the gas 72 73c 74c The above approach evaluates all intermediate quantities at the 75c quadrature point, then combines them to form the needed quantities there. 76c 1 alternative is calculating all quanties (only rho here) at the nodes and 77c then interpolating the result to the quadrature points. If this is done, 78c do not forget to do the same thing for rou in e3b!!! 79c ^^^^^^^^^^ 80c |||||||||| 81c WARNING 82c 83c.... dynamic model 84c 85 if (iLES .gt. 0 .and. iRANS.eq.0) then ! simple LES 86 rmu = rmu + xmudmi(:,intp) 87 else if (iRANS.lt.0) then 88 if (iRANS .eq. -1) then ! RANS (Spalart-Allmaras) 89 call AddEddyViscSA(yl, shape, rmu) 90 else if(iRANS.eq.-2) then ! RANS-KE 91 sigmaInv=1.0 ! use full eddy viscosity for flow equations 92 call AddEddyViscKE(yl, dwl, shape, rho, sigmaInv, rmu) 93 endif 94 if (iLES.gt.0) then ! this is DES so we have to blend in 95 ! xmudmi based on max edge length of 96 ! element 97 call EviscDESIC (xl,rmu,xmudmi) 98 endif 99 endif ! check for LES or RANS 100c 101 return 102 end 103 104 subroutine EviscDESIC(xl,xmut,xmudmi) 105 106 include "common.h" 107 real*8 xmut(npro),xl(npro,nenl,nsd),xmudmi(npro,ngauss) 108 109 110 do i=1,npro 111 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 112 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 113 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 114 emax=max(dx,max(dy,dz)) 115 if(emax.lt.eles) then ! pure les 116 xmut(i)=xmudmi(i,intp) 117 else if(emax.lt.two*eles) then ! blend 118 xi=(emax-eles)/(eles) 119 xmut(i)=xi*xmut(i)+(one-xi)*(xmudmi(1,intp)+datmat(1,2,2)) 120 endif ! leave at RANS value as edge is twice pure les 121 enddo 122 !this was made messy by the addEddyVisc routines Must clean up later. 123 124 125 126 return 127 end 128 129 subroutine getdiffsclr(shape, dwl, yl, diffus) 130 131 use turbSA 132 use turbKE ! access to KE model constants 133 include "common.h" 134 135 real*8 diffus(npro), rho(npro) 136 real*8 yl(npro,nshl,ndof), dwl(npro,nshl), shape(npro,nshl) 137 integer n, e 138 rho(:) = datmat(1,1,1) ! single fluid model, i.e., only 1 density 139 if(isclr.eq.0) then ! solving the temperature equation 140 diffus(:) = datmat(1,4,1) 141 else if(iRANS.eq.-1) then ! solving SA model 142 diffus(:) = datmat(1,2,1) 143 call AddSAVar(yl, shape, diffus) 144 else if(iRANS.eq.-2)then ! solving KE model 145 diffus(:) = datmat(1,2,1) 146 if(isclr.eq.2) then 147 sigmaInv=1.0/ke_sigma ! different eddy viscosity for epsilon 148 else 149 sigmaInv=1.0 ! full eddy viscosity for solving kay equation 150 endif 151 call AddEddyViscKE(yl, dwl, shape, rho, sigmaInv, diffus) 152 else ! solving scalar advection diffusion equations 153 diffus = scdiff(isclr) 154 endif 155c 156 return 157 end 158 159 function ev2sa(xmut,rm,cv1) 160 implicit none 161 real*8 err,ev2sa,rm,cv1,f,dfds,rat,efac 162 real*8 pt5,kappa,B,xmut,chi3,denom,cv1_3 163 integer iter 164 pt5=0.5 165 err=1.0d-6 166 ev2sa=rm*cv1*1.2599 ! inflection point chi=cv1*cuberoot(2) 167 kappa=0.4 168c$$$ B=5.5 169 efac=0.1108 ! exp(-kappa*B) 170 do iter=1,50 171 chi3=ev2sa/rm 172 chi3=chi3*chi3*chi3 173 cv1_3=cv1**3 174 denom=chi3+cv1_3 175 176 f=ev2sa*chi3/denom - xmut 177 dfds=chi3*(chi3+4.0*cv1_3)/(denom**2) 178 rat=-f/dfds 179 ev2sa=ev2sa+rat 180 if(abs(rat).le.err) goto 20 181 enddo 182 write(*,*)'ev2sa failed to converge' 183 write(*,*) 'dfds, rat, ev2sa, mu' 184 write(*,*) dfds,rat,ev2sa,rm 185 20 continue 186 return 187 end 188c 189 190 191 subroutine AddEddyViscSA(yl,shape,rmu) 192 use turbSA 193 include "common.h" 194c INPUTS 195 double precision, intent(in), dimension(npro,nshl,ndof) :: 196 & yl 197 double precision, intent(in), dimension(npro,nshl) :: 198 & shape 199c INPUT-OUTPUTS 200 double precision, intent(inout), dimension(npro) :: 201 & rmu 202c LOCALS 203 logical, dimension(nshl) :: 204 & wallnode 205 integer e, n 206 double precision xki, xki3, fv1, evisc 207c 208c Loop over elements in this block 209 do e = 1, npro 210c assume no wall nodes on this element 211 wallnode(:) = .false. 212 if(itwmod.eq.-2) then ! effective viscosity 213c mark the wall nodes for this element, if there are any 214 do n = 1, nshl 215 u1=yl(e,n,2) 216 u2=yl(e,n,3) 217 u3=yl(e,n,4) 218 if((u1.eq.zero).and.(u2.eq.zero).and.(u3.eq.zero)) 219 & then 220 wallnode(n)=.true. 221 endif 222 enddo 223 endif 224c 225 if( any(wallnode(:)) ) then 226c if there are wall nodes for this elt, then we are using effective- 227c viscosity near-wall modeling, and eddy viscosity has been stored 228c at the wall nodes in place of the spalart-allmaras variable; the 229c eddy viscosity for the whole element is taken to be the avg of the 230c wall values 231 evisc = zero 232 nwnode=0 233 do n = 1, nshl 234 if(wallnode(n)) then 235 evisc = evisc + yl(e,n,6) 236 nwnode = nwnode + 1 237 endif 238 enddo 239 evisc = evisc/nwnode 240 rmu(e) = rmu(e) + abs(evisc) 241c this is what we would use instead of the above if we were allowing 242c the eddy viscosity to vary through the element based on non-wall nodes 243c$$$ evisc = zero 244c$$$ Turb = zero 245c$$$ do n = 1, nshl 246c$$$ if(wallmask(n).eq.1) then 247c$$$ evisc = evisc + shape(e,n) * yl(e,n,6) 248c$$$ else 249c$$$ Turb = Turb + shape(e,n) * yl(e,n,6) 250c$$$ endif 251c$$$ enddo 252c$$$ xki = abs(Turb)/rmu(e) 253c$$$ xki3 = xki * xki * xki 254c$$$ fv1 = xki3 / (xki3 + saCv1P3) 255c$$$ rmu(e) = rmu(e) + fv1*abs(Turb) 256c$$$ rmu(e) = rmu(e) + abs(evisc) 257 else 258c else one of the following is the case: 259c using effective-viscosity, but no wall nodes on this elt 260c using slip-velocity 261c using no model; walls are resolved 262c in all of these cases, eddy viscosity is calculated normally 263 Turb = zero 264 do n = 1, nshl 265 Turb = Turb + shape(e,n) * yl(e,n,6) 266 enddo 267 xki = abs(Turb)/rmu(e) 268 xki3 = xki * xki * xki 269 fv1 = xki3 / (xki3 + saCv1P3) 270 rmu(e) = rmu(e) + fv1*abs(Turb) 271 endif 272 enddo ! end loop over elts 273 return 274 end subroutine AddEddyViscSA 275 276 277 278 subroutine AddSAVar(yl,shape,rmu) 279 use turbSA 280 include "common.h" 281c INPUTS 282 double precision, intent(in), dimension(npro,nshl,ndof) :: 283 & yl 284 double precision, intent(in), dimension(npro,nshl) :: 285 & shape 286c INPUT-OUTPUTS 287 double precision, intent(inout), dimension(npro) :: 288 & rmu 289c LOCALS 290 logical, dimension(nshl) :: 291 & wallnode 292 integer e, n 293 double precision savar, savarw 294c Loop over elements in this block 295 do e = 1, npro 296c assume no wall nodes on this element 297 wallnode(:) = .false. 298 if(itwmod.eq.-2) then ! effective viscosity 299c mark the wall nodes for this element, if there are any 300 do n = 1, nshl 301 u1=yl(e,n,2) 302 u2=yl(e,n,3) 303 u3=yl(e,n,4) 304 if((u1.eq.zero).and.(u2.eq.zero).and.(u3.eq.zero)) 305 & then 306 wallnode(n)=.true. 307 endif 308 enddo 309 endif 310c 311 savar=zero 312 do n = 1, nshl 313 if( wallnode(n) ) then 314c if wallmask was set, we're using effective-viscosity wall-model and 315c this node is on a wall. Eddy viscosity has been stored at the wall 316c nodes in place of the S-A variable, so we must convert it 317 savarw = ev2sa(yl(e,n,6),datmat(1,2,1),saCv1) 318 savar = savar + shape(e,n) * savarw 319 else 320c if wallmask wasn't set, then one of the following is the case: 321c using effective-viscosity, but this isn't a wall node 322c using slip-velocity 323c using no wall model; wall is resolved 324c in all these cases, the S-A variable is calculated normally 325 savar = savar + shape(e,n) * yl(e,n,6) 326 endif 327 enddo 328 rmu(e)=datmat(1,2,1) 329 rmu(e) = (rmu(e) + abs(savar)) * saSigmaInv 330 enddo ! end loop over elts 331 return 332 end subroutine AddSAVar 333 334 335 336 subroutine AddEddyViscKE(yl, dwl, shape, rho, sigmaInv, rmu) 337 use turbKE ! access to KE model constants 338 include "common.h" 339c INPUTS 340 double precision, intent(in), dimension(npro,nshl,ndof) :: 341 & yl 342 double precision, intent(in), dimension(npro,nshl) :: 343 & shape, dwl 344 double precision, intent(in), dimension(npro) :: 345 & rho 346 double precision sigmaInv 347c INPUT-OUTPUTS 348 double precision, intent(inout), dimension(npro) :: 349 & rmu 350c LOCALS 351 double precision eviscKE, kay, epsilon, dw, CmuKE 352 double precision epsInv, Rey, Ret, RetInv, tmp1, fmuKE 353 integer e,n 354c 355 do e = 1, npro 356 kay = 0.0 357 epsilon = 0.0 358 dw = 0.0 359 do n = 1, nshl 360 kay = kay + shape(e,n)*yl(e,n,6) 361 epsilon = epsilon + shape(e,n)*yl(e,n,7) 362 dw = dw + shape(e,n)*dwl(e,n) 363 enddo 364 kay = abs(kay) 365 if(kay.lt.1.0e-32) kay=0.0 366 epsInv = 0 367 if ( abs(epsilon) .gt.1.e-32) then 368 epsInv = 1. / abs(epsilon) 369 endif 370 371 Rey = sqrt(kay) * dw * rho(e) / rmu(e) 372 Ret = kay*kay * epsInv * rho(e) / rmu(e) 373 RetInv = 0 374 if(Ret.lt.1.d100.AND.Ret.gt.zero) RetInv=1./Ret 375 tmp1 = exp(-0.0165*Rey) 376 fmuKE = (1. -tmp1) ** 2 * (1.+20.5*RetInv) ! fmu of Lam-Bremhorst 377 378 eviscKE=rho(e)*ke_C_mu*fmuKE*kay*kay*epsInv 379 380 rmu(e) = rmu(e) + eviscKE*sigmaInv 381 enddo 382 return 383 end subroutine AddEddyViscKE 384 385 386