program flow * * driver for kirchhoff program kirch. * computes free-streamline flow over a polygonal obstacle. * l. n. trefethen, summer 1984 * implicit double precision (a-b,d-h,o-w,y) implicit complex*16 (c,x,z) common /dual/ sgx dimension z(20),x(20),betam(20),qwork(400) dimension x2(20),z2(20),betam2(20),qwork2(400) data zi,zoffst /(0.d0,1.d0),(2.d0,0.d0)/ sgx = 1. * * input control parameters: 1 print '(''n, nq?'')' read *,n,nq do 2 k = 1,n print '(''z('',i1,'') ?'')', k read *, aa,bb 2 z(k) = dcmplx(aa,bb) tol = 10.**(-min(12,nq+2)) * * solve mapping problem: call count0 call qinitx(n,z,betam,nq,qwork) call ksolv(0,tol,n,c,x,z,betam,nq,qwork) aw = abs(c) call count * * print results: write (10,'(''w'',4f10.5)') -3.,-3.,3.,3. print '(''stagnation point:'',2f10.7)', z(n+1) write (10,'(''m -2.5 -2.3''/''t 17 stagnation point:'')') write (10,'(''m -1. -2.3''/''t 20 '',2f10.6)') z(n+1) print '(''w:'',f10.7)', aw write (10,'(''m -2.5 -2.6''/''t 17 w:'')') write (10,'(''m -1. -2.6''/''t 10 '',f10.6)') aw * * force computation 1 -- along obstacle: sgx = -1. call qinitx(n,z,betam2,nq,qwork2) x2(1) = x(1) z2(1) = z(1) write (10,101) z2(1) + zoffst kmq = 1 do 12 k = 2,n km = k-1 if (betam2(k).ge.-.9999) then * --- non-sharp edge --- x2(k) = x(k) z2(k) = zx(x2(k),k,x2(km),z2(km),kmq,n,c,x,betam2,nq,qwork2) write (10,102) z2(k) + zoffst kmq = k else * --- sharp edge --- scale = min( abs(x(k)-x(km)) , abs(x(k+1)-x(k)) ) xx = x(k) - 1.d-2*scale zz = zx(xx,0,x2(km),z2(km),kmq,n,c,x,betam2,nq,qwork2) write (10,102) zz + zoffst xx2 = x(k) - .1*zi*scale zz = zx(xx2,0,xx,zz,0,n,c,x,betam2,nq,qwork2) x2(k) = x(k) + 1.d-2*scale z2(k) = zx(x2(k),0,xx2,zz,0,n,c,x,betam2,nq,qwork2) write (10,101) z2(k) + zoffst kmq = 0 endif 12 continue zf1 = zi * (z2(n)-z(n)) * * force computation 2 -- via point x = -i : ztmp = zx(-zi,0,x(1),z(1),1,n,c,x,betam2,nq,qwork2) ztmp = zx(x(n),n,-zi,ztmp,0,n,c,x,betam2,nq,qwork2) zf2 = zi * (ztmp-z(n)) * * force computation 3 -- via point x = 1.1-i : zz = (1.1d0,-1.d0) ztmp = zx(zz,0,x(1),z(1),1,n,c,x,betam2,nq,qwork2) ztmp = zx(x(n),n,zz,ztmp,0,n,c,x,betam2,nq,qwork2) zf3 = zi * (ztmp-z(n)) * * print forces; all three sets of numbers should agree: print '('' drag ='',3f12.6)', dreal(zf1),dreal(zf2),dreal(zf3) print '('' lift ='',3f12.6)', dimag(zf1),dimag(zf2),dimag(zf3) print '('' total ='',3f12.6)', abs(zf1),abs(zf2),abs(zf3) write (10,'(''m -2.5 -2.''/''t 24 drag, lift, total force:'')') write (10,'(''m 0. -2.''/''t 30 '',3f10.6)') dreal(zf2), & dimag(zf2), abs(zf2) sgx = 1. * * plot free streamlines: 23 write (10,101) z(1) write (10,102) (z(k),k=2,n) print '('' no. of streamlines ?'')' read *, nstr if (nstr.eq.0) stop space = .1d0*(2.d0/aw) npts = 30 do 4 nl = 0,1 write (10,101) z(n-nl*(n-1)) do 3 k = 0,npts xx = (-1.)**nl * (1. + k**3*sqrt(32.*space)/npts**3) call nearx(xx,x0,z0,k0,n,x,z,betam) zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork) 3 write (10,102) zz 4 continue print '(''finished free streamlines'')' * * plot internal streamlines: do 6 nl = -nstr,nstr do 5 k = -npts,npts if (nl.eq.0.and.k.gt.0) goto 6 xx = sqrt(dcmplx(space*k**3*54./npts**3,nl*space)) if (dimag(xx).lt.0.d0) xx = -xx xx = xx + x(n+1) call nearx(xx,x0,z0,k0,n,x,z,betam) zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork) if (k.eq.-npts) write (10,101) zz write (10,102) zz 5 continue 6 continue print '(''finished internal streamlines'')' * * plot upstream equipotential lines: do 8 nl = -16,-1 klim = 2*nstr if (nl.lt.-7) klim = nstr do 7 k = -klim,klim xx = sqrt(dcmplx(2.*nl*space,k*nstr*space/klim)) if (dimag(xx).lt.0.d0) xx = -xx xx = xx + x(n+1) call nearx(xx,x0,z0,k0,n,x,z,betam) zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork) if (k.eq.-klim) write (10,101) zz write (10,102) zz 7 continue 8 continue print '(''finished upstream equipotential lines'')' * * plot downstream equipotential lines: do 11 nn = 1,2 do 10 nl = 0,27 klim = 2*nstr if (nl.eq.0) klim = 10*klim if (nl.gt.5) klim = nstr do 9 k = 0,klim xx = sqrt(dcmplx(2.*nl*space,k*(-1)**nn*nstr*space/klim)) if (dimag(xx).lt.0.d0) xx = -xx if (nn.eq.1.and.k.eq.0) xx = -xx xx = xx + x(n+1) call nearx(xx,x0,z0,k0,n,x,z,betam) zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork) if (k.eq.0) write (10,101) zz write (10,102) zz 9 continue 10 continue 11 continue print '(''finished downstream equipotential lines'')' call count 101 format ('m ',2f9.4) 102 format ('c ',2f9.4) end * kirch1 * * computes the conformal map of the upper half plane (x) onto an * asymmetrical free-streamline flow domain (z) of kirchhoff type. * x is the square root of the velocity potential. * the obstacle is composed of n-1 straight line segments with * vertices z(k) indexed from 1 to n. z(n+1) is the stagnation point. * see monakhov, "boundary-value problems with free boundaries...", iv.1. * * adapted from scpack. l. n. trefethen, june 1984. subroutine qinitx(n,z,betam,nq,qwork) * * computes nodes and weights for gauss-jacobi quadrature * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension qwork(1),z(1),betam(1) call anglex(n,z,betam) iscr = nq*(2*n+2) + 1 do 1 k = 1,n+1 inodes = nq*(k-1) + 1 iwts = nq*(n+k) + 1 1 if (betam(k).ge.-.9999) call gaussj(nq,0.d0, & betam(k),qwork(iscr),qwork(inodes),qwork(iwts)) return end subroutine ksolv(iprint,tol,n,c,x,z,betam,nq,qwork) * * directs the solution of the modified schwarz-christoffel parameter * problem for kirchhoff flow over a polygonal obstacle * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) common /param1/ nq2,c2,x2(20),z2(20),qwork2(460),betam2(20) dimension x(1),z(1),betam(1),qwork(1) dimension ajinv(20,20),scr(900),fval(19),y(19) external kfun np = n+1 nm = n-1 nm2 = n-2 x(1) = -1. x(n) = 1. * * initial guess and ns01a control parameters: do 1 k = 2,nm 1 y(k-1) = 0. dstep = 1.d-6 dmax = 300. maxfun = 30*nm * * copy input data to /param1/ and solve nonlinear system: nq2 = nq do 2 k = 1,np x2(k) = x(k) z2(k) = z(k) 2 betam2(k) = betam(k) nwdim = nq*(2*n+3) do 3 i = 1,nwdim 3 qwork2(i) = qwork(i) call ns01a(nm2,y,fval,ajinv,dstep,dmax,tol,maxfun, & iprint,scr,kfun) * * copy output data from /param1/ and print results: c = c2 do 4 k = 2,np 4 x(k) = x2(k) call nearx(x(np),x0,z0,k0,n,x,z,betam) z(np) = zx(x(np),0,x0,z0,k0,n,c,x,betam,nq,qwork) if (iprint.ge.0) call koutp(n,c,x,z,betam,nq) return end function zx(xx,kxx,x0,z0,k0,n,c,x,betam,nq,qwork) * * computes map z(xx) * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension x(1),betam(1),qwork(1) zx = z0 + c * xquad(x0,k0,xx,kxx,n,x,betam,nq,qwork) return end subroutine yxtran(n,y,x,z,betam) * * transforms y(k) to x(k) * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension y(1),x(1),z(1),betam(1) nm = n-1 dx = 1. sumx = dx do 1 k = 2,nm dx = dx / exp(y(k-1)) 1 sumx = sumx + dx dx = 2./sumx x(2) = -1. + dx do 2 k = 3,nm dx = dx / exp(y(k-2)) 2 x(k) = x(k-1) + dx asum = 0. do 3 k = 2,nm 3 asum = asum + betam(k)*acos(dreal(-x(k))) x(n+1) = -cos( dimag(log(z(n)-z(nm))) + asum ) return end subroutine kfun(nm2,y,fval) * * this is the function whose zero must be found in ksolv * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension fval(nm2),y(nm2) common /param1/ nq,c,x(20),z(20),qwork(460),betam(20) nm = nm2+1 n = nm2+2 call yxtran(n,y,x,z,betam) zdenom = xquad(x(1),1,x(2),2,n,x,betam,nq,qwork) c = (z(2)-z(1))/zdenom do 1 k = 2,nm kp = k+1 xint = xquad(x(k),k,x(kp),kp,n,x,betam,nq,qwork) fval(k-1) = abs(z(kp)-z(k)) - abs(c*xint) 1 continue return end subroutine koutp(n,c,x,z,betam,nq) * * prints a table of k, z(k), betam(k), and x(k) * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension x(1),z(1),betam(1) print 101, n, nq do 1 k = 1,n 1 print 102, k,z(k),betam(k),x(k) print 103, z(n+1), x(n+1) print '(/'' c = ('',e15.8,'','',e15.8,'')''/)', c return 101 format (/' parameters defining map:',15x,'(n =', & i3,')',6x,'(nq =',i3,')'// & ' k', 9x,'z(k)',11x,'betam(k)',12x,'x(k)'/ & ' ---',8x,'----',11x,'--------',12x,'----'/) 102 format (i3,' (',f6.3,',',f6.3,')',f12.5, & 3x,'(',f11.8,',',f11.8,')') 103 format ('stagn. (',f6.3,',',f6.3,')',15x,'(',f11.8,',',f11.8,')') end function xquad(xa,ka,xb,kb,n,x,betam,nq,qwork) * * computes the complex line integral of xprod from xa to xb along a * straight line segment. function xquad1 is called twice. * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension x(1),betam(1),qwork(1) xmid = (xa + xb) / 2. xquad = xquad1(xa,xmid,ka,n,x,betam,nq,qwork) & -xquad1(xb,xmid,kb,n,x,betam,nq,qwork) return end function xquad1(xa,xb,ka,n,x,betam,nq,qwork) * * computes the complex line integral of xprod from xa to xb along a * straight line segment by compound one-sided gauss-jacobi quadrature. * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension x(1),betam(1),qwork(1) data resprm /1.d0/ * * step 1: one-sided gauss-jacobi quadrature for left endpoint: if (abs(xa-xb).gt.0.d0) goto 1 xquad1 = (0.,0.) return 1 r = min(1.d0,distx(xa,ka,n,x)*resprm/abs(xa-xb)) xaa = xa + r*(xb-xa) xquad1 = xqsum(xa,xaa,ka,n,x,betam,nq,qwork) * * step 2: adjoin intervals of pure gaussian quadrature if necessary: 2 if (r.eq.1.d0) return r = min(1.d0,distx(xaa,0,n,x)*resprm/abs(xaa-xb)) xbb = xaa + r*(xb-xaa) xquad1 = xquad1 + xqsum(xaa,xbb,0,n,x,betam,nq,qwork) xaa = xbb goto 2 end function distx(xx,ks,n,x) * * determines the distance from xx to the nearest singularity x(k) * other than x(ks) * implicit double precision (a-b,d-h,o-w,y) implicit complex*16(c,x,z) dimension x(1) distx = 9999. do 1 k = 1,n if (k.eq.ks) goto 1 distx = min(distx,abs(xx-x(k))) 1 continue return end function xqsum(xa,xb,ka,n,x,betam,nq,qwork) * * computes the integral of xprod from xa to xb by applying a * one-sided gauss-jacobi formula with possible singularity at xa. * ka=0: gauss-legendre * ka=1 or n: gauss-legendre after change of variables (separation pt.) * 1