想用2d三角形网格算一个台阶流动

shirazbj

我从http://nht.xjtu.edu.cn/nht/downpage.asp下载了delaunay.for & unstruct.for 程序。

用vba在word里编了一个画网格的程序。这个程序需要grid.dat保存在相同的目录下（和这个word文件）

生成了我的台阶流场。参看上传的附件。

那个顶盖流的例子，user.for里有这一段：
     DO I=NS-NB+51,NS-NB+100
      U(I)=1.0
      END DO
给顶盖速度赋初值。

不知道我的边界点怎么找？

[ 本帖最后由 shirazbj 于 2010-6-12 18:22 编辑 ]

shirazbj

打开该word文件
选tools\marco\visual basic editor打开vba程序
把光标点入源程序区
在工具栏里点三角执行图标执行该程序

执行结果网格绘制在文档里

[ 本帖最后由 shirazbj 于 2010-1-6 10:34 编辑 ]

shirazbj

我把那个三角形化的程序改成了VB的。想画出三角形化的过程，图视看看具体的过程。

在把每个entry变成sub时，发现一个kiss变量，在ENTRY BACKGROUND里这样
1400  CONTINUE
      KISS=NS
      RETURN

在ENTRY SORTING里，这样：

DO 1500 I=NS-NADD+1,NS
      NCIRCLE=I
      CALL CIRCLE
      XT=X0
      YT=Y0
      DO 1600 J=1,KISS

变成sub后，因为KISS没在公共块里，所以值传不过来。

问题是在fortran里没问题么？fortran好像是在第1次使用时就相当于声明了。对于不同的entry,可以直接使用在另外的entry里用到的变量？

这里（ENTRY INITIAL）
BACKGROUND是先被调用的，给kiss赋了值，然后是SORTING，所以没问题。BACKGROUND也写在SORTING的前面。这是没问题的原因么？

shirazbj

主程序里有这么一段初始化：
DO I=1,NS
      U(I)=0.
      V(I)=0.
     。。。  
END DO
NS是三角形的数量。
user.for里
DO I=NS-NB+51,NS-NB+100
nb是边界点的个数。
51点到100点是顶盖边界。

所以要自己数那条边对应那些点。

把我的台阶网格数好，改好user.for,开始计算，有个数学错误，除数为零。

发现生网格程序生成的网格文件和计算程序用的网格文件结构不一样。前者多了4个包容矩形的点。计算时，FIND NODE'S NEIBOR CENTER POINTER时，返回0。

shirazbj

编了一个小程序把grid.dat里的4个包容矩形点踢掉，并把受影响三角形的点数减4。终于可以算了。但是又来了一个溢出错误。Woo-la-la.

shirazbj

今天把这个程序在Linux下编译了一遍
在Ubuntu8.04下,还是用的g95编译器,参看http://www.g95.org/
适用Ubuntu的包是在这下的:http://www.gfd-dennou.org/library/cc-env/g95/g95/etch/g95_0.91_i386.deb,
把该文件保存在自己的目录后，浏览到该文件,右键菜单选Gdebi安装即可.
在终端里键入:
g95 2008del.for   
其中2008del.for是你的程序名，
就编译了,生成a.out运行文件.
执行程序在终端里键入:
./a.out
要注意linux的文件名区分大小写。

shirazbj

加了一个输出gmsh网格的entry如下：

附图只是一个插入10个内点例子。


        ENTRY GMSH
      OPEN(8,FILE='GMSH.MSH')
        WRITE(8,*) '$MeshFormat'
        WRITE(8,*) '2 0 8'
        WRITE(8,*) '$EndMeshFormat'
        WRITE(8,*) '$PhysicalNames'
        WRITE(8,*) '1'
        WRITE(8,*) '2 100 "My fancy surface label" '
        WRITE(8,*) '$EndPhysicalNames'
        WRITE(8,*) '$Nodes'
        WRITE(8,*) NPP
      DO 3800 I=1,NPP
      WRITE(8,*)I,APOINT(I,1),APOINT(I,2),0
3800  CONTINUE
        WRITE(8,*) '$EndNodes'
        WRITE(8,*) '$Elements'
        WRITE(8,*) NS
      DO 3900 I=1,NS
      WRITE(8,*)I," 2 3 0 6 0",KTSTACK(I,1),KTSTACK(I,2),KTSTACK(I,3)
