求助_涡流函数法

chtq22 · 发表于 2009-6-28 16:48:12

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？注册

x

小弟诚求涡流函数形式求解N-S方程的程序，欲计算圆柱绕流问题，与我现在的方法进行比较

liuhuafeifei · 发表于 2009-7-2 14:55:39

program cavity_sv
c
c This program is to be used only under the
c stipulations of the licensing agreement.
c----------------------------------------

c--------------------------------------------------
c  Steady flow in a two-dimension cavity driven
c  by a moving lid, computed by the
c  stream function / vorticity formulation
c
c  The cavity extends over 0<x<Lx, 0<y<Ly
c
c  SYMBOLS:
c  -------
c
c  Lx, Ly: cavity size
c  Nx, Ny: Number of divisions in the x and y directions
c
c  X,Y:  grid-point coordinates
c  U,V:  velocity components
c  S: streamfunction (psi)
c  O:    vorticity (omega)
c
c--------------------------------------------------

   Implicit Double Precision (a-h,o-z)

   Double Precision Lx,Ly

   Dimension X(129,129),Y(129,129)
   Dimension U(129,129),V(129,129)

   Dimension S(129,129),Snew(129,129)

   Dimension O(129,129),Onew(129,129)

   Dimension Source(129,129)

   Real vector(5,2)  ! to draw vectors

c----------
c constants
c----------

   pi  = 3.14159265358 D0
   pi4 = 4.0D0*pi
   pis = pi**2

c-----------
c parameters
c-----------

c    write (6,*)
c    write (6,*) " Choose the data input mode"
c    write (6,*)
c    write (6,*) " Enter 0 to quit"
c    write (6,*) "    1 to read data from file: cvt_sv.dat"
c    write (6,*) "    2 to type in data"
c    write (6,*) " ------------------------------------------"
c    read  (5,*) Ienrd

   Ienrd = 1

   If(Ienrd.eq.0) Go to 99

c------------------------
   If(Ienrd.eq.2) then
c------------------------

   write (6,*) " Please enter the lid velocity"
   write (6,*) " -----------------------------"
   read  (5,*) UUlid

   write (6,*) " Please enter cavity dimensions Lx and Ly "
   write (6,*) " ---------------------------------------- "
   read  (5,*) Lx,Ly

   write (6,*) " Enter discretization levels Nx and Ny "
   write (6,*) " --------------------------------------"
   read  (5,*) Nx,Ny

   write (6,*) " Enter the viscosity "
   write (6,*) " ------------------- "
   read  (5,*) visc

   write (6,*) " Enter the density "
   write (6,*) " ----------------- "
   read  (5,*) rho

   write (6,*) " Enter the relaxation parameters "
   write (6,*) " relax(psi), relax(vort) "
   write (6,*)
   write (6,*) " Optimal values are about 0.10, 0.10 "
   write (6,*) " ----------------------------------- "
   read  (5,*) rel_psi,rel_vort

   write (6,*) " How many local iterations "
   write (6,*) " for the stream function and the vorticity ? "
   write (6,*)
   write (6,*) " Better to choose less than 3 "
   write (6,*) " ---------------------------- "
   read  (5,*) Niter

   write (6,*) " How many global iterations before pausing ? "
   write (6,*) " ------------------------------------------ "
   read  (5,*) Niterg

   write (6,*) " Enter the inital guess for the vorticity "
   write (6,*) " at the grid points"
   write (6,*) " ------------------------------------ "
   read  (5,*) Vort_init

   write (6,*) " See transient streamfunction on screen ? "
   write (6,*)
   write (6,*) " Enter 0 for no; 1 for yes"
   write (6,*) " -------------------------"
   read  (5,*) IseeS

   write (6,*)
   write (6,*) " See transient vorticity on screen ? "
   write (6,*)
   write (6,*) " Enter 0 for no; 1 for yes"
   write (6,*) " -------------------------"
   read  (5,*) IseeV

