网格生成技术之无限插值法

合集下载

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

¤Grid Generation Series¤
网格生成•••适体坐标系•••代数方法•••无限插值法
Copyright © 2007 版权所有
目录
1. 概述 (1)
1.1前序 (1)
1.2名词解释 (1)
1.3映射关系 (1)
2. 二维无限插值法生成网格 (2)
2.1模型公式 (2)
2.2操作步骤 (3)
2.3编程实例 (3)
3. 三维无限插值法生成网格 (9)
3.1计算公式 (9)
3.2编程实例 (11)
4. 参考文献 (16)
5. 版权声明 (17)
1. 概述
1.1前序
网格生成技术的编程实现与应用曾是笔者感觉深奥而有趣的事情。

因对这方面并不熟悉，2006年夏初，笔者决定做网格生成方面的努力。

查看了若干资料后，虽对其数学原理未有涉足，但还是有幸获得或推出了生成二维网格和三维网格生成的TFI实现公式，并编写了测试程序进行验证，如封面图片的网格就是当时笔者用TFI方法编程生成的。

后来忙别的事情，就一直落在“纸堆”里。

2007年夏初，有网友询问TFI，又想起来，于是四处找了找，看着笔记，发现一年前的清晰思路都模糊了。

当时在图书馆借过一本书对我的帮助也很大，但书名已记不起来了，无法在后面的参考文献中列出。

现在决定用休息时间把笔记整理一下，以供需要的朋友查阅。

无限插值法（TFI）是结构化网格生成技术中属于适体坐标系的代数方法。

其优点是算法简单、生成网格速度很快，对于较规则区域，TFI法得到的网格效果也令人满意。

对于没有把握的复杂区域，笔者认为最好采用TFI方法生成初始网格场，然后采用PDE(偏微分网格生成技术)进行网格场优化。

1.2名词解释
（1）网格生成技术：对给定区域进行离散以生成计算网格的方法。

（2）结构化网格：排列有序、相邻节点位置关系明确的网格。

（3）适体坐标系：坐标轴与计算区域的边界一致的坐标系，又称贴体坐标系、附体坐标系。

（4）代数方法：通过代数关系式创建物理平面上的区域与计算平面上的区域的映射方法。

（5）无限插值法：把边界上规定的对应关系连续插值到区域内部，插值的点数是无限的，因而称为无限插值（transfinite interpolation，TFI）。

（6）物理空间：真实的求解区域。

通常不规则，不易进行网格节点剖分计算。

（7）计算空间：进行网格节点剖分的区域。

规则，最常见的为矩形区域或长方体，网格节点定位计算简单。

1.3映射关系
在计算空间内剖分得到的节点需要映射回物理空间，以便于进行物理求解。

适体坐标的网格生成方法的核心就在于，给出从计算空间到物理空间，节点位置的数学映射关系。

2. 二维无限插值法生成网格
2.1模型公式
这里以二维TFI方法为例进行说明。

参阅的资料中，都习惯以把计算区域“单位化”到[0,1]范围内，本文也遵守该约定，因此(ξ,η)的取值范围也在[0,1]内。

要通过计算空间的各网格节点位置(ξ,η)获得物理空间上对应的各网格节点的位置（x，y），相当于把x、y均看作变量ξ和η的函数，即x (ξ,η)，y (ξ,η)。

如何求得呢？
根据文献得到数学求解公式：
（Pξ⊕Pη）[x]=Pξ[x]+Pη[x]-PξPη[x]
（Pξ⊕Pη）[y]=Pξ[y]+Pη[y]-PξPη[y]
笔者对上述符号公式进行了推导（推导过程略），可整理得到如下具体的计算公式：
x(ξ,η)=(1-η)x b(ξ)+ηx t(ξ)+(1-ξ)x
(η)+ξx r(η)
-[ξηx t(1)+ξ(1-η)x b(1)+η(1-ξ)x t(0)+(1-ξ)(1-η)x b(0)]
y(ξ,η)=(1-η)y b(ξ)+ηy t(ξ)+(1-ξ)y L(η)+ξy r(η)
-[ξηy t(1)+ξ(1-η)y b(1)+η(1-ξ)y t(0)+(1-ξ)(1-η)y b(0)]
式中，
（1）变量x、y均带有4个下标t、b、L、r，是指该点为二维物理空间的上下左右四条边界上的点。

