program freeint c c free interaction in a parallel channel c c have to set the initial value of a c c internal channel problem c c implicit real*8 (a-h,o-z) implicit integer (i-n) parameter (iy=50,iyp=iy+1,iz=2*3*(iy+1)+3,ix=100,iyg=128) c c variables on upper wall c real *8 uu(0:ix,0:iy),duu(0:iyp) real *8 wu(0:ix,0:iy),dwu(0:iyp) real *8 psiu(0:ix,0:iy),dpsiu(0:iyp) real *8 vu(0:ix,0:iy) real *8 pu(0:ix),pxu(0:ix),dpxu c c variables on lower wall c real *8 u(0:ix,0:iy),du(0:iyp) real *8 w(0:ix,0:iy),dw(0:iyp) real *8 psi(0:ix,0:iy),dpsi(0:iyp) real *8 v(0:ix,0:iy) real *8 p(0:ix),px(0:ix),dpx real *8 f(0:ix),g(0:ix) real *8 a(-3:ix+2),ad(0:ix) real *8 cof(iz,iz),rhs(iz) real *8 x(0:ix),hx(0:ix),hy(0:iy),y(0:iy) real *8 xs(0:ix,0:iy),ys(0:ix,0:iy) real *8 dxi,xm,z,zm,third,hxi real *8 hs,h2,alpha,fac,fac1,xt,rho,nu,var real *8 mu0,gamma0,lambda0,zeta0,tolerance,hx0 c c need matrices to hold physical coordinates and velocities c real *8 yg(0:ix,0:iyg),dyg(0:ix) real *8 yl(0:ix,0:iy),yu(0:ix,0:iy) real *8 uo(0:ix,0:iyg),psio(0:ix,0:iyg) real *8 ucomp(0:ix,0:iyg),uinner(0:ix,0:iy),uloc(0:iy) real *8 upper(0:ix,0:iy),yloc(0:iy),flux(0:ix) c c yg= y coordinates in global domain c yl= y coordinates in lower wall layer } c yu= y coordinates in upper wall layer } in global coordinates, only known when R given c ys= Y coordinates c integer ix,ipx,ipy,iy,i,j,k,msize,msizep,row c set up initial guess c=5.73 tolerance=1.e-6 ipx=1 ipy=1 iterations=1 c c set up mesh c hy0=.01 hx0=0.013 facy=1.1 y(0)=0. hy(1)=hy0 y(1)=y(0)+hy0 do j=2,iy hy(j)=facy*hy(j-1) y(j)=y(j-1)+hy(j) enddo zm=y(iy) do i=0,ix hx(i)=hx0 x(i)=i*hx0 enddo xm=x(ix) write(7,7)(x(i),hx(i),i=0,ix) write(8,7)(y(j),hy(j),j=0,iy) 7 format(2e15.5) write(6,5)ix,iy,zm,xm,hx0,hy0,msize,msizep 5 format('Interactive channel solver',/, 'ix = ',i4,' iy = ',i4, + /,'zm =',f8.0, 'xm = ', f12.4,/, + ' hx0 = ',f10.4, ' hy0 = ',f10.4,2i10) write(6,10)k,var,x(0),x(ix/4), + x(ix/2),x(3*ix/4),x(ix),p(0) c set up initial guess for a(x),p(x) write(6,1) 1 format('Initialising a(x)') c do i=0,ix do j=0,iy xs(i,j)=x(i) ys(i,j)=y(j) enddo a(i)=0. p(i)=0 px(i)=0 pu(i)=0 pxu(i)=0 end do a(0)=-0.0001 a(-1)=a(0)*exp(-c*hx0) a(-2)=a(-1)*exp(-c*hx0) a(-3)=a(-2)*exp(-c*hx0) a(ix+1)=0. a(ix+2)=0. write(6,3) 3 format('Setting up plate start solution') do i=0,ix u(i,0)=0. psi(i,0)=0. w(i,0)=1./2 uu(i,0)=0. psiu(i,0)=0. wu(i,0)=1./2 do j=1,iy z=y(j) psi(i,j)=z*z/4 u(i,j)=z/2 w(i,j)=1./2 psiu(i,j)=z*z/4 uu(i,j)=z/2 wu(i,j)=1./2 end do end do c c set pressure perturbation c p(0)=0.0 c now have initial guess for all variables, c output intial values c write(9,21) ((xs(i,j),j=0,iy,ipy),i=0,ix,ipx) write(10,21)((ys(i,j),j=0,iy,ipy),i=0,ix,ipx) write(13,21)((psi(i,j),j=0,iy,ipy),i=0,ix,ipx) write(14,21)((u(i,j),j=0,iy,ipy),i=0,ix,ipx) write(15,65)(x(i),a(i), i=0,ix) write(16,65)(x(i),p(i), i=0,ix) write(17,65)(x(i),ad(i), i=0,ix) write(18,65)(x(i)-0.5*hx(i),px(i),i=0+1,ix) write(19,66)(y(j),u(0,j),j=0,iy) write(19,66)(y(j),u(0,j),j=0,iy) write(19,66)(y(j),u(ix/4,j),j=0,iy) write(19,66)(y(j),u(ix/2,j),j=0,iy) write(19,66)(y(j),u(3*ix/4,j),j=0,iy) write(19,66)(y(j),u(ix,j),j=0,iy) write(20,67)(x(i),w(i,0),i=0,ix) 67 format('vdata',/,(f12.3,f15.6)) 66 format('udata',/,(f12.3,f15.6)) 65 format('data',/, (f12.3,f15.6)) c-------------------------------------------------------------------- c main iteration | c-------------------------------------------------------------------- do k=1,iterations do i=1,ix do j=1,iy psi(i,j)=psi(i-1,j) u(i,j)=u(i-1,j) w(i,j)=w(i-1,j) psiu(i,j)=psiu(i-1,j) uu(i,j)=uu(i-1,j) wu(i,j)=wu(i-1,j) enddo c-------------------------------------------------------------------- c iterate along channel rows - u(0)=0 | c-------------------------------------------------------------------- do newton=1,5 do l=1,iz do m=1,iz cof(l,m)=0. enddo enddo c c lower layer c row=1 c set first two rows cof(row,row)=1 rhs(row)=0 row=row+1 cof(row,row)=1 rhs(row)=0 c loop through iy sets do j=1,iy wc=(w(i,j)+w(i-1,j)+ + w(i,j-1)+w(i-1,j-1))/4 uc=(u(i,j)+u(i-1,j)+ + u(i,j-1)+u(i-1,j-1))/4 phiy=psi(i,j)+ psi(i-1,j)- + psi(i,j-1)-psi(i-1,j-1) uy= u(i,j)+ u(i-1,j)- + u(i,j-1)- u(i-1,j-1) wy= w(i,j)+ w(i-1,j)- + w(i,j-1)- w(i-1,j-1) phix=(psi(i, j)+psi(i, j-1)- + psi(i-1,j)-psi(i-1,j-1))/(2*hx(i)) ux=( u(i, j)+ u(i, j-1)- + u(i-1,j)- u(i-1,j-1))/(2*hx(i)) row=row+1 cof(row,row-2)=-1 cof(row,row-1)=-2*hy(j)/4 cof(row,row+1)=1 cof(row,row+2)=-2*hy(j)/4 rhs(row)=2*hy(j)*uc-phiy row=row+1 cof(row,row-2)=-1 cof(row,row-1)=-2*hy(j)/4 cof(row,row+1)=1 cof(row,row+2)=-2*hy(j)/4 rhs(row)=2*hy(j)*wc-uy row=row+1 cof(row,row-4)= 2*hy(j)*wc/(2*hx(i)) cof(row,row-3)= -2*hy(j)*ux/4-2*hy(j)*uc/(2*hx(i)) cof(row,row-2)=-1+2*hy(j)*phix/4 cof(row,row-1)= 2*hy(j)*wc/(2*hx(i)) cof(row,row) = -2*hy(j)*ux/4-2*hy(j)*uc/(2*hx(i)) cof(row,row+1)= 1+2*hy(j)*phix/4 cof(row,iz-2) = -2*hy(j) rhs(row)=2*hy(j)*(px(i) + + uc*ux - wc*phix) -wy enddo row=row+1 cof(row,row-1)=-1. cof(row,iz)=0.5 rhs(row)=u(i,iy)-(y(iy)+a(i))/2 row=row+1 cof(row,row-1)=1. rhs(row)=0 c c now do upper layer c row=row+1 c set first two rows cof(row,row-1)=1 rhs(row)=0 row=row+1 cof(row,row-1)=1 rhs(row)=0 c loop through iy sets do j=1,iy wc=(wu(i,j)+wu(i-1,j)+ + wu(i,j-1)+wu(i-1,j-1))/4 uc=(uu(i,j)+uu(i-1,j)+ + uu(i,j-1)+uu(i-1,j-1))/4 phiy=psiu(i,j)+ psiu(i-1,j)- + psiu(i,j-1)-psiu(i-1,j-1) uy= uu(i,j)+ uu(i-1,j)- + uu(i,j-1)- uu(i-1,j-1) wy= wu(i,j)+ wu(i-1,j)- + wu(i,j-1)- wu(i-1,j-1) phix=(psiu(i, j)+psiu(i, j-1)- + psiu(i-1,j)-psiu(i-1,j-1))/(2*hx(i)) ux=( uu(i, j)+ uu(i, j-1)- + uu(i-1,j)- uu(i-1,j-1))/(2*hx(i)) row=row+1 cof(row,row-3)=-1 cof(row,row-2)=-2*hy(j)/4 cof(row,row)=1 cof(row,row+1)=-2*hy(j)/4 rhs(row)=2*hy(j)*uc-phiy row=row+1 cof(row,row-3)=-1 cof(row,row-2)=-2*hy(j)/4 cof(row,row)=1 cof(row,row+1)=-2*hy(j)/4 rhs(row)=2*hy(j)*wc-uy row=row+1 cof(row,row-5)= 2*hy(j)*wc/(2*hx(i)) cof(row,row-4)= -2*hy(j)*ux/4-2*hy(j)*uc/(2*hx(i)) cof(row,row-3)=-1+2*hy(j)*phix/4 cof(row,row-2)= 2*hy(j)*wc/(2*hx(i)) cof(row,row-1) = -2*hy(j)*ux/4-2*hy(j)*uc/(2*hx(i)) cof(row,row)= 1+2*hy(j)*phix/4 cof(row,iz-1) = -2*hy(j) rhs(row)=2*hy(j)*(pxu(i)+ uc*ux - wc*phix) -wy enddo row=row+1 cof(row,row-2)=-1. cof(row,iz)=-0.5 rhs(row)=uu(i,iy)-(y(iy)-a(i))/2 row=row+1 cof(row,row-2)=1. rhs(row)=0 c c pressure relationship: row should be equal to iz c addd=(a(i)-3*a(i-1)+3*a(i-2)-a(i-3))/(hx(i)**3) row=row+1 cof(row,row-2)=1. cof(row,row-1)=-1. cof(row,row)= 1./(120.*hx(i)**3) rhs(row)=pxu(i)-px(i)-addd/120 c c cof set up c call invert(cof,iz,iz,rhs) error=0 do j=0,iy dpsi(j)=rhs(3*j+1) psi(i,j)=psi(i,j)+dpsi(j) du(j)=rhs(3*j+2) u(i,j)=u(i,j)+du(j) dw(j)=rhs(3*j+3) w(i,j)=w(i,j)+dw(j) error=error+abs(dpsi(j))/psi(ix,iy) + +abs(du(j))/u(ix,iy) + +abs(dw(j))/w(ix,iy) enddo do j=0,iy key=3*(iy+1)+3*j dpsiu(j)=rhs(key+1) psiu(i,j)=psiu(i,j)+dpsiu(j) duu(j)=rhs(key+2) uu(i,j)=uu(i,j)+duu(j) dwu(j)=rhs(key+3) wu(i,j)=wu(i,j)+dwu(j) error=error+abs(dpsiu(j))/psi(ix,iy) + +abs(duu(j))/u(ix,iy) + +abs(dwu(j))/w(ix,iy) enddo c write(6,72)error/(3*iy) 72 format('plate',e12.3) dpx=rhs(iz-2) px(i)=px(i)+dpx dpxu=rhs(iz-1) pxu(i)=pxu(i)+dpxu a(i)=a(i)+rhs(iz) enddo enddo c c ---------------- impose periodicity c a(ix+1)=a(ix) a(ix+2)=a(ix) c c -------------- estimate a'' at left boundary c i=0 add=(a(i+1)-2*a(i)+a(i-1))/(hx(i+1)**2) c c calculate pressure on upper wall relative to lower wall c pu(0)=p(0)+add/120 do i=1,ix p(i)=p(i-1)+hx(i)*px(i) pu(i)=pu(i-1)+hx(i)*pxu(i) enddo suma=0. do i=0,ix suma=suma+a(i)*a(i) enddo wlmax=w(0,0) wlmin=w(0,0) wumax=wu(0,0) wumin=wu(0,0) do i=1,ix if(w(i,0).gt.wlmax)wlmax=w(i,0) if(w(i,0).lt.wlmin)wlmin=w(i,0) if(wu(i,0).gt.wumax)wumax=wu(i,0) if(wu(i,0).lt.wumin)wumin=wu(i,0) enddo write(6,10)k,wlmin,wlmax,wumin,wumax,p(0),pu(0),a(0),suma 10 format(i6,(10e12.3)) c----------------------------------------------------------------------------------- c end of iteration loop c----------------------------------------------------------------------------------- enddo write(22,21)((psi(i,j),j=0,iy,ipy),i=0,ix,ipx) write(23,21)((u(i,j),j=0,iy,ipy),i=0,ix,ipx) write(24,21)((w(i,j),j=0,iy,ipy),i=0,ix,ipx) write(25,65)(x(i), a(i), i=0,ix) write(26,65)(x(i),p(i), i=0,ix) write(27,65)(x(i), ad(i), i=0,ix) write(28,65)(x(i)-0.5*hx(i),px(i),i=0,ix) write(29,66)(y(j),u(0,j),j=0,iy) write(29,66)(y(j),u(0,j),j=0,iy) write(29,66)(y(j),u(ix/4,j),j=0,iy) write(29,66)(y(j),u(ix/2,j),j=0,iy) write(29,66)(y(j),u(3*ix/4,j),j=0,iy) write(29,66)(y(j),u(ix,j),j=0,iy) write(30,67)(x(i),w(i,0),i=0,ix) write(50,67)(x(i),wu(i,0),i=0,ix) write(46,67)(x(i),pu(i),i=0,ix) write(42,21)((psiu(i,j),j=0,iy,ipy),i=0,ix,ipx) write(43,21)((uu(i,j),j=0,iy,ipy),i=0,ix,ipx) write(44,21)((wu(i,j),j=0,iy,ipy),i=0,ix,ipx) 21 format(8e15.4) do kk=1,10 c c loop up to 10 times c write(6,45) 45 format('looping various composite expansions') write(6,46) 46 format('input length of furrows') read(5,*)flength if(flength.lt.0.)go to 101 epsilon=1./flength write(6,47) 47 format('input depth of furrow') read(5,*)fdepth fac=fdepth/(2*sfac) ain=log(fac)/log(epsilon) eps=epsilon**ain write(6,48)epsilon,ain,eps 48 format('epsilon ',f10.4, ' ain ',f15.4,' eps ',f15.4) re=flength**(1+3*ain) write(6,49)re 49 format('Reynolds number = ',f15.3) c c generate composite solution c write(6,50) 50 format('Generating composite expansion') eps2=eps*eps c c doing case of flat top wall, sinuous lower wall c do i=0,ix dyg(i)=(1.-eps*ys(i,0))/iyg do j=0,iyg yg(i,j)=eps*ys(i,0)+j*dyg(i) if(yg(i,j).ge. 0.)then vo=yg(i,j)*(1.-yg(i,j))/2 dvo=0.5-yg(i,j) c A(X)=-a(X) uo(i,j)=fluxcorrection*(vo-eps*a(i)*dvo) else uo(i,j)=fluxcorrection*(yg(i,j)-eps*a(i))/2 endif enddo do j=0,iy yl(i,j)=eps*ys(i,j) yu(i,j)=1.-eps*y(j) zu=1-yu(i,j) uinner(i,j)=eps*u(i,j)- + fluxcorrection*(yl(i,j)/2-eps*a(i)/2) upper(i,j)=eps*u(i,j)- + fluxcorrection*(zu/2+eps*a(i)/2) enddo enddo write(35,66)(yg(0,j),uo(0,j),j=0,iyg) write(35,66)(yl(0,j),eps*u(0,j),j=0,iy) write(35,66)(yu(0,j),eps*u(0,j),j=0,iy) do i=0,ix do j=0,iy yloc(j)=yl(i,j) uloc(j)=uinner(i,j) enddo do j=0,iyg ucomp(i,j)=uo(i,j)+ainterp(yloc,uloc,iy,yg(i,j)) enddo do j=0,iy yloc(j)=yu(i,j) uloc(j)=upper(i,j) enddo do j=0,iyg ucomp(i,j)=ucomp(i,j)+ainterp(yloc,uloc,iy,yg(i,j)) enddo flux(i)=0. do j=1,iyg flux(i)=flux(i)+(ucomp(i,j)+ucomp(i,j-1))*(yg(i,j)-yg(i,j-1))/2 enddo enddo write(36,66)(xs(i,0)/epsilon,yg(i,0),i=0,ix) write(36,66)(xs(i,0)/epsilon,yg(i,iyg),i=0,ix) do i=0,ix,ix/2 write(36,66)(xs(i,0)/epsilon+5*ucomp(i,j),yg(i,j),j=0,iyg) enddo c write(6,88)(xs(i,0)/epsilon,flux(i),i=0,ix) 88 format(' flux = ',(2f12.5)) write(37,66)(yg(0,j),ucomp(0,j),j=0,iyg) enddo 101 continue stop end real *8 function ainterp(y,f,iy,p) real *8 y(0:iy),f(0:iy),p,fac integer j,iy c c interpolate data given at (y,f) at point p c ainterp=0. do j=1,iy if((y(j)-p)*(y(j-1)-p).le.0.)then fac=(p-y(j-1))/(y(j)-y(j-1)) ainterp=(1-fac)*f(j-1)+fac*f(j) endif enddo return end c----------------------------------------------------------------------------------- c matrix inverter | c----------------------------------------------------------------------------------- subroutine invert(a,n,np,b) implicit real *8(a-h,o-z) real*8 a(np,np),b(np) c c routine to solve ax=b needed here (E.G. from Numerical Recipes) c c maximum dimension of a is np x np c c part being used is n x n c c solution returned in b c return end