c----------
   Else
c----------

      open (1,file="cvt_sv.dat")

      read (1,*) UUlid
      read (1,*)
      read (1,*) Lx,Ly
      read (1,*) Nx,Ny
      read (1,*)
      read (1,*) visc
      read (1,*) rho
      read (1,*)
      read (1,*) rel_psi,rel_vort
      read (1,*) Niter
      read (1,*) Niterg
      read (1,*) Vort_init
      read (1,*)
      read (1,*) IseeS
      read (1,*) IseeV

      close (1)

c-----------
   End If
c-----------

c----------
c constants
c----------

   Null  = 0
   None  = 1
   Ntwo  = 2
   Nfive = 5

   zero = 0.0D0

c--------
c prepare
c--------

   Nx1 = Nx+1
   Ny1 = Ny+1

   Dx = Lx/(Nx1-1.0D0)
   Dy = Ly/(Ny1-1.0D0)

   Dx2 = 2.0D0*Dx
   Dy2 = 2.0D0*Dy

   Dxs  = Dx**2
   Dys  = Dy**2

   beta = Dxs/Dys

   rnu = visc/rho  ! kinematic viscosity

c-----------------------------------------
c Initialize vorticity and stream function
c-----------------------------------------

   Do I=1,Nx1
   Do J=1,Ny1
      X(I,J) = (I-1.0D0)*Dx
      Y(I,J) = (J-1.0D0)*Dy
      S(I,J) = 0.0D0
      O(I,J) = - Vort_init
   End Do
   End Do

c-------------------------------
c Beginning of global iterations
c-------------------------------

   Kter = 0    ! global iterations counter
   Jter = 0    ! batch iterations counter

98 Continue

c-----------------------------------
c Iterate on Poisson equation
c for the stream function
c with Dirichlet boundary conditions
c-----------------------------------

   Iter = 0    ! stream function iteration counter

  93  Continue

   Iter = Iter + 1

c---------------------------------------
c Jacobi updating of the stream function
c at the interior nodes
c---------------------------------------

   Errp = 0.0 ! pointwise error register

   Do J=2,Ny
   Do I=2,Nx

   res =       S(I+1,J)+S(I-1,J)-2.0*S(I,J)
   +    + beta*(S(I,J+1)+S(I,J-1)-2.0*S(I,J))
   +    +    Dx**2 * O(I,J)

   Snew(I,J) = S(I,J) + rel_psi * res

   If(abs(res).gt.ERRP) Errp = abs(res)

   End Do
   End Do

c-------
c update
c-------

   Do J=2,Ny
   Do I=2,Nx
      S(I,J) = Snew(I,J)
   End Do
   End Do

   If(Iter.le.Niter) Go to 93

c-----------------
c printing session
c-----------------

   If(IseeS.eq.1) then

      write (6,*)
      write (6,*) " STREAM FUNCTION "
      write (6,*)

      Do j=1,Ny+1
      write (6,101) (S(i,j),i=1,Nx+1)
      End Do

   End If

c-----------------------------
c Compute the vorticity
c at the internal grid points
c----------------------------

c    Do I=2,Nx
c    Do J=2,Ny
c    O(I,J) = - (S(I+1,J)-2.0*S(I,J)+S(I-1,J))/Dxs
c +          - (S(I,J+1)-2.0*S(I,J)+S(I,J-1))/Dys
c    End Do
c    End Do

c-------------------------------------
c Compute the vorticity at boundary grid points
c using the velocity boundary conditions
c
c Lower-order boundary conditions are commented out
c-------------------------------------

c---
c top and bottom walls
c---

   Do I=2,Nx

