搜索
programmainimplicitnone!----------------integer,parameter::nx=101,ny=51real*8::p=0.3,q=0.3real*8::dks=1.0,dat=1.0real*8,dimension(nx,ny)::x=0,y=0,xp=0,yp=0real*8,allocatable::xb0(:),yb0(:),xb1(:),yb1(:)!---------------integer::nb0,nb1integer::i,ii,j,times,di,dj!----------------real*8::errmax=0.5,errreal*8::xks,xat,yks,yat,xksat,yksat,arfa,beta,gama,jjreal*8::a,ap,b,bp,c,cp,d,dp,e,epreal*8::l,lp,m,mp,n,npreal*8::xe,xw,xn,xs,ye,yw,yn,ysreal*8::xx,yy!---------------real*8::thta,thta1,thta2,k1,k2,dis,dis1,disa,disb,discreal*8::resinteger::iffind!------------------------------------------------------------------------------------------------!------------------------------------------------------------------------------------------------!边界坐标open(1,file='boundary.txt')read(1,*)nb0!左岸边界allocate(xb0(nb0),yb0(nb0))read(1,*)((xb0(i),yb0(i)),i=1,nb0)read(1,*)nb1!右岸岸边界allocate(xb1(nb1),yb1(nb1))read(1,*)((xb1(i),yb1(i)),i=1,nb1)close(1)!initial----------------------------------------------------------------------------------------callcboundary(nx,nb0,xb0,yb0,x(:,1),y(:,1))!(1,1)-(nx,1)河床左岸边界坐标callcboundary(nx,nb1,xb1,yb1,x(:,ny),y(:,ny))!(1,ny)-(nx,ny)河床右岸边界坐标!------------------------------------doi=1,nx,nx-1!(1,1)-(1,ny),(nx,1)-(nx,ny)河床上下边界坐标doj=2,ny-1x(i,j)=x(i,1)+(x(i,ny)-x(i,1))*(j-1)/(ny-1)y(i,j)=y(i,1)+(y(i,ny)-y(i,1))*(j-1)/(ny-1)enddoenddodoj=2,ny-1!(x,y)河床内部节点边界doi=2,nx-1x(i,j)=x(1,j)+(x(nx,j)-x(1,j))*(i-1)/(nx-1)y(i,j)=y(1,j)+(y(nx,j)-y(1,j))*(i-1)/(nx-1)enddoenddodoi=2,nx-1doj=2,ny-1x(i,j)=x(i,1)+(x(i,ny)-x(i,1))*(j-1)/(ny-1)y(i,j)=y(i,1)+(y(i,ny)-y(i,1))*(j-1)/(ny-1)enddoenddo!--------open(1,file='initial.dat')write(1,*)'variable=x,y'write(1,*)'zonei=',ny,'j=',nxdoi=1,nxdoj=1,nywrite(1,*)x(i,j),y(i,j)enddoenddoclose(1)!-------------------------------------------------------------------------------------------------!-------------------------------------------------------------------------------------------------times=0100continueerr=0doi=2,nx-1doj=2,ny-1xks=(x(i+1,j)-x(i-1,j))/(2*dks)xat=(x(i,j+1)-x(i,j-1))/(2*dat)yks=(y(i+1,j)-y(i-1,j))/(2*dks)yat=(y(i,j+1)-y(i,j-1))/(2*dat)xksat=(x(i+1,j+1)+x(i-1,j-1)-x(i+1,j-1)-x(i-1,j+1))/(4*dks*dat)yksat=(y(i+1,j+1)+y(i-1,j-1)-y(i+1,j-1)-y(i-1,j+1))/(4*dks*dat)!----------arfa=xat*xat+yat*yatgama=xks*xks+yks*yksbeta=xksat+yksatjj=xks*yat-xat*yks!考虑P、Q-----------------------------------------------------------------------!a=arfa*yks*yks/(xks*xks+yks*yks)!b=gama*yat*yat/(xat*xat+yat*yat)!c=-arfa*xks*yks/(xks*xks+yks*yks)!d=-gama*xat*yat/(xat*xat+yat*yat)!e=-2*beta*xksat!ap=c!bp=d!cp=arfa*xks*xks/(xks*xks+yks*yks)!dp=gama*xat*xat/(xat*xat+yat*yat)!ep=-2*beta*yksat!!--------------------------------!xe=x(i+1,j)!xw=x(i-1,j)!xn=x(i,j+1)!xs=x(i,j-1)!ye=y(i+1,j)!yw=y(i-1,j)!yn=y(i,j+1)!ys=y(i,j-1)!!--------------------------------!l=a*(xe+xw)/dks/dks+b*(xn+xs)/dat/dat+c*(ye+yw)/dks/dks+d*(yn+ys)/dat/dat+e!m=2*(a+b)!n=2*(c+d)!lp=ap*(xe+xw)/dks/dks+bp*(xn+xs)/dat/dat+cp*(ye+yw)/dks/dks+dp*(yn+ys)/dat/dat+ep!mp=2*(ap+bp)!np=2*(cp+dp)!!----------------------------!xp(i,j)=(l*np-lp*n)/(m*np-mp*n)!yp(i,j)=(l*mp-lp*m)/(mp*n-m*np)!if(abs(xp(i,j)-x(i,j))err)err=abs(xp(i,j)-x(i,j))!if(abs(yp(i,j)-y(i,j))err)err=abs(yp(i,j)-y(i,j))!x(i,j)=xp(i,j)!y(i,j)=yp(i,j)!不考虑P，Q----------------------------------------------------------------------p=0.q=0.!---------xp(i,j)=arfa*(x(i-1,j)+x(i+1,j))/dks/dks+gama*(x(i,j+1)+x(i,j-1))/dat/dat-2*beta*xksat+jj*(p*xks+q*xat)xp(i,j)=xp(i,j)/(2.0*(arfa/dks/dks+gama/dat/dat))if(abs(xp(i,j)-x(i,j))err)err=abs(xp(i,j)-x(i,j))x(i,j)=xp(i,j)!---------yp(i,j)=arfa*(y(i-1,j)+y(i+1,j))/dks/dks+gama*(y(i,j+1)+y(i,j-1))/dat/dat-2*beta*yksat+jj*(p*yks+q*yat)yp(i,j)=yp(i,j)/(2.0*(arfa/dks/dks+gama/dat/dat))if(abs(yp(i,j)-y(i,j))err)err=abs(yp(i,j)-y(i,j))y(i,j)=yp(i,j)enddoenddo!--------------------------------------------------------------------------------goto99!边界滑移------------------------------------------------------------------------!xboundary----!doj=1,ny,ny-1!if(j==1)dj=1!bottom!if(j==ny)dj=-1!top!doi=2,nx-1!!print*,i,j!!寻找边界控制点坐标!xx=0!yy=0!iffind=0!if(j==1)then!doii=1,nb0-1!if(x(i,j)=min(xb0(ii),xb0(ii+1)).and.x(i,j)max(xb0(ii+1),xb0(ii)))then!xx=xb0(ii)!yy=yb0(ii)!iffind=1!cycle!endif!enddo!else!doii=1,nb1-1!if(x(i,j)=xb1(ii).and.x(i,j)xb1(ii+1))then!xx=xb1(ii)!yy=yb1(ii)!iffind=1!cycle!endif!enddo!endif!if(iffind/=1)then!print*,'找不到边界控制点坐标'!stop!endif!!-----------------------------------------------------------------!!x(i,j+dj)!!/*!!//|!!c//|!!/a/|!!//|!!/b/thta|!!---------*---------------*--------------*----------!!(xx,yy)x(i,j)(xx1,yy2)!!------------------------------------------------------------------!disa=sqrt((x(i,j+dj)-x(i,j))**2.0+(y(i,j+dj)-y(i,j))**2.0)!disb=sqrt((xx-x(i,j))**2.0+(yy-y(i,j))**2.0)!disc=sqrt((x(i,j+dj)-xx)**2.0+(y(i,j+dj)-yy)**2.0)!if(disb/=0.and.disa/=0)then!thta=(disa*disa+disb*disb-disc*disc)/(2*disa*disb)!if(abs(thta)1)thta=thta/abs(thta)!thta=acos(thta)!thta=3.1415926-thta!x(i,j)=x(i,j)+(x(i,j)-xx)*disa*cos(thta)/disb!y(i,