gvs文件编写规范说明

说明：本gvs文件适用于N-S方程中对流项的基于有限体积法中迎风格式思想的处理。

gvs文件基本格式如下：

defi

disp u v w p

coor x,y,z

coef un,vn,wn

eqno 12

node 10

id p 1 1 1 1 0 0 0 0 0 0

$c6 dimension xc(3),uc(3),fx1(3),fx2(3),fx3(3),fx4(3)

mate ek fx fy 1.0d-2 0.0d0 0.0d0

equation

{

#

//

}

end

FORT

//

gvs文件包含defi,disp,coor,coef,eqno,node,id,mate,equation，FORT等关键字段。其中defi,disp,coor,coef,mate，FORT等字段其意义以及填写规则与pde文件填写规则相同，详见有关文献。

eqno字段表示单元内的方程的数目，node字段表示单元所包含的节点数，id字段表示压力节点。

注意，在gvs文件中数组定义应放在id字段和mate字段之间，语法同pde文件填写规则。

equation关键字标志有操作的程序段的开始，end字段标志结束。程序语句每行不得超过160列，且有相应关系的一段语句一般要放在一对{}里，前以#标识。

如：

equation

{

#

$c6 xc=(x(1)+x(2)+x(3))/3

}

例外情况是紧跟在equation字段后面的没有循环操作的语句，可不用加{}。如：equation

$c6 xc=(x(1)+x(2)+x(3))/3

有循环操作的语句必须用{}装入，如：

{

@ xm x 1 2 3 2 3 1 3 1 2 #

$c6 xm[3]=(x[1]+x[2])/2 }

@后跟用于循环的数组的名字（注意，这些数组需在前面的数组定义中定义，且必须为一维数组），下面的3行数字表示循环时的轮换的脚标顺序，循环设定结束时空一行，后以#开始程序语句。在程序语句中须用[]表示轮换的脚标。

循环设定还可以如下表示： {

&u u v {

&i 1 2 3 @ &u flux r 1 2 3 2 3 1 3 1 2 #

equation &i &u

%[(1.-r[3])*flux[3]%]&u[2]-%[r[2]*flux[2]%]&u[3] = 0.0d0 } }

这里以&u 和&i 定义变量，运行时对其分量进行循环操作（需要提请用户注意的是&后只能跟一个字符的变量名）。如上所示的是有嵌套的两重循环，&u 的u ，v 分量在外层循环，&i 的1，2，3分量在内层循环。@ &u flux r 这一行的&u 在内层循环中不替换。

程序语句中“%[”和“%]”取消其默认的轮换约定而变成实际的“[”和“]”。

程序的基本流程如下：

1. 计算体单元体心的坐标，速度（3D ）或面单元面心的坐标，速度（2D ）

2. 计算单元各条边的中点坐标，速度

3. 计算体单元的各个面的面心的坐标，速度（仅3D 下需要计算）

4. 计算各个节点的控制容积的流量

5. 计算对应流量的迎风系数r （子程序samar ）

6. 形成节点方程

对流项i 节点的节点方程的形成原理可以简单解释如下： N-S 方程的简化形式如(a ）式所示：

??

???=??=?+??+?-??0)(u f

p u u u t

u

； （a ）

由变分原理可得：

??

???=??=??-??+??-??0),()