如x b(ξ)表示在二维物理空间，对应于内部节点x(ξ,η) 的下边界点的横坐标值。

（2）x t(0)表示在二维物理空间，上边界最左侧点的横坐标值；y b(1)表示在二维物理空间，下边界最右侧点的纵坐标值，其它类似（略）。

可见，要得到二维物理空间内部节点划分，必须已知两个条件：
（1）计算空间中各个节点的位置(ξ,η)。

（2）物理空间中4条边界上所有节点的位置。

注：⊕符号，在文献中称为布尔和，在图形学中称为环和，不妨以概率中的事件和来理解。

2.2操作步骤
采用TFI方法生成网格的大致步骤如下：
（1）首先，根据需要（如采用等间距还是渐疏/渐密等策略）手工或编程生成二维物理空间求解区域的t、b、L、r（即上下左右）4条边界上的节点（位置，（x，y））。

（2）接着，根据物理空间的网格划分需求，确定节点划分策略（如矩形网格形式），手工或编程生成计算空间的各个节点（位置，（ξ,η））。

（3）采用上节给出的计算公式，注意节点对应关系，求出物理空间上内部区域的各个网格节点位置（x(ξ,η)，y(ξ,η)）。

（4）把物理空间的网格节点信息存储为需要的格式，进行显示验证或供物理求解调用。

2.3编程实例
程序需要读取物理平面的4条边界上的节点信息和计算平面的全部节点信息。

由于篇幅所限，不再贴上输入文件数据。

本例代码在Visual Fortran6.5编译环境下运行测试通过。

东岸线 2006-6-21 ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～
!///////////////////////////////////主程序////////////////////////////////////////////////
program main
integer i,j,imax,jmax
!*******打开物理平面边界坐标文件，先读取网格范围编号。

i从1～imax，j从1～jmax**** open(1,file='physicsb.dat',err=555)
read(1,*) imax,jmax
!************调用子程序**************************************
call compute(imax,jmax)
!************返回主程序**************************************
goto 888
555 print *,'ERROR OCCURED 读取文件出错！请检查无误后，重新进行计算。

'
pause
888 return
end subroutine
!////////////////////////////////////////////////////////////////////////////////////////////// ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～
注意：程序中，i、j为节点编号；x，y为物理平面节点坐标；xx，yy为计算平面节点坐标。

算例说明
1.本例考虑的二维环形物理空间：
如图示，共有1-1、1-2、2-2、2-1（即①②③④）四条边界。

2.对上图物理空间取如下对应的计算空间：
边界网格节点分布划分，采用均分正方形网格单元节点。

每行、每列均为10等份（刚发现示意图中少画了一行^_^!），即有11个节点（0、0.1…、0.9、1），
3.这是程序读入的计算空间所有节点坐标（xx，yy）的文件（也可在程序中生成）：
computec.dat
4.把物理空间①②③④四边也均分为10等份，为11个节点。

这是程序读入的物理空间4条边界上所有节点坐标（x，y）的文件：
physicsb.dat
5.这是程序生成的物理空间所有网格节点的坐标（x，y）文件：
2Dmesh.grd
6.这是2Dmesh.grd读入Gridgen的显示效果：（grid pts―import―domain－surfaces）
行文至此，2007-6-24下午17：05。

忽闻外界笛声嘈杂，方感身体微凉，略有不安，赶紧遥控关了空调。