3900  CONTINUE
        WRITE(8,*) '$EndElements'
      CLOSE(8)
        RETURN

shirazbj

开源的gmsh在这下载http://geuz.org/gmsh/

在unstruct的user里加了个输出矢量图的entry

附图是顶盖流的结果

      ENTRY GMSHvect
        OPEN (3,FILE='GMSHFEILD.POS')          ! gmsh RESULT
        WRITE(3,*)'View "A vector map" {'
        DO I=1,NS-NB
        WRITE(3,*)'VP(',CENTER(I,1),',',CENTER(I,2),',0){'
        &,U(I),',',V(I),',',0,'};'
        ENDDO
        WRITE(3,*)'};'
        CLOSE(3)
      RETURN

[ 本帖最后由 shirazbj 于 2010-5-17 08:05 编辑 ]

qxzzcc

shirazbj

终于得到了个收敛解，方法是减小入口的速度。

我的台阶是10宽1高，台阶在中间5处，0.5高。原来缺省的入口速度是1，如果是空气的话，取密度1.205，特征长度1，室温粘度0.00001983，算出
re = 1.205 * 1 * 0.5 / 0.00001983=30000多，这个程序没湍流模型，应该是算不了。

另一个问题，原程序算的是顶盖驱动流，上边一个u速度，其他3边固壁条件。我的台阶进口好办，u=0.01,v=0,但出口想只给v=0.  u应该从内节点推出，这要改程序了。麻烦。

weixing1531

楼主强

现在我连那个程序都还没有看懂

shirazbj

写了一小段程序，把每次迭代的结果存在了一个文件里，然后用后处理软件看计算过程的图像。附件是矢量图的例子，只存了100来步。这个lid driven例子算到800多，图像就变化不大了。因为程序只有cell的数据，简单地把它赋到节点上了。准确的得再加段插值的程序。

shirazbj

这是没改出口条件的图示。流体最终流不出去。

shirazbj

终于改好了出口条件。

用了个笨办法找到和出口三角形有限体共边的内三角形有限体，把内三角形的U值简单赋给出口有限体。代码如下，加在计算程序添加边界有限体代码的后边.

*-----find outlet's neighbor cell------------------------------------------
      DO I=210,210+10
           do J=1,NS
            N1=TRIANGLE(j,1)
            N2=TRIANGLE(j,2)
            N3=TRIANGLE(j,3)
              if (N1.eq.I) then
                 if ((N2.eq.(I+1)).or.(N3.eq.(I+1))) then
                    NBC(NS+I)=J
                 endif
                  end if
              if (N2.eq.I) then
                 if ((N1.eq.(I+1)).or.(N3.eq.(I+1))) then
                    NBC(NS+I)=J
                 endif
                  end if
              if (N3.eq.I) then
                 if ((N1.eq.(I+1)).or.(N2.eq.(I+1))) then
                    NBC(NS+I)=J
                 endif
                  end if
           end do
           write(*,*) I,NS+I,NBC(NS+I)
      END DO

在每次计算后,更新出口的U值.

*-----CORRECT outlet--------------------
      DO I=NS-NB+210,NS-NB+210+10
           U(I)=U(NBC(I))
        end do
加在这句的前面.
      ITER=ITER+1

这些代码是网格相关的, 象NS-NB+210中的210就是出口边界的起点.所以要算别的网格,要自己改.

附图是出口局部处迭代的图示. 可以看到进口边界速度通过迭代传过来的过程.因为gif文件大小的限制，只录了前面的迭代，还没到完全收敛。粗一看，出口的速度抛面形状还可以，应该还是层流流动吧。内节点还应该加密点。

贝贝daisy

为什么当确定某点在此三角形内，却要：
IF(LT) THEN
      MPOINT(1)=N1
      MPOINT(2)=N2
      MPOINT(3)=N3
      NPOINT=3
      DO 3200 JJ=J,NS-1
      DO 3200 JM=1,3
      KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM) ？？？？？？？？？？将第二个点的坐标赋予第一个点？
3200  B(JJ)=B(JJ+1)
      NS=NS-1
      GOTO 3300
      END IF