,(),(),(),(),(q u f p u u u t

u

υυυυυ； （b ） 考察（b ）式中的对流项),(υu u ??。对于形如),(υu b ??

的项在二维情况下有如下表

达：

∑?=Ω??=??N j i

dxdy u b u b 1)(),(υυ

ds u u r n b j i N i i j ij

ij i )))(1)(((1

--?≈∑∑?=Λ∈τυ

（c ）

其中ij r 表示相邻节点i,j 之间的迎风系数，u 为待求的速度场。关于(c)式的详细推导过程我们不作更详细的说明，可以证明在网格尺寸足够小并且保持弱锐剖分的话，由(c)式所得到的方程有收敛的解。而i 节点的节点方程（二维）由下式给出：

∑?Λ∈=--?N i j ij

j i ij ds u u r n b τ0)))(1)(((

(d ）

三维情况下形式与之相似，在此略去。

下面分别就三角形6节点单元、四边形9节点单元、四面体10节点单元和六面体27节点单元四种情况分别说明如下： 1． 三角形6节点单元

如图1所示，我们考察三角形单元123，其节点1、2、3如图逆时针布置。我们先来看怎样形成1节点的节点方程，其他节点的节点方程可通过轮换节点编号获得。对单元123

图1

而言，我们只须计算控制容积在该单元内的(d)式即可（因为形函数在其他单元的值为0）。如图所示，单元的重心为c ；i-j 边的中点编号由该边所对应的顶点编号而命名，，比如12节点形成的边的中点我们编号为m3，其中的3对应于边1-2所对应的顶点3；对于i 节点的控

2

制容积，通过ij 节点之间界面的流量编号我们也按照同样的约定来命名，如1节点的控制容积中通过m3-c 的流量flux3。这里需要特别说明的是流量的计算。对于i 节点的控制容积，考虑通过ij 节点之间的界面的流量flux ，我们总这样约定：flux 等于以ij 之间界面上的平均速度叉乘由i-j 边中点指向重心的矢量，并且我们总假定平均速度方向由i 指向j 。如果叉乘之后的方向按右手定则向上，则其符号为正，反之则为负。以1节点的控制容积为例，考虑通过12节点之间的界面的流量（即通过边m3-c 的流量）flux3， flux3等于以边m3-c 的平均速度叉乘由1-2边中点m3指向重心c 的矢量，在程序段中如下表示：flux[3]=0.5*((um[3]+uc)*(yc-ym[3])-(vm[3]+vc)*(xc-xm[3]))。

迎风方案有这样的特点：1=+ji ij r r ， （1） 其中ij r 为i 节点的控制容积中ij 节点之间界面的迎风系数，它由单元内通过ij 节点之间界面的流量决定，在程序段中以子程序给出： subroutine samar(rflux,ek,rr) implicit real*8(a-h,o-z) if(rflux.ge.0.) then

rr=(ek+rflux)/(2.*ek+rflux) else

rr=ek/(2.*ek-rflux) end if end

其他三种单元类型下的迎风系数都按这个子程序计算得到。（1）式在任何情况下均成立．所以根据这一性质在三角形单元的程序编制中可以做适当的简化来形成最后的节点方程。 1节点的控制容积各个界面上的流量和迎风系数都计算出来之后，由(d)式便可分别得到1节点关于u 、v 的节点方程。

其他节点的节点方程只须轮换一下节点编号便可以同理得到。

2． 四边形9节点单元

图2

四边形单元的节点方程形成过程与三角形基本一样，不同的地方在于边与顶点不存在三角形中的对应关系。以图2所示的四边形单元为例，重心为c ，i-j 边的中点编号由该边的起始顶点编号而命名，如1-2边的中点编号m1。对于i 节点的控制容积（图2中没有全部画出），通过i 节点与j 节点之间的界面的流量编号也按照同样的规则命名，如通过1节点与2节点之间的界面的流量编号flux1。其他规则与三角形单元相同。

2

3．四面体10节点单元

图3

为了形成单元1234的1节点的节点方程，我们只须计算1节点的控制容积在该单元内的(d)式即可（因为形函数在其他单元的值为0）。如图所示，单元的重心为c ；i-j 边的中点编号由该边所对如图3所示，单元的重心为c ；i-j 边的中点编号由该边的两个顶点编号而命名，，比如12节点形成的边的中点我们编号为m12；面心的编号按面所对应的顶点编号来命名，如图中面124的重心编号f3；对于i 节点的控制容积，通过ij 节点之间界面的流量编号我们也按照ij 来命名，如1节点的控制容积中通过面m12-f4-c-f3的流量为flux12。对于i 节点的控制容积，考虑通过ij 节点之间的界面的流量flux_ij ，ij 节点之间的界面可看成两个三角形拼接而成，如图中1节点的控制容积12之间的界面看成由三角形f4-m12-c 和三角形c-m12-f3拼接而成；我们总这样约定：flux_ij 等于以ij 之间形成界面的两个三角形各自的平均速度点乘各自三角形的面积矢量获得（这里我们对平均速度的方向不做假定），对于面积矢量我们作出如下规定：按右手规则，大拇指由i 指向j ，在ij 之间界面上的三条分别以边中点为起点的矢量按手指的绕向逆时针安排，三角形也逆时针安排，其面积矢量为相应的两个矢量逆时针叉乘获得，以图3所示的1节点的控制容积为例，在12之间界面上的三个矢量逆时针安排为m12-f4、m12-c 、m12-f3，然后三角形f4-m12-c 的面积矢量为m12-f4?m12-c ，三角形c-m12-f3的面积矢量为m12-c ?m12-f3。上面流量的计算在程序段中如下表示： flux=ufu1(1)*(aa(2)*bb(3)-aa(3)*bb(2)) & -ufu1(2)*(aa(1)*bb(3)-aa(3)*bb(1)) & +ufu1(3)*(aa(1)*bb(2)-aa(2)*bb(1)) & +ufu2(1)*(bb(2)*cc(3)-bb(3)*cc(2)) & -ufu2(2)*(bb(1)*cc(3)-bb(3)*cc(1)) & +ufu2(3)*(bb(1)*cc(2)-bb(2)*cc(1)) flux=0.5*flux

其中aa 、bb 、cc 分别为逆时针排列的矢量，而ufu1、ufu2分别为两个三角形的平均速度（平均速度也可以该边中点的速度来代替，这种近似带来的影响并不十分显著，由用户考虑决定）。

1节点的控制容积各个界面上的流量和迎风系数（按前述子程序计算）都计算出来之后，由(d)式便可分别得到1节点关于u 、v 、w 的节点方程。

其他节点的节点方程只须轮换一下节点编号便可以同理得到。

４．六面体27节点单元

3

1

图4

六面体的计算与四面体基本相同，不同的地方在于面心的编号规则不一样。如图4所示，重心为c ，面心的编号按面的对角线上的顶点取得，如面1234的重心编号为f13，面1265的重心编号为f16。其他的规则与四面体相同。

1

2

3

5

7

m12

m15