砰然门窗翕动，接着巨响隆起，霎时楼外瓢泼如注，其势如…约半小时后，风停雨歇，雷声遁去…
3. 三维无限插值法生成网格
3.1计算公式
根据文献，得到数学求解公式：
(P
ξ⊕Pη⊕Pζ)[x]=Pξ[x]+Pη[x]+Pζ[x]-PξPη[x]-PηPζ[x]-PζPξ[x]+PξPηPζ[x]
(Pξ⊕Pη⊕Pζ)[y]=Pξ[y]+Pη[y]+Pζ[y]-PξPη[y]-PηPζ[y]-PζPξ[y]+PξPηPζ[y]
(Pξ⊕Pη⊕Pζ)[z]=Pξ[z]+Pη[z]+Pζ[z]-PξPη[z]-PηPζ[z]-PζPξ[x]+PξPηPζ[z] 笔者对上述符号公式进行了推导（推导过程略），可整理得到如下具体的计算公式：（这里给的计算公式正是下一小节程序的核心部分，这里直接使用了其中的变量）
px1=(1-xx(i,j,k))*x(1,j,k)+xx(i,j,k)*x(imax,j,k)
px2=(1-yy(i,j,k))*x(i,1,k)+yy(i,j,k)*x(i,jmax,k)
px3=(1-zz(i,j,k))*x(i,j,1)+zz(i,j,k)*x(i,j,kmax)
px1px2=(1-xx(i,j,k))*(1-yy(i,j,k))*x(1,1,k)+xx(i,j,k)*(1-yy(i,j,k))*x(imax,1,k)+yy(i,j,k)*(1-xx(i,j,k)) *x(1,jmax,k)+xx(i,j,k)*yy(i,j,k)*x(imax,jmax,k)
px2px3=(1-yy(i,j,k))*(1-zz(i,j,k))*x(i,1,1)+yy(i,j,k)*(1-zz(i,j,k))*x(i,jmax,1)+zz(i,j,k)*(1-yy(i,j,k))* x(i,1,kmax)+yy(i,j,k)*zz(i,j,k)*x(i,jmax,kmax)
px3px1=(1-zz(i,j,k))*(1-xx(i,j,k))*x(1,j,1)+zz(i,j,k)*(1-xx(i,j,k))*x(1,j,kmax)+xx(i,j,k)*(1-zz(i,j,k))* x(imax,j,1)+zz(i,j,k)*xx(i,j,k)*x(imax,j,kmax)
px2px3i1=(1-yy(i,j,k))*(1-zz(i,j,k))*x(1,1,1)+yy(i,j,k)*(1-zz(i,j,k))*x(1,jmax,1)+zz(i,j,k)*(1-yy(i,j,k) )*x(1,1,kmax)+yy(i,j,k)*zz(i,j,k)*x(1,jmax,kmax)
px2px3imax=(1-yy(i,j,k))*(1-zz(i,j,k))*x(imax,1,1)+yy(i,j,k)*(1-zz(i,j,k))*x(imax,jmax,1)+zz(i,j,k)* (1-yy(i,j,k))*x(imax,1,kmax)+yy(i,j,k)*zz(i,j,k)*x(imax,jmax,kmax)
px1px2px3=(1-xx(i,j,k))*px2px3i1+xx(i,j,k)*px2px3imax
py1=(1-xx(i,j,k))*y(1,j,k)+xx(i,j,k)*y(imax,j,k)
py2=(1-yy(i,j,k))*y(i,1,k)+yy(i,j,k)*y(i,jmax,k)
py3=(1-zz(i,j,k))*y(i,j,1)+zz(i,j,k)*y(i,j,kmax)
py1py2=(1-xx(i,j,k))*(1-yy(i,j,k))*y(1,1,k)+xx(i,j,k)*(1-yy(i,j,k))*y(imax,1,k)+yy(i,j,k)*(1-xx(i,j,k)) *y(1,jmax,k)+xx(i,j,k)*yy(i,j,k)*y(imax,jmax,k)
py2py3=(1-yy(i,j,k))*(1-zz(i,j,k))*y(i,1,1)+yy(i,j,k)*(1-zz(i,j,k))*y(i,jmax,1)+zz(i,j,k)*(1-yy(i,j,k))* y(i,1,kmax)+yy(i,j,k)*zz(i,j,k)*y(i,jmax,kmax)
py3py1=(1-zz(i,j,k))*(1-xx(i,j,k))*y(1,j,1)+zz(i,j,k)*(1-xx(i,j,k))*y(1,j,kmax)+xx(i,j,k)*(1-zz(i,j,k))* y(imax,j,1)+zz(i,j,k)*xx(i,j,k)*y(imax,j,kmax)
py2py3i1=(1-yy(i,j,k))*(1-zz(i,j,k))*y(1,1,1)+yy(i,j,k)*(1-zz(i,j,k))*y(1,jmax,1)+zz(i,j,k)*(1-yy(i,j,k) )*y(1,1,kmax)+yy(i,j,k)*zz(i,j,k)*y(1,jmax,kmax)
py2py3imax=(1-yy(i,j,k))*(1-zz(i,j,k))*y(imax,1,1)+yy(i,j,k)*(1-zz(i,j,k))*y(imax,jmax,1)+zz(i,j,k)* (1-yy(i,j,k))*y(imax,1,kmax)+yy(i,j,k)*zz(i,j,k)*y(imax,jmax,kmax)
py1py2py3=(1-xx(i,j,k))*py2py3i1+xx(i,j,k)*py2py3imax
pz1=(1-xx(i,j,k))*z(1,j,k)+xx(i,j,k)*z(imax,j,k)
pz2=(1-yy(i,j,k))*z(i,1,k)+yy(i,j,k)*z(i,jmax,k)
pz3=(1-zz(i,j,k))*z(i,j,1)+zz(i,j,k)*z(i,j,kmax)
pz1pz2=(1-xx(i,j,k))*(1-yy(i,j,k))*z(1,1,k)+xx(i,j,k)*(1-yy(i,j,k))*z(imax,1,k)+yy(i,j,k)*(1-xx(i,j,k)) *z(1,jmax,k)+xx(i,j,k)*yy(i,j,k)*z(imax,jmax,k)
pz2pz3=(1-yy(i,j,k))*(1-zz(i,j,k))*z(i,1,1)+yy(i,j,k)*(1-zz(i,j,k))*z(i,jmax,1)+zz(i,j,k)*(1-yy(i,j,k))*z (i,1,kmax)+yy(i,j,k)*zz(i,j,k)*z(i,jmax,kmax)
pz3pz1=(1-zz(i,j,k))*(1-xx(i,j,k))*z(1,j,1)+zz(i,j,k)*(1-xx(i,j,k))*z(1,j,kmax)+xx(i,j,k)*(1-zz(i,j,k))* z(imax,j,1)+zz(i,j,k)*xx(i,j,k)*z(imax,j,kmax)
pz2pz3i1=(1-yy(i,j,k))*(1-zz(i,j,k))*z(1,1,1)+yy(i,j,k)*(1-zz(i,j,k))*z(1,jmax,1)+zz(i,j,k)*(1-yy(i,j,k) )*z(1,1,kmax)+yy(i,j,k)*zz(i,j,k)*z(1,jmax,kmax)
pz2pz3imax=(1-yy(i,j,k))*(1-zz(i,j,k))*z(imax,1,1)+yy(i,j,k)*(1-zz(i,j,k))*z(imax,jmax,1)+zz(i,j,k)*( 1-yy(i,j,k))*z(imax,1,kmax)+yy(i,j,k)*zz(i,j,k)*z(imax,jmax,kmax)
pz1pz2pz3=(1-xx(i,j,k))*pz2pz3i1+xx(i,j,k)*pz2pz3imax
px=px1+px2+px3-px1px2-px2px3-px3px1+px1px2px3
py=py1+py2+py3-py1py2-py2py3-py3py1+py1py2py3
pz=pz1+pz2+pz3-pz1pz2-pz2pz3-pz3pz1+pz1pz2pz3
式中，
（1）px即(Pξ⊕Pη⊕Pζ)[x]；py即(Pξ⊕Pη⊕Pζ)[y]；pz即(Pξ⊕Pη⊕Pζ)[z]。