c    O(I,1  ) = 2.0*(S(I,1)  -S(I,2)) /Dys
c    O(I,Ny1) = 2.0*(S(I,Ny1)-S(I,Ny))/Dys - 2.0*UUlid/Dys

      O(I,1  ) = (7.0*S(I,  1)-8.0*S(I,2) +S(I,  3))/(2.0*Dys)
      O(I,Ny1) = (7.0*S(I,Ny1)-8.0*S(I,Ny)+S(I,Ny-1))/(2.0*Dys)
   +          - 3.0*UUlid/Dy

   End Do

liuhuafeifei · 发表于 2009-7-2 14:56:05

c---
c left and right walls
c---

   Do J=2,Ny

c    O(1  ,J) = 2.0*(S(1,J) - S(2,J) )/Dxs
c    O(Nx1,J) = 2.0*(S(Nx+1,J) - S(Nx,J))/Dxs

      O(1  ,J) = (7.0*S(1  ,J)-8.0*S(2,J) +S(3, J))/(2.0*Dxs)
      O(Nx1,J) = (7.0*S(Nx1,J)-8.0*S(Nx,J)+S(Nx-1,J))/(2.0*Dxs)

   End Do

c--------------------------------
c Compute the velocity at the grid points
c       by central differences
c--------------------------------

   Do J=2,Ny
   Do I=2,Nx
      U(I,J) = (S(I,J+1)-S(I,J-1))/Dy2
      V(I,J) = - (S(I+1,J)-S(I-1,J))/Dx2
   End Do
   End Do

c--------------------------------------------
c Iterate on Poisson's equation for the vorticity
c
c The source term is the rhs of equation (13.1.8)
c of Pozrikidis (1997)
c--------------------------------------------

   Iter = 0

  92  Continue

   Iter = Iter + 1

c---
c Jacobi updating
c---

   ErrO = 0.0D0 ! error register

   Do J=2,Ny
   Do I=2,Nx

   source(I,J) =  U(I,J)*(O(I+1,J)-O(I-1,J))/Dx2
   +             + V(I,J)*(O(I,J+1)-O(I,J-1))/Dy2

   source(I,J) = source(I,J)/rnu

   res =       O(I+1,J)+O(I-1,J)-2.0*O(I,J)
   +    + beta*(O(I,J+1)+O(I,J-1)-2.0*O(I,J))
   +    -    Dxs * source(I,J)

   Onew(I,J) = O(I,J) + rel_vort * res

   If(abs(res).GT.ErrO) ErrO = abs(res)

   End Do
   End Do

c---
c replace
c---

   Do J=2,Ny
   Do I=2,Nx
      O(I,J) = Onew(I,J)
   End Do
   End Do

   If(Iter.le.Niter) Go to 92

c-----------------
c printing session
c-----------------

   If(IseeV.eq.1) then

      write (6,*)
      write (6,*) " Vorticity"
      write (6,*)

      Do J=1,Ny1
      write (6,101) (O(I,J),I=1,Nx1)
      End Do

   End If

   write (6,103) Kter,ERRP,ERRO

c-----------------
c advance counters
c-----------------

   Jter = Jter + 1
   Kter = Kter + 1

   If(Jter.le.Niterg) Go to 98

c------------------------
c End of batch iterations
c------------------------

   Jter = 1

   write (6,*)
   write (6,*) " One more set of global iterations ?"
   write (6,*)
   write (6,*) " Enter 0 for No, 1 for YES "
   write (6,*) " ---------------------------"
   read  (5,*) more

   If(more.ne.0) Go to 98

c-------------------------
c-------------------------
c End of global iterations
c-------------------------
c-------------------------

99 Continue

c--------------------------------
c print the final stream function
c--------------------------------

   open (4,file="cvt_sv.str")

      write (4,*) Nx,Ny
      write (4,*) Dx,Dy

      Do j=1,Ny1
      write (4,101) (S(i,j),i=1,Nx1)
      End Do

      write (4,*)
      write (4,400) UUlid
      write (4,401) Lx,Ly
      write (4,402) Nx,Ny
      write (4,403) visc
      write (4,404) rho
      write (4,405) rel_psi
      write (4,406) rel_vort
      write (4,407) Niter
      write (4,408) Niterg

   close (4)

