program deck2 c c Triple deck solver for c internal channel problem c solves for two layers simultaneously c implicit real*8 (a-h,o-z) implicit integer (i-n) parameter (iy=32,iyp=iy+1,iz=2*3*(iy+1)+3,ixm=40,ixp=40,iyg=128) c c iy is number of points in each layer, iz is resulting Newton matrix size c c c variables on upper wall c real *8 uu(-ixm:ixp,0:iy),duu(0:iyp) real *8 wu(-ixm:ixp,0:iy),dwu(0:iyp) real *8 psiu(-ixm:ixp,0:iy),dpsiu(0:iyp) real *8 vu(-ixm:ixp,0:iy) real *8 pu(-ixm:ixp),pxu(-ixm:ixp),dpxu c c variables on lower wall c 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 f(-ixm:ixp),g(-ixm:ixp) real *8 a(-ixm-2:ixp+2),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 mu0,gamma0,lambda0,zeta0,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=1. c set up initial guess tolerance=1.e-6 ixu=40 ixd=40 ipx=1 ipy=1 msize=ixm msizep=ixp iterations=50 c c set up mesh c hy0=.1 c hy0=0.025 hx0=1./80 c hy0=hx0 facy=1.1 c facy=1.05 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 facd=1 c facu=1.05 c facd=1.05 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, mu' 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 c for asymmetric pressure distribution, start a(x) as -(g+f)/2 f(i)=ys(i,0) g(i)=0. a(i)=-(g(i)+f(i))/2 c a(i)=0.5*a(i) p(i)=0 px(i)=0 pu(i)=0 pxu(i)=0 end do a(-ixu)=a(ixd) a(-ixu-1)=a(ixd-1) a(-ixu-2)=a(ixd-2) a(ixd+1)=a(-ixu+1) a(ixd+2)=a(-ixu+2) 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 uu(i,0)=0. psiu(i,0)=0. wu(i,0)=fluxcorrection/2 do j=1,iy z=y(j) psi(i,j)=fluxcorrection*(z*z/2+z*(a(i)+f(i)))/2 u(i,j)=fluxcorrection*(z+a(i)+f(i))/2 w(i,j)=fluxcorrection/2 psiu(i,j)=fluxcorrection*(z*z/2-z*(a(i)+g(i)))/2 uu(i,j)=fluxcorrection*(z-(a(i)+g(i)))/2 wu(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 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*fluxcorrection rhs(row)=u(i,iy)-fluxcorrection*(y(iy)+a(i)+f(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*fluxcorrection rhs(row)=uu(i,iy)-fluxcorrection*(y(iy)-a(i)-g(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+1)-3*a(i)+3*a(i-1)-a(i-2))/(hx(i)**3) row=row+1 cof(row,row-2)=1. cof(row,row-1)=-1. cof(row,row)= -3.*fluxcorrection**2/(120.*hx(i)**3) rhs(row)=pxu(i)-px(i)-fluxcorrection**2*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(ixd,iy) + +abs(du(j))/u(ixd,iy) + +abs(dw(j))/w(ixd,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(ixd,iy) + +abs(duu(j))/u(ixd,iy) + +abs(dwu(j))/w(ixd,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 do j=0,iy psi(-ixu,j)=psi(ixd,j) u(-ixu,j)=u(ixd,j) w(-ixu,j)=w(ixd,j) enddo do j=0,iy psiu(-ixu,j)=psiu(ixd,j) uu(-ixu,j)=uu(ixd,j) wu(-ixu,j)=wu(ixd,j) enddo a(-ixu)=a(ixd) a(-ixu-1)=a(ixd-1) a(-ixu-2)=a(ixd-2) a(ixd+1)=a(-ixu+1) a(ixd+2)=a(-ixu+2) c c -------------- estimate a'' at left boundary c i=-ixu 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 p(-ixu)=0. pu(-ixu)=p(-ixu)+fluxcorrection**2*add/120 do i=-ixu+1,ixd p(i)=p(i-1)+hx(i)*px(i) pu(i)=pu(i-1)+hx(i)*pxu(i) enddo c write(6,669)a(-ixu-1),a(-ixu),a(-ixu+1),p(-ixu),add,pu(-ixu) 669 format(6e14.4) write(6,10)k,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) write(50,67)(x(i),wu(i,0),i=-ixu,ixd) write(46,67)(x(i),pu(i),i=-ixu,ixd) write(42,21)((psiu(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(43,21)((uu(i,j),j=0,iy,ipy),i=-ixu,ixd,ipx) write(44,21)((wu(i,j),j=0,iy,ipy),i=-ixu,ixd,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=-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 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