program channel c c implement newton method to solve for interactive channel flow c implicit real*8 (a-h,o-z) implicit integer (i-n) parameter (iy=64,iyp=iy+1,iz=3*(iy+1)+1,ixm=128,ixp=128,iyg=128) real *8 u(-ixm:ixp,0:iy),du(0:iyp) real *8 w(-ixm:ixp,0:iy),dw(0:iyp) real *8 psi(-ixm:ixp,0:iy),dpsi(0:iyp) real *8 v(-ixm:ixp,0:iy) real *8 p(-ixm:ixp),px(-ixm:ixp),dpx real *8 a(-ixm:ixp),ad(-ixm:ixp) real *8 cof(iz,iz),rhs(iz) real *8 x(-ixm:ixp),hx(-ixm:ixp),hy(0:iy),y(0:iy) real *8 xs(-ixm:ixp,0:iy),ys(-ixm:ixp,0:iy) real *8 dxi,xm,z,zm,third,hxi real *8 hs,h2,alpha,fac,fac1,xt,rho,nu,var real *8 tolerance,hx0 c c need matrices to hold physical coordinates and velocities c real *8 yg(-ixm:ixp,0:iyg),dyg(-ixm:ixp) real *8 yl(-ixm:ixp,0:iy),yu(-ixm:ixp,0:iy) real *8 uo(-ixm:ixp,0:iyg),psio(-ixm:ixp,0:iyg) real *8 ucomp(-ixm:ixp,0:iyg),uinner(-ixm:ixp,0:iy),uloc(0:iy) real *8 upper(-ixm:ixp,0:iy),yloc(0:iy),flux(-ixm:ixp) 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 ixu,ixd,ipx,ipy,iy,i,j,k,msize,msizep,row c c specify flux through channel c fluxin=1. fluxcorrection=fluxin*12 fluxcorrection=1. c set up initial guess tolerance=1.e-6 ixu=ixm ixd=ixp ipx=1 ipy=1 msize=ixm msizep=ixp iterations=200 c c set up mesh c hx0=1./(ixu+ixd) hy0=0.05 facy=1.04 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) facu=1.0 facd=1.0 x(0)=0. x(1)=hx0 hx(1)=hx0 hx(0)=hx0 x(-1)=-hx0 do i=2,msizep hx(i)=facd*hx(i-1) x(i)=x(i-1)+hx(i) enddo do i=2,msize hx(-i+1)=facu*hx(-i+2) x(-i)=x(-i+1)-hx(-i+1) enddo xm=x(ixd) write(7,7)(x(i),hx(i),i=-ixu,ixd) write(8,7)(y(j),hy(j),j=0,iy) 7 format(2e15.5) write(6,5)ixu,ixd,iy,zm,xm,hx0,hy0,msize,msizep 5 format('Interactive channel solver',/, 'ixu = ',i4, + ' ixd = ',i4, ' iy = ',i4, + /,'zm =',f8.0, 'xm = ', f12.4,/, + ' hx0 = ',f10.4, ' hy0 = ',f10.4,2i10) write(6,10)k,var,x(-ixu/2),x(0),x(ixd/4), + x(ixd/2),x(3*ixd/4),x(ixd),p(0) c set up initial guess for a(x),p(x) write(6,1) 1 format('Initialising a(x)') c c++++++++++++++++++++++++++++++++++ c + write(6,*)' Input wall size factor' read(5,*)sfac c + c++++++++++++++++++++++++++++++++++ c do i=-ixu,ixd t=(x(i)-x(-ixu))/(x(ixd)-x(-ixu)) do j=0,iy xs(i,j)=x(i) ys(i,j)=-sfac*(1-cos(2*3.14159*t))/2+y(j) enddo c c note that a(X) here is B(X)=A(X)+f(X)=f/2=-A if g=0 c a(i)=ys(i,0)/2 p(i)=0 px(i)=0 end do write(6,3) 3 format('Setting up plate start solution') do i=-ixu,ixp u(i,0)=0. psi(i,0)=0. w(i,0)=fluxcorrection/2 do j=1,iy z=y(j) psi(i,j)=fluxcorrection*(z*z/2+z*a(i))/2 u(i,j)=fluxcorrection*(z+a(i))/2 w(i,j)=fluxcorrection/2 end do end do c now have initial guess for all variables, c output intial values c write(9,21) ((xs(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(10,21)((ys(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(13,21)((psi(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(14,21)((u(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(15,65)(x(i),a(i), i=-ixu,ixd) write(16,65)(x(i),p(i), i=-ixu,ixd) write(17,65)(x(i),ad(i), i=-ixu,ixd) write(18,65)(x(i)-0.5*hx(i),px(i),i=-ixu+1,ixd) write(19,66)(y(j),u(-ixu,j),j=0,iy) write(19,66)(y(j),u(0,j),j=0,iy) write(19,66)(y(j),u(ixd/4,j),j=0,iy) write(19,66)(y(j),u(ixd/2,j),j=0,iy) write(19,66)(y(j),u(3*ixd/4,j),j=0,iy) write(19,66)(y(j),u(ixd,j),j=0,iy) write(20,67)(x(i),w(i,0),i=-ixu,ixd) 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=-ixu+1,ixp c-------------------------------------------------------------------- c iterate along plate rows - u(0)=0 | c-------------------------------------------------------------------- do newton=1,5 do l=1,iz do m=1,iz cof(l,m)=0. enddo enddo 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*hy(j) rhs(row)=2*hy(j)*(px(i) + + uc*ux - wc*phix) -wy enddo row=row+1 cof(row,row-1)=1. rhs(row)=0 row=row+1 cof(row,row-1)=1. rhs(row)=0 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(ixd,iy) + +abs(du(j))/u(ixd,iy) + +abs(dw(j))/w(ixd,iy) enddo c write(6,72)error/(3*iy) 72 format('plate',e12.3) dpx=rhs(iz) px(i)=px(i)+dpx enddo enddo do i=-ixu+1,ixd p(i)=p(i-1)+hx(i)*px(i) enddo c c ---------------- impose periodicity c do j=0,iy psi(-ixu,j)=psi(ixd,j) u(-ixu,j)=u(ixd,j) w(-ixu,j)=w(ixd,j) enddo wmin=100. do i=-ixu,ixd if(w(i,0).lt.wmin)wmin=w(i,0) enddo write(6,10)k,wmin,w(-ixu/2,0),w(0,0),w(ixd/4,0), + w(ixd/2,0),w(3*ixd/4,0),w(ixd,0),p(0),px(0) 10 format(i6,(10e12.3)) c----------------------------------------------------------------------------------- c end of iteration loop c----------------------------------------------------------------------------------- enddo write(22,21)((psi(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(23,21)((u(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(24,21)((w(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(25,65)(x(i), a(i), i=-ixu,ixd) write(26,65)(x(i),p(i), i=-ixu,ixd) write(27,65)(x(i), ad(i), i=-ixu,ixd) write(28,65)(x(i)-0.5*hx(i),px(i),i=-msize+1,msizep) write(29,66)(y(j),u(-ixu,j),j=0,iy) write(29,66)(y(j),u(0,j),j=0,iy) write(29,66)(y(j),u(ixd/4,j),j=0,iy) write(29,66)(y(j),u(ixd/2,j),j=0,iy) write(29,66)(y(j),u(3*ixd/4,j),j=0,iy) write(29,66)(y(j),u(ixd,j),j=0,iy) write(30,67)(x(i),w(i,0),i=-ixu,ixd) 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=-ixu,ixd 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=-ixu,ixd 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=-ixu,ixd) write(36,66)(xs(i,0)/epsilon,yg(i,iyg),i=-ixu,ixd) do i=-ixu,ixd,ixu/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=-ixu,ixd) 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 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 c**************************************************** c YOU MUST PROVIDE THIS ROUTINE **** c for instance from Numerical Recipies **** c**************************************************** c return end