c--------------------------------
c print the final vorticity field
c--------------------------------

   open (4,file="cvt_sv.vort")

      write (4,*) Nx,Ny
      write (4,*) Dx,Dy

      Do j=1,Ny1
      write (4,101) (O(i,j),i=1,Nx1)
      End Do

      write (4,*)
      write (4,400) UUlid
      write (4,401) Lx,Ly
      write (4,402) Nx,Ny
      write (4,403) visc
      write (4,404) rho
      write (4,405) rel_psi
      write (4,406) rel_vort
      write (4,407) Niter
      write (4,408) Niterg

   close (4)

c--------------------------------------------------
c print the final upper wall vorticity distribution
c along with the local solution for Stokes flow
c--------------------------------------------------

   open (4,file="cvt_sv.vortw")

      write (6,*) " upper wall vorticity distribution"
      write (6,*)

      j=Ny1
      write (4,*) Nx1

      Do i=1,Nx1
      xh = (i-1.0D0)/Nx
      Ocorner = 0.0D0

      If(xh.ne.0) then
      Ocorner = - pi4/(pis-4.0D0) * UUlid/xh
      End If

      write (4,102) i,xh,O(i,j),Ocorner
      write (6,102) i,xh,O(i,j),Ocorner
      End Do

      write (4,*) Null

   close (4)

c-----------------------------
c Plot velocity vectors
c reduced by maximum magnitude
c for visualization
c-----------------------------

   open (4,file="cvt_sv.vel")

c---
c print cavity walls
c---

   write (4,104) Nfive
   write (4,104) None,None,zero,zero
   write (4,104) None,None,Lx  ,zero
   write (4,104) None,None,Lx  ,Ly
   write (4,104) None,None,zero,Ly
   write (4,104) None,None,zero,zero

c-----
c Plot velocity vectors
c reduced by maximum magnitude
c for visualization
c----

   Umax = 0.0D0

   Do J=2,Ny
   Do I=2,Nx
      Umag = U(I,J)**2+V(I,J)**2
      If(Umag.gt.Umax) Umax = Umag
   End Do
   End Do

   Umax = Dsqrt(Umax)

   write (6,*) " Maximum velocity: ", Umax

   Do J=2,Ny
   Do I=2,Nx

      x1 = X(I,J)
      y1 = Y(I,J)
      Darx = 0.10*U(I,J)/Umax
      Dary = 0.10*V(I,J)/Umax

      call draw_arrow_2d
   +
   + (x1,y1
   + ,Darx,Dary
   + ,vector
   + )

      write (4,104) Nfive

      Do l=1,5
      write (4,104) J,I,vector(l,1),vector(l,2)
      End Do

   End Do
   End Do

   write (4,104) Null

   write (4,*)
   write (4,400) UUlid
   write (4,401) Lx,Ly
   write (4,402) Nx,Ny
   write (4,403) visc
   write (4,404) rho
   write (4,405) rel_psi
   write (4,406) rel_vort
   write (4,407) Niter

   close (4)

c-----
c Done
c-----

  100 Format (2(1x,f15.10))
  101 Format (200(1x,f10.5))
  102 Format (1X,I4,3(1x,f15.8))
  103 Format ("Iter: ",i4,", Max error in psi and omega ",
   +    2(1x,f15.10))
  104 Format (2(1x,i2),20(1x,f10.5))

400  Format (1x," UUlid    ",F10.5)
401  Format (1x," Lx,Ly    ",F10.5,1x,F10.5)
402  Format (1x," Nx,Ny    ",I3,1x,I3)
403  Format (1x," Viscosity ",f10.8)
404  Format (1x," Density    ",f10.8)
405  Format (1x," rel_psi    ",f10.8)
406  Format (1x," rel_vort ",f10.8)
407  Format (1x," Niter    ",I5)
408  Format (1x," Niterg    ",I5)

   Stop
   End

