- ( -

• (side, axis) if(boundaryCondition(side, axis).eq.noSlipWall) then is1=0 is2=0 is3=0 if(axis.eq.0) then is1=1-2 *side n1a=indexRange(side, axis)!-is1!boundary is 1 pt outside n1b=n1a else if(axis.eq.1) then is2=1-2 *side n2a=indexRange(side, axis)!-is2 n2b=n2a else is3=1-2 *side n3a=indexRange(side, axis)!-is3 n3b=n3a end if io(1)=0 io(2)=0 io(3)=0 io(axis+1)=1-2 *side ibb=indexRange(0, axis) ibe=indexRange(1, axis)-1c write(*,*) ibb, ibe do ii3=n3a, n3b do ii2=n2a, n2b do ii1=n1a, n1b if(ii3.ge.ktrip.and.ii2.ge.jtrip.and.ii1.ge.itrip) then i1=ii1 i2=ii2 i3=ii3 if(nd.eq.2) then if(axis.eq.0) then ditrip=ii2-jtrip else ditrip=ii1-itrip endif else if(axis.eq.0) then ditrip=min((ii3-ktrip),(ii2-jtrip)) else if(axis.eq.1) then ditrip=min((ii1-itrip),(ii3-ktrip)) else ditrip=min((ii1-itrip),(ii2-jtrip)) endif endif ctrans=(1-exp(-ditrip/3.))**2c ctrans=1c write(*,*) i1, i2, i3, ctrans norm(1)=0 norm(2)=0 norm(3)=0 norm(1)=rxi(axis, 0) norm(2)=rxi(axis, 1) if(nd.eq.3) norm(3)=rxi(axis, 2) nmag=sqrt(norm(1)*norm(1)+norm(2)*norm(2)+norm(3)*norm(3)) norm(1)=norm(1)/nmag norm(2)=norm(2)/nmag norm(3)=norm(3)/nmag ftan(1)=0 ftan(2)=0 ftan(3)=0 if(nd.eq.2) then if(gridType.eq.rectangular) then ftan(1)=2 *norm(1)*ux2(uc)+norm(2)*(ux2(vc)+uy2(uc)) ftan(2)=norm(1)*(uy2(uc)+ux2(vc))+2 *norm(2)*uy2(vc) else ftan(1)=2 *norm(1)*ux2c(uc)+norm(2)*(ux2c(vc)+uy2c(uc)) ftan(2)=norm(1)*(uy2c(uc)+ux2c(vc))+2 *norm(2)*uy2c(vc) end if else if(gridType.eq.rectangular) then ftan(1)=2 *norm(1)*ux2(uc)+norm(2)*(ux2(vc)+uy2(uc))+norm(3)*(ux2(wc)+uz2(uc)) ftan(2)=norm(1)*(ux2(vc)+uy2(uc))+2 *norm(2)*uy2(vc)+norm(3)*(uy2(wc)+uz2(vc)) ftan(3)=norm(1)*(ux2(wc)+uz2(uc))+norm(2)*(uy2(wc)+uz2(vc))+2 *norm(3)*uz2(wc) else ftan(1)=2 *norm(1)*ux3c(uc)+norm(2)*(ux3c(vc)+uy3c(uc))+norm(3)*(ux3c(wc)+uz3c(uc)) ftan(2)=norm(1)*(ux3c(vc)+uy3c(uc))+2 *norm(2)*uy3c(vc)+norm(3)*(uy3c(wc)+uz3c(vc)) ftan(3)=norm(1)*(ux3c(wc)+uz3c(uc))+norm(2)*(uy3c(wc)+uz3c(vc))+2 *norm(3)*uz3c(wc) end if end if fdotn=ftan(1)*norm(1)+ftan(2)*norm(2)+ftan(3)*norm(3) ftan(1)=ftan(1)-norm(1)*fdotn ftan(2)=ftan(2)-norm(2)*fdotn ftan(3)=ftan(3)-norm(3)*fdotn tauw=nu *sqrt(ftan(1)*ftan(1)+ftan(2)*ftan(2)+ftan(3)*ftan(3)) c yplus=y *yscale yscale=sqrt(tauw)/nu!assuming density=1 here...ymax=0 lmixmax=0 lmix2max=0 maxumag=0 ulmax=0 do i=ibb, ibe i1=ii1+io(1)*i i2=ii2+io(2)*i i3=ii3+io(3)*i u(i1, i2, i3, nc)=0 if(gridType.eq.rectangular) then if(nd.eq.2) then vort=abs(ux2(vc)-uy2(uc)) else vort=sqrt((uy2(wc)-uz2(vc))*(uy2(wc)-uz2(vc))-(ux2(wc)-uz2(uc))*(ux2(wc)-uz2(uc))+(ux2(vc)-uy2(uc))*(ux2(vc)-uy2(uc))) end if else if(nd.eq.2) then vort=abs(ux2c(vc)-uy2c(uc)) else vort=sqrt((uy3c(wc)-uz3c(vc))*(uy3c(wc)-uz3c(vc))-(ux3c(wc)-uz3c(uc))*(ux3c(wc)-uz3c(uc))+(ux3c(vc)-uy3c(uc))*(ux3c(vc)-uy3c(uc))) end if end if yplus=dw(i1, i2, i3)*yscale lmixw=vort *kbl *kbl *dw(i1, i2, i3)*dw(i1, i2, i3)*(1.-exp(-yplus/a0p))**2 c write(*,*)"yplus, vort ", yplus, vortc write(*,*)"dw, yscale, yplus, lmixw is ", dw(i1, i2, i3)," ", yscale," ", yplus," ", lmixw magu=u(i1, i2, i3, uc)*u(i2, i2, i3, uc)+u(i1, i2, i3, vc)*u(i1, i2, i3, vc) if(nd.eq.3) magu=magu+u(i1, i2, i3, wc)*u(i1, i2, i3, wc) magumax=max(magu, maxumag) if((vort *kbl *dw(i1, i2, i3)*(1.-exp(-yplus/a0p))).gt.lmixmax) then ymax=dw(i1, i2, i3) ulmax=magu lmixmax=vort *kbl *dw(i1, i2, i3)*(1.-exp(-yplus/a0p)) lmix2max=lmixwc write(*,*)"--", i, ymax, lmixmax, lmix2max end if u(i1, i2, i3, nc)=lmixw end do!i=ibb, ibe c now that we know lmixmax, ulmax and maxumag we can compute the eddy viscosity magumax=sqrt(magumax) ulmax=sqrt(ulmax) c write(*,*)"ymax is ", ymax," lmix2max ", lmix2max iswitch=0 do i=ibb, ibe i1=ii1+io(1)*i i2=ii2+io(2)*i i3=ii3+io(3)*i vto=alpha *ccp *min(ymax *lmixmax/kbl, cwk *ymax *(maxumag-ulmax)*(maxumag-ulmax)*kbl/lmixmax)/(1+5.5 *(dw(i1, i2, i3)*ckleb/ymax)**6) c vto=alpha *ccp *ymax *lmixmax/kbl/(1+5.5 *(dw(i1, i2, i3)*ckleb/ymax)**6) c write(*,*) ymax, dw(i1, i2, i3) c write(*,* : lineSolveBL.h