1 subroutine genscale(y, x, iBC) 2c 3c---------------------------------------------------------------------- 4c This subroutine calculate the y^+ and eta at inflow and internal face. 5c From these generate the scaling for the inflow data. 6c 7c input: 8c iBC (numnp) : boundary condition code 9c x (numnp,nsd) : node coordinates 10c 11c output: 12c y (numnp,ndof) : initial values of Y variables 13c 14c 15c Elaine Bohr december 2001 16c---------------------------------------------------------------------- 17c 18 use spebc 19c use pointer_data 20 include "common.h" 21 include "mpif.h" 22 include "auxmpi.h" 23c 24 dimension y(numnp,ndof), iBC(numnp), 25 & x(numnp,nsd), velbarR(nfint,nflow) 26 dimension ifath(numnp), velbarl(nelint,nshl,nflow), 27 & v1(nfint), ymapped(numnp,ndof), 28 & shapef(nshl), shgradl(nshl,nsd), 29 & xsi(nsd), yintl(nelint,nshl,nflow), 30 & flucl(nelint,nshl,nflow), 31 & ubarintl(nelint,nshl,nflow), 32 & fluc1(npin,nflow), fluc2(npin,nflow), 33 & ubar1(npin,nflow), ubar2(npin,nflow) 34 integer element, dir 35 36 real*8 ymax, displTi, displTr, correction 37 real*8 freestream(nflow) 38 save deltaint 39 40 41c return 42 ymapped(:,2:4)=y(:,1:3) 43 ymapped(:,1)=y(:,4) 44 ymapped(:,5)=y(:,5) 45 46 ubar2 = 0 47 fluc2 = 0 48 49 ymax = xyn(nfint) 50 51c 52c .... Localizing the solution vector on virtual plane 53c 54 55 do i = 1, nelint 56 do j = 1, 3 57 yintl(i,:,j+1) = y(ien2D(i,:),j) 58 enddo 59 yintl(i,:,1) = y(ien2D(i,:),4) 60 if(nflow.gt.4) then 61 do j = 5, nflow 62 yintl(i,:,j) = y(ien2D(i,:),j) 63 enddo 64 endif 65 enddo 66 67c 68c .... Finding averaged velocity in spanwise direction 69c for the virtual plane 70c 71 72 do i=1,nfint 73 velbarR(i,:)=0 74 do j=1,imax(i)+1 75 call shptet(ipord,xsinfin(i,j,:),shapef(:),shgradl(:,:)) 76 do k=1,nshl 77 velbarR(i,:)=velbarR(i,:) 78 & + yintl(elcnfin(i,j),k,:)*shapef(k) 79 enddo 80 enddo 81 velbarR(i,:)=velbarR(i,:) / (imax(i)+1) 82 enddo 83 84c 85c .... Label the nodes that near the BL thickness 86c 87 88 if (thetag.eq.0.0) then 89 dir = 2 90 else 91 dir = 4 92 endif 93 94 v1(1)=10.0 95 do i=2,nfint+1 96 v1(i)=velbarR(i-1,dir)-0.99*vel 97 if((v1(i).gt.0).and.(v1(i-1).le.0)) then 98 label=i-1 99 go to 200 100 endif 101 enddo 102 label=i-1 103 104c 105c.... Find the BL thickness by means of finding the y coord 106c 107 108 200 continue 109 dv=velbarR(label,dir)-velbarR(label-1,dir) 110 dy=xyn(label)-xyn(label-1) 111 112c 113c .... Current calculation of bl thickness at recycle plane 114c 115 116 if(istep.ne.0) then 117 dlast=deltaint 118 deltaint=xyn(label-1) 119 & + dy*(0.99*vel-velbarR(label-1,dir))/dv 120 121c 122c .... Early transients cause jumpy delta, smooth it. 123c 124 125 deltaint=min(1.05*dlast,max(deltaint,0.95*dlast)) 126 else 127 deltaint=xyn(label-1) 128 & + dy*(0.99*vel-velbarR(label-1,dir))/dv 129 endif 130 131c 132c .... Deltaint is now the ratio of BL thickness at the interior plane 133c to the BL thickness at the inlet plane 134c 135 136 deltaint=min(two*rbltin,max(deltaint,pt5*rbltin)) 137 rdelta = deltaint/rbltin 138 139c 140c .... Finding freestream solutions 141c 142 143 freestream = 0 144 icount = 0 145 do i=1, nfint 146 if (xyn(i).ge.deltaint) then 147 freestream(:) = freestream(:) + velbarR(i,:) 148 icount = icount + 1 149 endif 150 enddo 151 freestream = freestream / icount 152 153c 154c .... Putting the freestream values into the average outside the BLT 155c 156 157 do i=1, nfint 158 if (xyn(i).ge.deltaint) then 159 velbarR(i,:) = freestream(:) 160 endif 161 enddo 162 163c 164c .... Localizing the averaged velocity found above 165c 166 167 do i=1,nelint 168 do k=1,nshl 169 do j=1,nfint-1 170 if (thetag.eq.0.0) then 171 if ((x(ien2D(i,k),2).ge.xyn(j)) .and. 172 & (x(ien2D(i,k),2).le.(xyn(j+1)+0.000001))) then 173 tmp = (x(ien2D(i,k),2) - xyn(j)) / 174 & (xyn(j+1) - xyn(j)) 175 do l=1,nflow 176 velbarl(i,k,l) = 177 & (velbarR(j+1,l) - velbarR(j,l)) * 178 & tmp + velbarR(j,l) 179 enddo 180 endif 181 else 182 if ((xcyl(ien2D(i,k),1).ge.xcyl(nrint(j+1),1)) .and. 183 & (xcyl(ien2D(i,k),1).le.xcyl(nrint(j),1))) then 184 tmp = (xcyl(ien2D(i,k),1) - xcyl(nrint(j+1),1)) / 185 & (xcyl(nrint(j),1) - xcyl(nrint(j+1),1)) 186 do l=1,nflow 187 velbarl(i,k,l) = 188 & (velbarR(j,l) - velbarR(j+1,l)) * 189 & tmp + velbarR(j+1,l) 190 enddo 191 endif 192 endif 193 enddo 194 enddo 195 enddo 196 197c 198c --- For now only Blasius is coded --- 199c 200 201c 202c .... Calculate fluctuations on elements of internal plane 203c 204 205c flucl = yintl - velbarl 206 207c 208c .... Calculate mean values on elements of internal plane 209c 210 211 ubarintl = velbarl 212 213c 214c .... Calculating the coordinates of the point from where the 215c solution will be projected to the inlet plane 216c 217 218 219 do i=1,npin 220 221c 222c .... Cartesian coodinate system 223c 224 225 if (thetag.eq.0.0) then 226 xts1 = x(nen1(i),1) + plandist 227 if (xynin(i)*rdelta.gt.ymax) then 228 xts2 = ymax 229 else 230 xts2 = xynin(i)*rdelta 231 endif 232 xts3 = x(nen1(i),3) 233 234c 235c .... Cylindrical coordinate system 236c 237 238 else 239 if (xynin(i).le.0.00001) then 240 xts1 = (radcyl-xynin(i)*rdelta*sang-tolerence) 241 & *cos(xcyl(nen1(i),2)) 242 xts2 = (radcyl-xynin(i)*rdelta*sang-tolerence) 243 & *sin(xcyl(nen1(i),2)) 244 xts3 = (aR-(radcyl-xynin(i)*rdelta*sang-tolerence) 245 & * (xnrml*cos(xcyl(nen1(i),2)) 246 & + ynrml*sin(xcyl(nen1(i),2))))/znrml 247 elseif (xynin(i)*rdelta.gt.ymax) then 248 xts1 = (radcyl-ymax*sang) 249 & *cos(xcyl(nen1(i),2)) 250 xts2 = (radcyl-ymax*sang) 251 & *sin(xcyl(nen1(i),2)) 252 xts3 = (aR-(radcyl-ymax*sang) 253 & * (xnrml*cos(xcyl(nen1(i),2)) 254 & + ynrml*sin(xcyl(nen1(i),2))))/znrml 255 else 256 xts1 = (radcyl-xynin(i)*rdelta*sang) 257 & *cos(xcyl(nen1(i),2)) 258 xts2 = (radcyl-xynin(i)*rdelta*sang) 259 & *sin(xcyl(nen1(i),2)) 260 xts3 = (aR-(radcyl-xynin(i)*rdelta*sang) 261 & * (xnrml*cos(xcyl(nen1(i),2)) 262 & + ynrml*sin(xcyl(nen1(i),2))))/znrml 263 endif 264 endif 265 266c 267c .... Searching for the appropriate element 268c 269 270 call elem_search(xintl, xts1, xts2, xts3, 271 & xsi(:), element, 2) 272 call shptet(ipord,xsi(:),shapef(:),shgradl(:,:)) 273 274c 275c .... Calculating the average velocity and fluctuations 276c for the inlet plane 277c 278 279 do k=1,nshl 280 fluc2(i,:)= 0 !fluc2(i,:) + flucl(element,k,:)*shapef(k) 281 ubar2(i,:)=ubar2(i,:) + ubarintl(element,k,:)*shapef(k) 282 enddo 283 enddo 284 285 286c$$$c 287c$$$c keep freestream values set through averages 288c$$$c 289c$$$ ubaro=0 290c$$$ tbaro=0 291c$$$ icount=0 292c$$$ do i=1,nfin 293c$$$ if(yin(i).ge.rbltin) then 294c$$$ nzl=nsons(i) !Elaine 295c$$$ nzb=ienson1(i,1) 296c$$$ nze=nzb+nzl-1 297c$$$ tbaro=tbaro+2.0*ubar2(i,5)+sum(fluc2(nzb:nze,5)) 298c$$$ ubaro=ubaro +sum(fluc2(nzb:nze,2)) 299c$$$ icount=icount+nzl 300c$$$ endif 301c$$$ enddo 302c$$$ 303c$$$c alternative to myway 304c$$$c 305c$$$ ubaro=ubaro/icount 306c$$$ tmeaninflow=0.0625097048890964 307c$$$ fact= tmeaninflow/(tbaro/icount) 308c$$$ if (fact.ge. 0.9999999 .and. fact.le.1.0000001) fact = 1.0 309c$$$ 310c$$$ do i=1,nfin 311c$$$ if(yin(i).ge.rbltin) then 312c$$$ ubar2(i,2)=1.0-ubaro 313c$$$ endif 314c$$$ enddo 315 fact = 1.0 316 rvscal = 1.0 317 318c 319c .... Putting the freestream value outside the BLT into ubar2 320c 321 322 do i = 1, npin 323 if (xynin(i).ge.rbltin) then 324 ubar2(i,:) = freestream(:) 325c ubar2(i,dir) = 1.0 326 fluc2(i,:) = 0 327 endif 328 enddo 329 330c$$$c 331c$$$c .... For the cylindrical case the freestream velocity needs 332c$$$c to be corrected for the blockage 333c$$$c 334c$$$ 335c$$$ if (thetag.ne.0.0) then 336c$$$ displTi = 0.0 337c$$$ displTr = 0.0 338c$$$ do i = 2, nfint 339c$$$ 340c$$$c 341c$$$c .... Displacement thickness for inlet plane 342c$$$c 343c$$$ 344c$$$ displTi = displTi + (1 - y(nrint(i),3)) 345c$$$ & * (xyn(i) - xyn(i-1)) * (radcyl - xyn(i)) 346c$$$ 347c$$$c 348c$$$c .... Displacement thickness for recycle plane 349c$$$c 350c$$$ 351c$$$ displTr = displTr + (1 - velbarR(i,4)) 352c$$$ & * (xyn(i) - xyn(i-1)) * (radcyl - xyn(i)) 353c$$$ enddo 354c$$$c displTi = radcyl - sqrt(radcyl*radcyl - displTi) 355c$$$c displTr = radcyl - sqrt(radcyl*radcyl - displTr) 356c$$$ correction = (radcyl*radcyl - displTr) 357c$$$ & / (radcyl*radcyl - displTi) 358c$$$ else 359 correction = 1.0 360c$$$ endif 361 362c 363c .... Scaled plane extraction boundary condition 364c 365 366 ymapped(nen1(1:npin),1)= correction * (ubar2(:,1)+fluc2(:,1)) 367 ymapped(nen1(1:npin),2)= correction * 368 & (ubar2(:,2)+fluc2(:,2)) !myway *factu 369 ymapped(nen1(1:npin),3)= correction * 370 & (ubar2(:,3)+fluc2(:,3))*rvscal 371 ymapped(nen1(1:npin),4)= correction * (ubar2(:,4)+fluc2(:,4)) 372c & freestream(4) 373 ymapped(nen1(1:npin),5)= correction * fact*(ubar2(:,5)+fluc2(:,5)) 374 375 376 377c 378c .... Ready to put the solution on the inflow plane 379c 380 381 if(intpres.eq.1) then !interpolating pressure at inflow 382 where (btest(iBC,11)) 383 y(:,1) = ymapped(:,2) 384 y(:,2) = ymapped(:,3) 385 y(:,3) = ymapped(:,4) 386 y(:,4) = ymapped(:,1) 387 y(:,5) = ymapped(:,5) 388 endwhere 389 else ! not interpolating pressure at inflow 390 where (btest(iBC,11)) 391 y(:,1) = ymapped(:,2) 392 y(:,2) = ymapped(:,3) 393 y(:,3) = ymapped(:,4) 394 y(:,5) = ymapped(:,5) 395 endwhere 396 endif 397 398c 399c debugging variables 400c 401c if(iter.eq.nitr) then 402c write(555,556)lstep+1,deltaint,label,nfint 403c write(554,556)lstep+1,yplusi(2),ypluso(2),factu,factt,gamt 404c call flush(554) 405c call flush(555) 406c endif 407c 556 format(i6,5(2x,e14.7)) 408c 409c.... return 410c 411 412 return 413 end 414