求2维贴体正交网格生成算例

ak47long

对解泊松方程构造贴体正交网格很感兴趣，希望各位大牛提供这方面的算例，最好有步骤有程序，或者给小弟推荐几本介绍网格生成技术的书，小弟的论文要自己编程，所以软件运用方面的书就算了，先谢谢各位

PS：THOMPSON的几篇经典文献小弟已经读过，可惜资质驽钝，难以应用，希望能够和大牛交流下

weixing1531

这是我自己编的程序

高手不要嘲笑哈


program mygrid
  implicit none

  integer::imax,jmax,imax1,jmax1,n,i,j
  real(8),parameter::pi=3.14159265358979324d0,eps0=1.0d-7 ! eps0精度控制值
  real(8)::xet,xxi,yet,yxi,xxiet,yxiet,af,bt,gm,fx,fy,xyj,xad,yad,eps,w,dx,dy
!   real(8)::an,as,aw,ae,ac,bx,by ! 第二种方法时用到的变量
  real(8),allocatable::q( : ),p( : )
  real(8),allocatable::x( : , : ),y( : , : ),xz( : , : ),xzz( : , : ),xe( : , : ),xee( : , : ),yz( : , : ),&
  yzz( : , : ),ye( : , : ),yee( : , : ),xze( : , : ),yze( : , : ),arfa( : , : ),beta( : , : ),gama( : , : ),binj( : , : ),&
  errx( : , : ),erry( : , : )
  character(10)::name
  logical::torf

  write(*,*)'请输入文件名(不包含后缀名，10个字符之内):'
  read(*,*)name
  open(9,file=trim(adjustl(name))//'.dat',form='unformatted') ! 该文件已经由grid.exe生成
  read(9)imax,jmax  ! 网格的长与宽
  imax1=imax-1 ! 网格的长划分份数
  jmax1=jmax-1 ! 网格的宽划分份数
  allocate(x(imax,jmax),y(imax,jmax),xz(imax,jmax),xzz(imax,jmax),xe(imax,jmax),&
  xee(imax,jmax),yz(imax,jmax),yzz(imax,jmax),ye(imax,jmax),yee(imax,jmax),&
  xze(imax,jmax),yze(imax,jmax),arfa(imax,jmax),beta(imax,jmax),gama(imax,jmax),&
  binj(imax,jmax),errx(2:imax1,2:jmax1),erry(2:imax1,2:jmax1))
  allocate(q(jmax),p(imax))

  write(*,*)'是否自己设置松弛因子的值?'
  read(*,*)torf   ! T代表设置,F代表不设置
  w=4.0d0/(2.0d0+sqrt(4.0d0-(cos(pi/imax1)+cos(pi/jmax1)**2.0d0))) ! 默认使用公式
  if(torf.eqv..true.)then ! 当为T时
    write(*,*)'请输入自己设置松弛因子的值(0至2之间):'
    read(*,*)w
  end if

  do j=1,jmax
    do i=1,imax
      read(9)x(i,j),y(i,j) ! 读取网格的初始坐标
    end do
  end do

  open(11,file=trim(adjustl(name))//'.out') ! 一些有用的结果文件
  open(12,file=trim(adjustl(name))//'_tecplot.dat') ! tecplot格式文件

! 坐标变换的计算
  do j=2,jmax1
    q(j)=0.0 ! 调节因子Q,只迭加一项(从1到1)
  write(11,200)q(j)
  end do
  do i=2,imax1
    p(i)=0.0  ! 调节因子P
    write(11,200)p(i)
  end do

  n=0 ! 迭代次数计数器
  eps=0.0d0
  do ! 与原文不同,这是我改写的标准FORTRAN,下同
    do j=2,jmax1
      do i=2,imax1
        xet=x(i,j+1)-x(i,j-1) ! 详见原书第158页说明
        xxi=x(i+1,j)-x(i-1,j)
                yet=y(i,j+1)-y(i,j-1)
                yxi=y(i+1,j)-y(i-1,j)
                xxiet=x(i+1,j+1)-x(i+1,j-1)+x(i-1,j-1)-x(i-1,j+1)
                yxiet=y(i+1,j+1)-y(i+1,j-1)+y(i-1,j-1)-y(i-1,j+1)
                af=xet*xet+yet*yet ! 4α
                bt=xxi*xet+yxi*yet ! 4β
                gm=xxi*xxi+yxi*yxi ! 4γ

! 这是第二种计算方法,详见《计算流体力学基础》
!         xyj=.25d0*(xxi*yet-xet*yxi)
!         ae=af/4.0d0+xyj*xyj*p(i)/2.0d0
!                 aw=af/4.0d0-xyj*xyj*p(i)/2.0d0
!                 an=gm/4.0d0+xyj*xyj*q(j)/2.0d0
!                 as=gm/4.0d0-xyj*xyj*q(j)/2.0d0
!                 ac=(af+gm)/2.0d0
!                 bx=-bt/8.0d0
!                 by=-bx*yxiet
!                 bx=-bx*xxiet
!                 fx=(ae*x(i+1,j)+aw*x(i-1,j)+an*x(i,j+1)+as*x(i,j-1)+bx)/ac
!                 fy=(ae*y(i+1,j)+aw*y(i-1,j)+an*y(i,j+1)+as*y(i,j-1)+by)/ac
!                 dx=fx-x(i,j)
!                 dy=fy-y(i,j)
!                 errx(i,j)=abs(dx)
!         erry(i,j)=abs(dy)
!         x(i,j)=x(i,j)+dx*w  ! SOR迭代
!                  y(i,j)=y(i,j)+dy*w

                fx=(af*(x(i-1,j)+x(i+1,j))-(bt*xxiet/2.d0)+gm*(x(i,j+1)+x(i,j-1)))/(2.d0*(af+gm))
                fy=(af*(y(i-1,j)+y(i+1,j))-(bt*yxiet/2.d0)+gm*(y(i,j+1)+y(i,j-1)))/(2.d0*(af+gm))
                xyj=.25d0*(xxi*yet-xet*yxi)
                xad=((xyj)**2*(xet*q(j)+xxi*p(i)))/(af+gm)
                yad=((xyj)**2*(yet*q(j)+yxi*p(i)))/(af+gm)
                fx=fx+xad ! 公式(7-91)
                fy=fy+yad
                dx=fx-x(i,j)
                dy=fy-y(i,j)
                errx(i,j)=abs(dx) ! 迭代中变化值
        erry(i,j)=abs(dy)
                x(i,j)=x(i,j)+dx*w  ! SOR迭代
                y(i,j)=y(i,j)+dy*w

      end do
    end do
    eps=max(maxval(errx),maxval(erry)) ! 找出迭代中变化最大的值
    n=n+1  
    write(*,"('迭代次数及精确度分别为:',i4,f13.7)")n,eps  ! 迭代次数及精确度
    if(eps<eps0)then
      write(*,*)'迭代收敛!'
          exit   ! 精度满足,跳出do循环
        end if
        if(n>1000.and.eps>0.01d0)then ! 最大循环数为1000
          write(*,*)'迭代不收敛!'  
          stop  ! 避免死循环
        end if
  end do

  write(11,500)n
  do i=1,imax
    write(11,100)i
    write(11,200)(x(i,j),j=1,jmax)
    write(11,200)(y(i,j),j=1,jmax)
  end do

  write(12,'(A)')'Variables="x","y"'
  write(12,"('Zone i=',I3,' j=',I3,' f=point')")imax,jmax
  do j=1,jmax
    do i=1,imax
          write(12,'(g15.5,2x,g15.5)')x(i,j),y(i,j)
        end do
  end do

100 format('I=',i2,50x,'**X-COORDINATE**')
200 format(8f10.5)
400 format(9f5.2)
500 format(50x,'***** N=',i4,'*****')
! calculation of transformed coefficients
! 转化系数计算
  do i=2,imax1
    do j=2,jmax1
          xz(i,j)=0.5d0*(x(i+1,j)-x(i-1,j))
          xzz(i,j)=x(i+1,j)-2.0d0*x(i,j)+x(i-1,j)
          xe(i,j)=0.5d0*(x(i,j+1)-x(i,j-1))
          xee(i,j)=x(i,j+1)-2.0d0*x(i,j)+x(i,j-1)
          yz(i,j)=0.5d0*(y(i+1,j)-y(i-1,j))
          yzz(i,j)=y(i+1,j)-2.0d0*y(i,j)+y(i-1,j)
          ye(i,j)=0.5d0*(y(i,j+1)-y(i,j-1))
          yee(i,j)=y(i,j+1)-2.0d0*y(i,j)+y(i,j-1)
          xze(i,j)=0.25d0*(x(i+1,j+1)-x(i+1,j-1)+x(i-1,j-1)-x(i-1,j+1))
          yze(i,j)=0.25d0*(y(i+1,j+1)-y(i+1,j-1)+y(i-1,j-1)-y(i-1,j+1))
          arfa(i,j)=xe(i,j)*xe(i,j)+ye(i,j)*ye(i,j)
          beta(i,j)=xz(i,j)*xe(i,j)+yz(i,j)*ye(i,j)
          gama(i,j)=xz(i,j)*xz(i,j)+yz(i,j)*yz(i,j)
          binj(i,j)=xz(i,j)*ye(i,j)-xe(i,j)*yz(i,j)
    end do
  end do
  write(11,*)xz
  write(11,*)xzz
  write(11,*)xe
  write(11,*)xee
  write(11,*)yz
  write(11,*)yzz
  write(11,*)ye
  write(11,*)yee
  write(11,*)xze
  write(11,*)yze

  write(11,*)arfa
  write(11,*)beta
  write(11,*)gama
  write(11,*)binj

end program

! function subroutine to calculation the contraction factor
! 计算收缩因子的子程序
real(8) function con(i,j1,j2)  ! 采用Thompson建议形式,j1,j2分别为求和的起点和终点
  implicit none
  integer,intent(in)::i,j1,j2
  integer::j,m,l
  real(8)::ak(j1:j2),ck(j1:j2),con1,c1,c2,sig
  con=0.d0
  c1=100.d0 ! a=100.
  c2=0.5d0  ! b=0.5,详见陶文铨《数值传热学》第二版第448页介绍
  do j=j1,j2
    ak(j)=c1-(j-j1)*1.d0
    ck(j)=c2
  end do
  do j=j1,j2
    m=i-j
!         if(m)1,2,3
! 1   sig=-1.
!     l=-m
!     go to 4
! 3   sig=1.
!     l=m
! 4   gp=ck(j)*l
!     xm=dexp(gp)
!     con1=(-1.*ak(j)*sig)/(xm)
!     go to 6
! 2   con1=0.
! 6   con=con+con1
    if(m==0)then  ! 这是我改写的标准FORTRAN,以上是原文
      con1=0.d0
        else
          if(m<0)then
            sig=-1.d0
                l=-m
          else
            sig=1.d0
                l=m
          end if
      con1=(-1.d0*ak(j)*sig)/exp(ck(j)*l)
        end if
    con=con+con1
  end do
  return
end function con

[ 本帖最后由 weixing1531 于 2009-4-27 21:19 编辑 ]

ak47long

呵呵 我市菜鸟 学习学习

ak47long

请问你提到的计算流体力学基础是苏铭德的那本么？  程序中有些参数我不知道啥意思 我去把那本书找来看看

ak47long

你提到的原书，原文分别是什么书什么文呢？  我在苏铭德的《计算力体力学基础》上没找到阿

weixing1531

原书指的是陈仁景的《湍流模型及有限分析法》

程序是在第156页，是在此基础上修改的

网格的初始坐标采用PERIC写的 grid.f 程序生成

《计算力体力学基础》是清华大学出版的，作者是  任玉新 陈海昕

ak47long

谢谢了 能加个QQ交流下么  我的QQ112165626

jsccl1988

过了一年了，才看到你的帖子 不知道你还打理这个页面不。我正在做这方面的事情，做了一点点，想请教你一些问题。