1 subroutine settauw (y, x, BC, 2 & ifath, velbar) 3c 4c---------------------------------------------------------------------- 5c 6c This routine computes the time varying viscous flux for turbulence wall 7c boundary elements. 8c 9c Zdenek Johan, Winter 1991. (Fortran 90) 10c---------------------------------------------------------------------- 11c 12 use pointer_data 13 use turbSA 14 include "common.h" 15c 16 dimension y(nshg,ndofl), x(numnp,nsd), 17 & BC(nshg,ndofBC), 18 & ifath(numnp), velbar(nfath,nflow), 19 & ull(nsd), trx(numnp,nsd), 20 & ullb(nsd), dull(nsd), 21 & evisc(numnp) 22 real*8 outvec(nshg,nsd+nsd+nsd) 23c 24c calculate the traction magnitude for all nodes on the wall 25c 26 trx=zero 27 evisc=zero 28 rm=datmat(1,2,1) 29 do nodw = 1, numnp ! loop over nodes 30 if ( otwn(nodw).ne.nodw ) then ! wall node check 31 if (iLES.gt.0) then 32 ull(:)=velbar(ifath(otwn(nodw)),2:4) ! recall old #ing 33 endif 34 if (iRANS.lt.0) then 35 ull(:)=y(otwn(nodw),1:3) 36 endif 37 ub=ull(1)*wnrm(nodw,1) ! 38 & +ull(2)*wnrm(nodw,2) ! store u.n here for now 39 & +ull(3)*wnrm(nodw,3) ! 40c u_parallel_to_boundary=u-(u.n)n : 41 ull(:)=ull(:)-ub*wnrm(nodw,:) ! ull in flow 42c u_b is || u_ll || 43 ub=sqrt(ull(1)**2+ull(2)**2+ull(3)**2) 44c perp d2wall is (x2-x1).n : 45 dw=abs( 46 & (x(otwn(nodw),1)-x(nodw,1))*wnrm(nodw,1) 47 & +(x(otwn(nodw),2)-x(nodw,2))*wnrm(nodw,2) 48 & +(x(otwn(nodw),3)-x(nodw,3))*wnrm(nodw,3) 49 & ) 50 if(ub.eq.0) then 51 ut=zero 52 twoub=zero 53 else 54 ut=utau(ub,dw,rm,nodw,x(nodw,:)) ! find utau 55 twoub=-ut*ut/ub ! find -tau_w/ub 56 endif 57c 58c for LES ull has the mean-parallel velocity vector. We want the 59c instantaneous-parallel velocity vector 60c 61 if (iLES.gt.0) then 62 ullb=ull 63 ull(:)=y(otwn(nodw),1:3) 64 ubn=ull(1)*wnrm(nodw,1) ! 65 & +ull(2)*wnrm(nodw,2) ! store u.n here for now 66 & +ull(3)*wnrm(nodw,3) ! 67c 68c u_parallel_to_boundary=u-(u.n)n : 69c 70 ull(:)=ull(:)-ubn*wnrm(nodw,:) ! ull in flow 71c 72c hack a limiter into this fluctuating vector. Early transients have 73c huge differences from mean values 74c 75 dull=ull-ullb ! the current vector difference 76 dullm=sqrt(dull(1)*dull(1)+dull(2)*dull(2)+dull(3)*dull(3)) 77 & + 1.0e-9 78 ullbm=ub ! mag of ullb still there 79c 80c limit the magnitude of the difference to a 40% change from the mean. 81c if less than that already we will take the whole difference, otherwise 82c only take a 40% change. 83c 84 85 dullmod=min(one,0.4*ullbm/dullm) 86 ull=ullb+dullmod*dull 87 endif 88 89 trx(nodw,:)=twoub*ull(:) 90 if(itwmod.eq.-2) then ! effective-viscosity 91 tauw=ut*ut 92 BC(nodw,7)=tauw*dw/ub-rm 93 endif 94 if(itwmod.eq.2) then ! effective-viscosity 95c 96c mag of u instantaneous 97c 98 ullm=sqrt(ull(1)*ull(1)+ull(2)*ull(2)+ull(3)*ull(3)) 99 tauw=ut*ut*ullm/ub 100 evisc(nodw)=tauw*dw/ub-rm 101 endif 102 if((itwmod.eq.-1)) then ! slip-velocity RANS 103 up=sqrt( 104 & y(nodw,1)**2 ! 105 & +y(nodw,2)**2 ! flow is auto-|| at boundary 106 & +y(nodw,3)**2 ! 107 & )/ut 108 BC(nodw,7)=savarw(up,rm,saCv1,nodw,x(nodw,:)) 109 endif 110 endif ! wallnode check 111 enddo ! loop over nodes 112c 113c Write the traction vectors to a file restart.4077.n 114c 115c$$$ ilstep=4077 116c$$$ outvec(:,1:3)=trx(:,1:3) 117c$$$ outvec(:,4:6)=0 118c$$$ do i=1,numnp 119c$$$ if(otwn(i).ne. i ) outvec(i,4:6)=y(otwn(i),1:3) 120c$$$ enddo 121c$$$ outvec(:,7:9)=wnrm(:,1:3) 122c$$$ call write_restart(myrank,ilstep,numnp,nsd*3,outvec,outvec) 123c$$$ write(*,*) 'Traction dumped to restart.4077.*' 124c 125c Put traction calculations into BCB 126c 127 do iblk = 1, nelblb 128 iel = lcblkb(1,iblk) 129 nenl = lcblkb(5,iblk) ! no. of vertices per element 130 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 131 ndofl = lcblkb(8,iblk) 132 npro = lcblkb(1,iblk+1) - iel 133c For all elements in this block that lie on a wall, assign the traction 134 do i=1,npro 135 if(btest(miBCB(iblk)%p(i,1),4)) then ! wall elt 136 do j = 1, nenbl 137 if(itwmod.eq.-2) then ! effective-viscosity 138 mBCB(iblk)%p(i,j,3:5)=0.0 139 endif 140 if((itwmod.eq.-1).or.(itwmod.eq.1)) then ! slip-velocity 141 mBCB(iblk)%p(i,j,3:5)=trx(mienb(iblk)%p(i,j),:) 142 endif 143 enddo 144 endif 145 enddo 146 enddo 147 148 149 if(itwmod.eq.2) then !effective viscosity 150c 151c For the elements which touch a modeled wall, 152c modify the eddy viscosity at all quadrature points to be the average 153c of the "optimal" nodal values. This is an element-wise constant. 154c 155 do iblk = 1,nelblk 156 nenl = lcblk(5,iblk) ! no. of vertices per element 157 iel = lcblk(1,iblk) 158 npro = lcblk(1,iblk+1) - iel 159 lcsyst = lcblk(3,iblk) 160 ngauss = nint(lcsyst) 161 do i=1,npro 162 xnave=zero 163 avevisc=zero 164 do j=1,nenl 165 ev = evisc(mien(iblk)%p(i,j)) 166 avevisc=avevisc + ev 167 if(ev.ne.zero) xnave=xnave+1 168 enddo 169 if(xnave.ne.0) mxmudmi(iblk)%p(i,1:ngauss)=avevisc/xnave 170 enddo 171 enddo 172 endif 173c 174c.... end 175c 176 return 177 end 178 179 function utaul(u,y,rm) 180 real*8 utaul 181 utaul=.2 182 return 183 end 184 185 function utau(u,y,rm,nodw,x) 186 implicit none 187 real*8 u,err,utau,yrmi,y,rm,yp,up,kup,f,dfds,rat,efac 188 real*8 kup2,kup3,pt5,sxth,kappa 189 integer iter 190 integer nodw 191 real*8 x(3) 192 pt5=0.5 193 sxth=0.1666666666666667 194 err=1.0d-6 195 utau=0.04 196 yrmi=y/rm 197 kappa=0.4 198c$$$ B=5.5 199 efac=0.1108 ! exp(-kappa*B) 200 do iter=1,500 201 yp=yrmi*utau 202 up=u/utau 203 kup=kappa*up 204 kup2=kup*kup 205 kup3=kup*kup2 206 f= up-yp+efac*(exp(kup)-1.0-kup - kup2*pt5 - kup3*sxth) 207 dfds=up+yp+efac*(exp(kup)*kup-kup - kup2 - kup3*pt5) !/-utau in rat 208 rat=f*utau/dfds ! this is -f/dfds but we did not do 1/(-utau) yet. 209 utau=utau+rat 210 if(abs(rat).le.err) goto 20 211 enddo 212 write(*,*)'utau failed to converge at ',nodw,x(1),x(2),x(3) 213 write(*,*) 'u, dwallperp, mu' 214 write(*,*) u, y, rm 215 write(*,*) 'dfds, rat, utau' 216 write(*,*) dfds,rat,utau 217c stop 218c utau=0.0 219c 220c if the above fails then try a simple no-slip traction 221c 222 utau=sqrt(u/y*rm) 223 224 20 continue 225 return 226 end 227 228c$$$ function utau_log_layer_only(u,y,rm) 229c$$$ real*8 u,err,utau,yrmi,y,rm,lnyp,f,dfds,rat 230c$$$ err=1.0d-6 231c$$$ utau=0.04 232c$$$ yrmi=y/rm 233c$$$ do iter=1,50 234c$$$ lnyp=log(yrmi*utau) 235c$$$ f=u-utau*(2.5*lnyp+5.2) 236c$$$ dfds=-2.5*(lnyp+6.2) 237c$$$ rat=-f/dfds 238c$$$ utau=utau+rat 239c$$$ if(utau.gt.0.5) then !flow is laminar for now give the parabolic 240c$$$ utau=sqrt(2.0*rm*u)/y 241c$$$ goto 20 242c$$$ endif 243c$$$ 244c$$$ if(abs(rat).le.err) goto 20 245c$$$ enddo 246c$$$ write(*,*)'utau failed to converge',r,dfds,rat,utau,y,u,rm 247c$$$ stop 248c$$$ 20 continue 249c$$$ return 250c$$$ end 251 252 function savarw(up,rm,cv1,nodw,x) 253 implicit none 254 real*8 err,savarw,rm,up,cv1,f,dfds,rat,efac 255 real*8 pt5,kappa,B,xmut,chi3,denom,cv1_3 256 integer iter 257 integer nodw 258 real*8 x(3) 259 pt5=0.5 260 err=1.0d-6 261 savarw=rm*cv1*1.2599 ! inflection point chi=cv1*cuberoot(2) 262 kappa=0.4 263c$$$ B=5.5 264 efac=0.1108 ! exp(-kappa*B) 265 xmut=rm*kappa*efac*(exp(kappa*up)-1.0-kappa*up 266 & -pt5*(kappa**2)*(up**2)) 267 do iter=1,50 268 chi3=savarw/rm 269 chi3=chi3*chi3*chi3 270 cv1_3=cv1**3 271 denom=chi3+cv1_3 272 273 f=savarw*chi3/denom - xmut 274 dfds=chi3*(chi3+4.0*cv1_3)/(denom**2) 275 rat=-f/dfds 276 savarw=savarw+rat 277 if(abs(rat).le.err) goto 20 278 enddo 279 write(*,*)'savarw failed to converge at ',nodw,x(1),x(2),x(3) 280 write(*,*) 'uplus, mu' 281 write(*,*) up,rm 282 write(*,*) 'dfds, rat, savarw' 283 write(*,*) dfds,rat,savarw 284 savarw=10000.0*rm 285 20 continue 286 return 287 end 288