（2）px1即Pξ[x]；px2即Pη[x]；px3即Pζ[x]，其它类同。

（3）px1px2即PξPη[x]；px2px3即PηPζ[x] ；px3px1即PζPξ[x]；其它类同。

（4）px1px2px3即PξPηPζ[x]；py1py2py3即PξPηPζ[y]；pz1pz2pz3即PξPηPζ[z]。

（5）xx(i,j,k)为计算空间节点横坐标；x(i,j,k)为物理空间节点横坐标；其它类同。

3.2编程实例
程序需要读取物理空间的6个边界面的节点信息和计算空间的全部节点信息。

由于篇幅所限，不再贴上输入文件数据。

本例代码在Visual Fortran6.5编译环境下运行测试通过。

东岸线 2006-6-26 ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～
!/////////////////////////////////主程序//////////////////////////////////////////////////
program main
integer i,j,k,imax,jmax,kmax
!***打开物理平面边界坐标文件，先读取网格范围编号。

i从1～imax，j从1～jmax****** open(1,file='pbottom.dat',err=555)
read(1,*)
read(1,*) jmax,imax
close(1)
open(1,file='pleft.dat',err=555)
read(1,*)
read(1,*) kmax,jmax
close(1)
!************调用子程序******************************
call compute(imax,jmax,kmax)
!************返回主程序******************************
goto 888
555 print *,'ERROR OCCURED 读取文件出错！请检查无误后，重新进行计算。

