program tdeckn c c Triple deck solver using Newton iteration c implicit real*8 (a-h,o-z) implicit integer (i-n) parameter (iy=64,iyp=iy+1,iz=3*(iy+1)+1,ixm=200,ixp=200, + ixm2=2*ixm,ixp2=2*ixp) common/wakedata/h0(0:2000),dh0(0:2000),ddh0(0:2000),hwake,nwake 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(-ixm2:ixp2),px(-ixm2:ixp2),dpx real *8 a(-ixm2:ixp2),ad(-ixm2:ixp2) real *8 a0(-ixm:ixp) real *8 cof(iz,iz),rhs(iz) real *8 x(-ixm2:ixp2),hx(-ixm2:ixp2),hy(0:iy),y(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 integer ixu,ixd,ipx,ipy,iy,i,j,k,msize,msizep,row,nits c set up initial guess real *8 gamma mu0 = 1.132138751 gamma0 = 0.244547191 lambda0= 0.87890145 zeta0=mu0/lambda0 third=1./3. nits=5 a1=.3321 xscale=a1**(1.25) pscale=a1**(-.5) ascale=a1**.75 tolerance=1.e-6 ixu=200 ixd=200 ipx=5 ipy=1 msize=ixm2 msizep=ixp2 iterations=30 c c set up y mesh c hy0=.05 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) c c set up x mesh c hx0=0.02 facu=1.05 facd=1.05 x(0)=0. x(1)=hx0 hx(1)=hx0 hx(0)=hx0 x(-1)=-hx0 do i=2,msizep if(i.le.40)then hx(i)=facd*hx(i-1) else hx(i)=1.01*hx(i-1) endif x(i)=x(i-1)+hx(i) enddo do i=2,msize if(i.le.40)then hx(-i+1)=facu*hx(-i+2) else hx(-i+1)=1.01*hx(-i+2) endif x(-i)=x(-i+1)-hx(-i+1) enddo xm=x(ixd) write(7,7)(x(i),hx(i),i=-ixu,ixd,ipx) write(8,7)(y(j),hy(j),j=0,iy,ipy) 7 format(2e15.5) write(6,5)ixu,ixd,iy,zm,xm,hx0,hy0,msize,msizep 5 format('Interactive bl 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) gamma=zeta0/(3.*3.**.5) alpha=2*lambda0*lambda0*lambda0*gamma0 c c set up initial guess for a(x),p(x) c write(6,1) 1 format('Initialising a(x), p(x)') do i=-msize,-1 ad(i)=0. a(i)=0. xt=(-x(i))**third p(i)=-2*gamma/(xt*xt) xc=-(x(i)-0.5*hx(i)) xt=xc**third px(i)=-2*2*gamma/(3*xc*xt*xt) end do ad(0)=0 a(0)=0 p(0)=0 xc=0.5*hx(0) xt=xc**third px(0)=-2*2*gamma/(3*xc*xt*xt) do i=1,msizep xt=(x(i))**third ad(i)=zeta0/(3*xt*xt) a(i)=zeta0*xt p(i)= gamma/(xt*xt) xc=x(i)-0.5*hx(i) xt=xc**third px(i)=2*gamma/(3*xc*xt*xt) end do c c now have asymptotic form for a,p c write(6,2) 2 format('Calling Hilbert transform') call hilb(ad,msize,msizep,p,ixu,ixd,hx,x) c c calculate how numerical hilbert transform sees a(x) from p(x) c a(-ixu)=0 do i=-ixu+1,ixd a(i)=a(i-1)+hx(i)*(ad(i-1)+ad(i))/2 enddo write(6,3) 3 format('Setting up plate start solution') do i=-ixu,0 u(i,0)=0. psi(i,0)=0. w(i,0)=alpha do j=1,iy z=y(j) psi(i,j)=alpha*(z*z/2+z*a(i)) u(i,j)=alpha*(z+a(i)) w(i,j)=alpha end do end do sm=zm/(xm**third) write(6,4) 4 format('Setting up wake start solution') c c intialise wake solution c call wake0 do i=1,ixd xt=x(i)**third do j=0,iy z=y(j) s=z/xt if(s.le.sm)then call wake(s,ho,dho,ddho) psi(i,j)=xt*xt*ho u (i,j)=xt*dho w(i,j)=ddho else u(i,j)=alpha*(z+a(i)) psi(i,j)=alpha*(z*z/2+z*a(i)) w(i,j)=alpha endif end do end do c c now have initial guess for all variables, c output intial values c 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)(xscale*x(i),ascale*a(i), i=-msize,msizep) write(16,65)(xscale*x(i),pscale*p(i), i=-msize,msizep) write(17,65)(x(i),ad(i), i=-msize,msizep) write(18,65)(x(i)-0.5*hx(i),px(i),i=-msize+1,msizep) 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,0),(x(i),u(i,0),i=1,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,ixd u(i,iy)=alpha*(zm+a(i)) enddo do i=-ixu+1,0 c-------------------------------------------------------------------- c iterate along plate rows - u(0)=0 | c-------------------------------------------------------------------- do newton=1,nits 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 c ---------------------------------------------------------------- c iterate along wake | c------------------------------------------------------------------ do i=1,ixd do newton=1,10 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)=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 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,71)error/(3*iy) 71 format('wake ',e12.5) dpx=rhs(iz) px(i)=px(i)+dpx enddo enddo c c get pressures from gradients, gradients at midpoints c do i=-ixu+1,ixd p(i)=p(i-1)+hx(i)*px(i) enddo c c ---------------- use Hilbert transform to get a' from p(x) c call hilb(ad,msize,msizep,p,ixu,ixd,hx,x) c c calculate a(x) from a'(x) c a(-ixu)=0. omega=.05 a0(-ixu)=0. do i=-ixu+1,ixd a0(i)=a0(i-1)+hx(i)*(ad(i-1)+ad(i))/2 a(i)=(1-omega)*a(i)+omega*a0(i) enddo write(6,10)k,a(-ixu/2),a(0),a(ixd/4), + a(ixd/2),a(3*ixd/4),a(ixd),p(0),px(1) 10 format(i6,(10e12.3)) 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)(xscale*x(i), ascale*a(i), i=-ixu,ixd) write(26,65)(xscale*x(i),pscale*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,0),(x(i),u(i,0),i=1,ixd) 21 format(8e15.4) c c calculate drag c drag=0. do i=-ixu,-1 drag=drag+(w(i+1,0)+w(i,0)-2*alpha)*(x(i+1)-x(i))/2 enddo write(6,25)drag 25 format('Calculated drag coefficient', f15.6) stop end c c-------------------------------------------------------------------------------hilb------------- c subroutine hilb(ad,msize,msizep,p,ixu,ixd,hx,x) real *8 p(-msize:msizep),ad(-msize:msizep) real *8 x(-msize:msizep),hx(-msize:msizep) integer msize,msizep,ixu,ixd real *8 mu0,gamma0,lambda0,zeta0 real *8 phx,k1,sum,b,t1,t2,dx mu0 = 1.132138751 gamma0 = 0.244547191 lambda0= 0.87890145 zeta0=mu0/lambda0 phx=1./3.14159 do i=-ixu,ixd c pressure at grid points sum1=0. sum=hx(i+1)*p(i+1)/(x(i)-x(i+1)) do j=i+1,2*ixd-1 sum1=sum1+0.5*(hx(j)+hx(j+1))*p(j)/(x(i)-x(j)) enddo sum1=sum1+0.5*p(2*ixd)*hx(2*ixd)/(x(i)-x(2*ixd)) sum1=sum1*phx b=x(2*ixd) t2=tail(x(i),b,100)*zeta0/(3**1.5) j=-2*ixu sum2=0.5*hx(j+1)*p(j)/(x(i)-x(j)) do j=-2*ixu+1,i-1 sum2=sum2+0.5*(hx(j)+hx(j+1))*p(j)/(x(i)-x(j)) enddo sum2=sum2*phx b=-x(-2*ixu) t1=tail(-x(i),b,100)*2*zeta0/(3**1.5) ad(i)=-(sum1+sum2+t1+t2) enddo return end real*8 function tail(x,b,m) integer i,k,m real *8 b1,bx,x,b b1=3.*b**(-2./3.)/3.14159 k=3*m+2 bx=x/b tail=1./k do i=1,m k=k-3 tail=1./k+bx*tail enddo tail=-b1*tail return end subroutine wake(s,ho,dho,ddho) implicit real*8(a-h,o-z) implicit integer (i-n) common/wakedata/h0(0:2000),dh0(0:2000),ddh0(0:2000),h,n real *8 a0,b0,df,e real *8 fa,fb,fc,fd,dfa,dfb,dfc,dfd,ddfa,ddfb,ddfc,ddfd real *8 sp,sm,mu0,mu1,gamma0,gamma1,lambda0 integer i,j,k,l,m,n lambda0=0.87890145 z=lambda0*s h=0.01 n=1000 if(z.gt.n*h)write(6,1)s,z 1 format('panic in wake', 2f15.6) fc=n*h gamma0=ddh0(n)/2 fa=h0(n)/gamma0 fb=fa**0.5 mu0=fb-fc np=z/h dt=z-np*h dt1=1-dt ho=h0(np)*dt1+dt*h0(np+1) dho=dh0(np)*dt1+dt*dh0(np+1) ddho=ddh0(np)*dt1+dt*ddh0(np+1) ho=lambda0*ho dho=lambda0*lambda0*dho ddho=lambda0*lambda0*lambda0*ddho return end subroutine wake0 implicit real*8(a-h,o-z) implicit integer (i-n) common/wakedata/h0(0:2000),dh0(0:2000),ddh0(0:2000),h,n real *8 a0,b0,df,e real *8 fa,fb,fc,fd,dfa,dfb,dfc,dfd,ddfa,ddfb,ddfc,ddfd real *8 sp,sm,mu0,mu1,gamma0,gamma1,lambda0 integer i,j,k,l,m,n h=0.01 n=1000 h0(0)=0. dh0(0)=1. ddh0(0)=0. do i=1,n fa= h*dh0(i-1) dfa= h*ddh0(i-1) ddfa= h*(dh0(i-1)*dh0(i-1)-2*h0(i-1)*ddh0(i-1))/3 fb= h*(dh0(i-1)+dfa/2) dfb= h*(ddh0(i-1)+ddfa/2) ddfb= h*((dh0(i-1)+dfa/2)*(dh0(i-1)+dfa/2)- + 2*(h0(i-1)+fa/2)*(ddh0(i-1)+ddfa/2))/3 fc= h*(dh0(i-1)+dfb/2) dfc= h*(ddh0(i-1)+ddfb/2) ddfc= h*((dh0(i-1)+dfb/2)*(dh0(i-1)+dfb/2)- + 2*(h0(i-1)+fb/2)*(ddh0(i-1)+ddfb/2))/3 fd= h*(dh0(i-1)+dfc) dfd= h*(ddh0(i-1)+ddfc) ddfd= h*((dh0(i-1)+dfc)*(dh0(i-1)+dfc)- + 2*(h0(i-1)+fc)*(ddh0(i-1)+ddfc))/3 h0(i)=h0(i-1)+(fa+2*fb+2*fc+fd)/6 dh0(i)=dh0(i-1)+(dfa+2*dfb+2*dfc+dfd)/6 ddh0(i)=ddh0(i-1)+(ddfa+2*ddfb+2*ddfc+ddfd)/6 enddo return end subroutine invert(a,n,np,b) implicit real *8(a-h,o-z) real*8 a(np,np),b(np) integer indx(500) c c---------------------------------------------------------------- c c You must provide matrix inversion here. For instance c using routine from Numerical Recipes c c---------------------------------------------------------------- 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 return end