liuhuafeifei · 发表于 2009-7-2 14:56:55

1.0  UUlid

1.0 0.5    Lx, Ly
32 16 64 32 128 64 128 128       Nx, Ny

0.01       viscosity
1.00       density

0.1 0.1 0.20 0.20 relaxation factors for psi and omega
05 10 03 02       number of inner iterations
100       number of global batch iterations
0.0       initial vorticity

0          IseeS
0          IseeV

chtq22 · 发表于 2009-7-2 21:52:09

非常感谢，要是能有圆柱绕流的例子就更好了

liuhuafeifei · 发表于 2009-7-3 08:30:41

兄弟，你也太懒了，涡量流函数法是CFD中比较简单的一种，写程序没几句。
直交坐标改柱坐标，很直接，一天估计就搞定，不就是要考虑r吗

liuhuafeifei · 发表于 2009-7-3 08:39:18

哦，不是柱坐标，估计需要适体坐标。还又可能涉及多块，看你怎么作了。
有点难度。参考一下。

不过我以前有适体的单块涡量流函数的程序，回去找找，好多年了。

onesupeng · 发表于 2009-7-3 08:45:14

真的有那么复杂吗？

onesupeng · 发表于 2009-7-3 08:47:43

原帖由 chtq22 于 2009-6-28 08:48 发表

 登录/注册后可看大图

小弟诚求涡流函数形式求解N-S方程的程序，欲计算圆柱绕流问题，与我现在的方法进行比较

两种途径，一是直接找参考文献就是了，圆柱绕流涡量流函数算的，这会少么？
二是自己写一段

liuhuafeifei · 发表于 2009-7-3 10:31:10

先想想，常规的坐标系能处理吗（直角，柱，球）？
如要分析分离等现象，单块结构化网格可以吗？

liuhuafeifei · 发表于 2009-7-3 10:57:28

随便看了一篇，是单块直交坐标，有限差分，圆柱边界处采用阶梯网格，陶老师的红皮书上介绍过，还算简单，将上面程序修改一下，细心一点即可。哎，很古老的方法，近些年都搞非结构的，差点忘了。对付这个问题可以。
另一种是多重坐标，太繁，不推荐。

chtq22 · 发表于 2009-7-3 20:22:43

我看到的文献多是用极坐标，这样可以描述圆柱的曲线边界，然后经坐标转换，变成直角坐标，利用差分法求解。
其实，问题是这样的，我发现两篇文献用涡流函数形式的N-S方程所得时均升力为正值，而其它方法结果为负值。其所模拟流动Re=50-150，属于层流状态，涡已经脱离。
我用LBM计算，时均升力也是负值，顾想找个涡流函数形式程序，比较一下。

onesupeng · 发表于 2009-7-4 11:51:46

这个肯定地说，没有意义

chtq22 · 发表于 2009-7-4 20:33:21

onesupeng君为什么说这没有意义呢，请指教。
我是想看看究竟是边界条件，阻塞比，网格粗细等引起不同学者的结果不一样，还是其它什么原因引起的，按理说升力均值的方向在特定Re下应该是一定的。

onesupeng · 发表于 2009-7-4 21:45:01

我说的就是，你说得上面讲要研究的这些，是没有意义的。按理说，应该是零，所谓正负由于数值误差和不同初值条件引起的而已。所有你提到的，肯定有非常多的人做过，高兴多调研一下就好了。

类似网格什么的，任何一种算法的基本功而已，不足以当个话题研究什么的

		自动登录	找回密码
密码			注册

求助_涡流函数法

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

p1

p2

dat

回复 4# liuhuafeifei 的帖子

re

回复 11# liuhuafeifei 的帖子

回复 13# onesupeng 的帖子

站点信息