'
888 print *,'Good Luck! ^_^'
print *,'----------东部动力网～倾情奉献～'
pause
end program
!//////////////////////////////////////////////////////////////////////////////////////////////
!/////////////////////////////////子程序//////////////////////////////////////////////////
subroutine compute(imax,jmax,kmax)
real
x(imax,jmax,kmax),y(imax,jmax,kmax),z(imax,jmax,kmax),xx(imax,jmax,kmax),yy(imax,jmax,km ax),zz(imax,jmax,kmax)
real
px1,px2,px3,px1px2,px2px3,px3px1,px1px2px3,px,py1,py2,py3,py1py2,py2py3,py3py1,py1py2py 3,py,pz1,pz2,pz3,pz1pz2,pz2pz3,pz3pz1,pz1pz2pz3,pz
!******打开物理平面文件，读入各边界面上节点坐标**************
!******从“水平的”bottom全边到“水平的”top全边*************
!******再到“垂直的”left中边到“垂直的”right中边*************
!******水平全边从左向右，垂直中边从下到上********************
open(1,file='pbottom.dat',err=555)
!************读入z=1面***********************
do 10 j=1,jmax
do 10 i=1,imax
read(1,*) x(i,j,1),y(i,j,1),z(i,j,1)
10 continue
!************读入z=kmax面********************
do 11 j=1,jmax
do 11 i=1,imax
read(1,*) x(i,j,kmax),y(i,j,kmax),z(i,j,kmax)
11 continue
!************读入y=1面***********************
do 12 k=1,kmax
do 12 i=1,imax
read(1,*) x(i,1,k),y(i,1,k),z(i,1,k)
12 continue
!************读入y=jmax面********************
do 13 k=1,kmax
do 13 i=1,imax
read(1,*) x(i,jmax,k),y(i,jmax,k),z(i,jmax,k)
13 continue
!************读入x=1面***********************
do 14 k=1,kmax
do 14 j=1,jmax
read(1,*) x(1,j,k),y(1,j,k),z(1,j,k)
14 continue
!************读入x=imax面********************
do 15 k=1,kmax
do 15 j=1,jmax
read(1,*) x(imax,j,k),y(imax,j,k),z(imax,j,k)
15 continue
close(1)
!*****从最下面逐行向上扫描读取计算区域的6个曲面各节点坐标值******* !*****从底往顶逐面、从下往上逐行、从左往右逐点扫描***************** open(1,file='computec.dat',err=555)
do 50 k=1,kmax
do 50 j=1,jmax
do 50 i=1,imax
read(1,*) xx(i,j,k),yy(i,j,k),zz(i,j,k)
50 continue
!***根据资料自行猜测、拓展到三维区域，网格内部节点坐标的计算公式***
!***注意理解该公式各项的具体含义和具体表达式***********************
!***如0、1在本例中则为1/1/1、imax/jmax/kmax************************
do 80 i=2,imax-1
do 80 j=2,jmax-1
do 80 k=2,kmax-1
!***********计算x坐标所需的各项具体表达式************************
px1=(1-xx(i,j,k))*x(1,j,k)+xx(i,j,k)*x(imax,j,k)
px2=(1-yy(i,j,k))*x(i,1,k)+yy(i,j,k)*x(i,jmax,k)
px3=(1-zz(i,j,k))*x(i,j,1)+zz(i,j,k)*x(i,j,kmax)
px1px2=(1-xx(i,j,k))*(1-yy(i,j,k))*x(1,1,k)+xx(i,j,k)*(1-yy(i,j,k))*x(imax,1,k)+yy(i,j,k)*(1-xx(i,j,k)) *x(1,jmax,k)+xx(i,j,k)*yy(i,j,k)*x(imax,jmax,k)
px2px3=(1-yy(i,j,k))*(1-zz(i,j,k))*x(i,1,1)+yy(i,j,k)*(1-zz(i,j,k))*x(i,jmax,1)+zz(i,j,k)*(1-yy(i,j,k))* x(i,1,kmax)+yy(i,j,k)*zz(i,j,k)*x(i,jmax,kmax)
px3px1=(1-zz(i,j,k))*(1-xx(i,j,k))*x(1,j,1)+zz(i,j,k)*(1-xx(i,j,k))*x(1,j,kmax)+xx(i,j,k)*(1-zz(i,j,k))* x(imax,j,1)+zz(i,j,k)*xx(i,j,k)*x(imax,j,kmax)
px2px3i1=(1-yy(i,j,k))*(1-zz(i,j,k))*x(1,1,1)+yy(i,j,k)*(1-zz(i,j,k))*x(1,jmax,1)+zz(i,j,k)*(1-yy(i,j,k) )*x(1,1,kmax)+yy(i,j,k)*zz(i,j,k)*x(1,jmax,kmax)
px2px3imax=(1-yy(i,j,k))*(1-zz(i,j,k))*x(imax,1,1)+yy(i,j,k)*(1-zz(i,j,k))*x(imax,jmax,1)+zz(i,j,k)* (1-yy(i,j,k))*x(imax,1,kmax)+yy(i,j,k)*zz(i,j,k)*x(imax,jmax,kmax)
px1px2px3=(1-xx(i,j,k))*px2px3i1+xx(i,j,k)*px2px3imax
!**************计算y坐标所需的各项具体表达式**********************
py1=(1-xx(i,j,k))*y(1,j,k)+xx(i,j,k)*y(imax,j,k)
py2=(1-yy(i,j,k))*y(i,1,k)+yy(i,j,k)*y(i,jmax,k)
py3=(1-zz(i,j,k))*y(i,j,1)+zz(i,j,k)*y(i,j,kmax)
py1py2=(1-xx(i,j,k))*(1-yy(i,j,k))*y(1,1,k)+xx(i,j,k)*(1-yy(i,j,k))*y(imax,1,k)+yy(i,j,k)*(1-xx(i,j,k)) *y(1,jmax,k)+xx(i,j,k)*yy(i,j,k)*y(imax,jmax,k)
py2py3=(1-yy(i,j,k))*(1-zz(i,j,k))*y(i,1,1)+yy(i,j,k)*(1-zz(i,j,k))*y(i,jmax,1)+zz(i,j,k)*(1-yy(i,j,k))* y(i,1,kmax)+yy(i,j,k)*zz(i,j,k)*y(i,jmax,kmax)
py3py1=(1-zz(i,j,k))*(1-xx(i,j,k))*y(1,j,1)+zz(i,j,k)*(1-xx(i,j,k))*y(1,j,kmax)+xx(i,j,k)*(1-zz(i,j,k))* y(imax,j,1)+zz(i,j,k)*xx(i,j,k)*y(imax,j,kmax)
py2py3i1=(1-yy(i,j,k))*(1-zz(i,j,k))*y(1,1,1)+yy(i,j,k)*(1-zz(i,j,k))*y(1,jmax,1)+zz(i,j,k)*(1-yy(i,j,k) )*y(1,1,kmax)+yy(i,j,k)*zz(i,j,k)*y(1,jmax,kmax)
py2py3imax=(1-yy(i,j,k))*(1-zz(i,j,k))*y(imax,1,1)+yy(i,j,k)*(1-zz(i,j,k))*y(imax,jmax,1)+zz(i,j,k)* (1-yy(i,j,k))*y(imax,1,kmax)+yy(i,j,k)*zz(i,j,k)*y(imax,jmax,kmax)
py1py2py3=(1-xx(i,j,k))*py2py3i1+xx(i,j,k)*py2py3imax
!**************计算z坐标所需的各项具体表达式***********************
pz1=(1-xx(i,j,k))*z(1,j,k)+xx(i,j,k)*z(imax,j,k)
pz2=(1-yy(i,j,k))*z(i,1,k)+yy(i,j,k)*z(i,jmax,k)
pz3=(1-zz(i,j,k))*z(i,j,1)+zz(i,j,k)*z(i,j,kmax)
pz1pz2=(1-xx(i,j,k))*(1-yy(i,j,k))*z(1,1,k)+xx(i,j,k)*(1-yy(i,j,k))*z(imax,1,k)+yy(i,j,k)*(1-xx(i,j,k)) *z(1,jmax,k)+xx(i,j,k)*yy(i,j,k)*z(imax,jmax,k)
pz2pz3=(1-yy(i,j,k))*(1-zz(i,j,k))*z(i,1,1)+yy(i,j,k)*(1-zz(i,j,k))*z(i,jmax,1)+zz(i,j,k)*(1-yy(i,j,k))*z (i,1,kmax)+yy(i,j,k)*zz(i,j,k)*z(i,jmax,kmax)
pz3pz1=(1-zz(i,j,k))*(1-xx(i,j,k))*z(1,j,1)+zz(i,j,k)*(1-xx(i,j,k))*z(1,j,kmax)+xx(i,j,k)*(1-zz(i,j,k))* z(imax,j,1)+zz(i,j,k)*xx(i,j,k)*z(imax,j,kmax)
pz2pz3i1=(1-yy(i,j,k))*(1-zz(i,j,k))*z(1,1,1)+yy(i,j,k)*(1-zz(i,j,k))*z(1,jmax,1)+zz(i,j,k)*(1-yy(i,j,k) )*z(1,1,kmax)+yy(i,j,k)*zz(i,j,k)*z(1,jmax,kmax)
pz2pz3imax=(1-yy(i,j,k))*(1-zz(i,j,k))*z(imax,1,1)+yy(i,j,k)*(1-zz(i,j,k))*z(imax,jmax,1)+zz(i,j,k)*( 1-yy(i,j,k))*z(imax,1,kmax)+yy(i,j,k)*zz(i,j,k)*z(imax,jmax,kmax)
pz1pz2pz3=(1-xx(i,j,k))*pz2pz3i1+xx(i,j,k)*pz2pz3imax
!**********获得最终的物理坐标系中各点坐标值**************************
px=px1+px2+px3-px1px2-px2px3-px3px1+px1px2px3
py=py1+py2+py3-py1py2-py2py3-py3py1+py1py2py3
pz=pz1+pz2+pz3-pz1pz2-pz2pz3-pz3pz1+pz1pz2pz3
x(i,j,k)=px
y(i,j,k)=py
z(i,j,k)=pz
80 continue
!**********输出Gridgen的PLOT3D格式的网格文件.grd*******************
open(1,file='3Dmesh.grd')
write(1,*) 1
write(1,*) imax,jmax,kmax
do 100 k=1,kmax
do 100 j=1,jmax
do 100 i=1,imax
write(1,*) x(i,j,k)
100 continue
do 200 k=1,kmax
do 200 j=1,jmax
do 200 i=1,imax
write(1,*) y(i,j,k)
200 continue
do 300 k=1,kmax
do 300 j=1,jmax
do 300 i=1,imax
write(1,*) z(i,j,k)
300 continue
close(1)
goto 888
555 print *,'ERROR OCCURED 读取文件出错！请检查无误后，重新进行计算。

'
pause
888 return
end subroutine ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～本例生成的TFI体网格效果参见封面图片。

4. 参考文献
（1）数值传热学（第二版）。

（2）Mesh Generation。

（来自网络）
（3）Numerical Grid Generation Foundations and Applications。

（来自网络）
5. 版权声明
未经作者许可严禁将本文内容以任何形式用于商业用途!
官方论坛:/cnbbs
E-MAIL：ecengine@。