基于VB与Fortran混合编程的混凝土坝温度应力仿真分析软件开发

郭祥伟,陈国荣,刘银芳

(河海大学力学与材料学院,江苏南京 210098)

随着计算机技术的快速发展和混凝土坝在我国的建设进程,有限元软件在混凝土坝温度应力仿真分析中的应用越来越普遍。近年来,以Fortran为核心的有限元仿真程序为解决许多工程技术问题提供了极大的便利[1]。但有限元前后处理的可视化存在较大困难,因此,应用软件的可视化日益得到工程界的重视。在Windows操作平台下,VB是用于开发和创建具有图形用户界面应用程序的强有力工具之一[2],而且可以方便地调用OpenGL图形库,用户可以使用OpenGL快速开发出高质量的三维图形。本文基于Fortran源程序结合VB仿真平台,采用混合编程的方法实现混凝土坝温度应力仿真分析的可视化。

1 有限元计算前后期处理

混凝土坝温度应力有限元计算分析[3-4]以Fortran源程序为核心,前后处理采用VB软件设计开发,全面实现混凝土坝有限元网格的自动剖分,有限元网格结点与单元信息、边界条件与荷载信息的获取,后期处理包括生成各种计算成果的曲线分布图等,所有模块均为子程序定义的程序单元。

1.1 空间有限元网格剖分算法

当空间结构的边界面比较复杂或结构内部区域的材料不均匀,或者要考虑网格的布置要求时,一般都将空间结构划分为若干空间子区域(超单元),并给出各超单元的单元信息、结点坐标以及各个方向的剖分结点数、剖分比例因子和材料类型号。超单元的定义一般要根据结构的几何特征以及内部区域材料类型、网格的布置要求来确定,并将其变换到相应的规则区域中,求出规则区域中网格结点的局部坐标,然后利用参数变换求出单元内各网格结点的整体坐标,并相应地定义结点和单元编号以及单元信息,对于边界上的公共点还必须进行重复结点的消除。

1.1.1空间网格结点的局部坐标和整体坐标

由于混凝土坝结构边界的复杂性和内部结构材料的差异性,将其划分为若干个超单元。下面对八结点六面体超单元做详细描述,当结构的外形轮廓线呈曲线形状时,可在边界上增加结点,从而采取高次空间超单元进行剖分。空间子区域的几何形状由8个基本结点p1～p8确定,其结点坐标为(xα,yα,zα),α=1,2,…,8,单元信息定义为 M=(p1,p2,…,p8)。将该超单元变换到规则的正方体(标准单元或母单元)区域中,此时,超单元的各条棱边为对 应于 正 方体 的,…,假定 l,m,n为对应于正方体方向生成的结点数,现用等值面将正方体在方向上划分为l-1,m-1和 n-1个单元,从而得到一个规则的空间立体网格。正方体网格中任意一点的局部坐标为[5]

式中 :λξ,λη,λζ为剖分比例因子。

而超单元中对应网格结点的整体坐标如下:

其中

式中:φα为局部坐标插值函数;(xα,yα,zα)为超单元的8个已知结点坐标。

1.1.2空间结点、单元编号及单元信息

在空间超单元变换剖分法中,网格的结点编号、单元编号和单元信息的定义可在规则区域网格中进行,对于不同的空间单元类型,其结点、单元的自动编号及单元信息的定义有所不同。当生成的单元类型为空间八结点六面体单元时,假定结点的编号累计顺序是先从 ζ=1到 ζ=-1,再从 η=1到 η=-1,最后从ξ=1到 ξ=-1,那么网格内任一点的一维编号为[5]如果假定单元的自动编号顺序和结点编号顺序一致,则网格内任一个单元的编号为而对应的单元信息可定义为

1.2 等值线生成算法

等值线法是对大量离散的具有一定规律的物理量值进行处理,并用图形表达这些量值的分布规律的方法,其特性主要表现为:等值线上的点对应的物理量值相等;等值线是连续曲线,可能是封闭曲线,也可能由于目标区域的限定,从区域的边界开始到区域的另一边界结束;也可能某一标称值对应有多条等值线。

本文结合混凝土坝温度应力有限元模拟的特点,提出一种简单易行的等值线生成算法,该算法的思想是将整个混凝土坝的等值线划分转化为对单元进行等值线划分;体单元可以通过剖面或体单元表面转化为平面单元,对每个单元进行判断,根据判断结果决定是否需要细分[6];然后在单元内部采用直接网格法进行等值点的计算,再按一定的方位判别法连接各等值点得到等值线。由于单元内部等值点成对出现,对于四边形单元和三角形单元,遍历所有单元,当单元内等值点数目不超过2个时,直接将各单元的等值点连接起来。这样生成的等值线要么是闭合的,要么终止于边界上,且不相互交错。该算法的实现步骤如下:计算单元的物理量数值范围,求出最大和最小物理量;计算最大、最小物理量对应的等值线范围区间;对单元进行细分,对每一次细分得到的小单元重新进行上述步骤,对所有单元进行循环计算。

按上述方法生成等值线,直接计算生成一系列直线段,绘出所有线段便组成了各条等值线,但是这样生成的等值线表现为折线图,当网格比较稀疏而且剖切面不垂直于坐标轴时,不光滑程度较高,甚至表现为锯齿状,因此应使用样条曲线拟合方法[7]使其光滑。

2 OpenGL在VB环境中的应用

2.1 Vbog1.tlb的使用

在VB环境中应用OpenGL进行三维模型的设计操作大多通过第三方函数库VB OpenGL type library(Vbog1.tlb)来进行,可省去大量的底层编程工作,在应用程序设计中可以起到事半功倍的效果,TLB是一种OLE或ActiveX定义文件,它包括常数、类、接口等的定义。将Vbog1.tlb应用于Microsoft平台的具体操作方法为:①将Vbog1.tlb安装在适当的工作目录;②注册Vbog1.tlb,可以使用regsvr32.exe进行注册,也可以在VB集成环境的“工程”菜单下通过“引用”子菜单将TLB文件加入项目;③在“对象浏览器”中查看该函数库;④在代码窗口中调用OpenGL函数。

2.2 VB环境下OpenGL程序框架的建立

基于VB程序框架,在VB集成环境的“工程”菜单下通过“引用”子菜单将VB OpenGL API加入到项目中,在代码窗口中调用OpenGL函数,运用网格剖分、等值线和云图等算法,可以在VB环境中应用OpenGL进行三维模型的设计操作[8-9]。

a.编写窗体的Load事件过程。此过程主要设置像素格式,建立渲染描述表,设置投影矩阵的模式,利用显示列表建立三维图形。具体代码如下:

Private Sub Form_load()

Dim hGLRC As long

Dim pfd As PIXELFOR MATDESCRIPTOR '定义像素格式

pfd.nSize=Len(pfd) '设置像素格式大小

pfd.nVersion=1 '版本号

PixelFormat:ChoosePixelFormat hDC.pfd '匹配设备上下文的像素与指定的像素

SetPixelFormat hDC,PixelFormat,pfd '设置设备上下文的像素格式

hGLRC=wglCreateContext hDC '创建渲染描述表

wglMakeCurrent hDC,hGLRC '设置为当前渲染描述表

glViewport…… '设置视区

Form.Paint '窗体重绘

End Sub

b.编写窗体的Paint事件过程。当窗体重绘时发生Paint事件,在此过程中设置视点,调用显示列表并显示所绘图形。具体代码如下:

Private Sub Form_Paint()

glLoadIdentity '重置当前指定的矩阵为单位矩阵,即初始化矩阵

gluLookAt…… '设置视点

glCallList…… '调用显示列表

SwapBuffers hDC '设置双缓冲

End Sub

c.编写窗体的Unload事件过程。在此过程中主要删除hGLRC和hPalete。具体代码如下:

Private Sub Form.Unload(Cancel As Integer)

If hGLRC＜＞0 Then

wglMakeCurrent 0,0

wglDeleteContext hGLRC

End If

If hPalette＜＞0 Then

DeleteObject hPalette

End If

End Sub

3 VB和Fortran混合编程的实现方法

实现VB和Fortran2种语言的混合编程,首先要解决相关接口问题,包括2种语言之间的调用约定,数据的传递和2种语言的命名约定等。在正确处理和充分考虑相关接口后,采用call()函数实现VB调用Fortran生成的动态链接库DLL文件的目的。动态链接库DLL文件是一个库函数,可以是一个函数也可以是多个函数,且这些函数独立与主程序进行编译、链接和储存,可为程序提供良好的设计平台。Fortran动态链接库建立的具体步骤参见文献[10-11]。结合混凝土坝温度应力数值仿真Fortran源程序,VB和Fortran混合编程的具体步骤包括:

步骤1创建通用的Fortran计算源程序[12],代码如下:

SUBROUTINEWENDU(da,re,num)

1 ！DEC＄ATTRIBUTES DLLEXPOR T::WENDU

INTEGER*4::da,re,i……

REAL*8::num(12)……

2 OPEN(20,file='data_direct1.txt',status='old',access='direct',form='formatted',recl=158)

3 CALL STRESS(u,v,w,sx,sy,sz,sxy,syz,szx,s1,s2,s3)

4 WRITE(20,100,rec=i)u(i),v(i),w(i),sx(i),sy(i),sz(i),sxy(i),sxz(i),syz(i),s1(i),s2(i),s3(i)

RETUR N

END SUBROUTINE

语句1用于建立VB和Fortran接口的引入点,ATTRIBUTES是Fortran的元命令,用于声明微软扩展属性,DLLEXPORT属性表明该函数或子过程能被其他程序或DLL调用,WENDU为该程序的过程名。语句2的作用是打开(建立)一个名为'data_direct1.txt'的文件,文件号为20,文件中每个记录长度为158B,采用直接(随机)有格式的存取方式。语句3表示调用温度应力计算模块。语句4表示在20号文件中按照格式输出每一个记录的数据。RETURN语句把Fortran计算结果返回VB主程序。

步骤2建立VB标准工程[12-13],代码如下:

5 Private Declare Sub WENDU lib″WENDU.dll″Alias″_WENDU@12″(da as Long,re as Long,num as Double)

Private Sub Readdata_Click()

6 CommonDialog1.Filter=″text|*.txt|All Files|*.*″CommonDialog1.ShowOpen fname=CommonDialog1.FileName

Dim day as Long,re as Long,filenum as Long……

Dim num(1 to 12)as Double

7 Call WENDU(da,re,num(1))

8 Open fname For Random As filenum Len=158

9 Get#filenum,sumnum,numrec

End Sub

语句5为VB调用Fortran动态链接库的声明语句,Alias属性表示子过程的标识符,以防子过程名中出现VB不可识别的字符,括号内为定义的参数。语句6表示将从打开文件对话框中选定的结果数据文件及其路径赋给fname变量。语句7体现数据的传递,可以看出只要将数组的第1个元素传递到DLL过程中,就可以访问数组的所有元素。语句8表示打开结果数据文件。语句9为VB读取结果数据文件。

4 观音岩水电站混凝土坝温度应力仿真分析

观音岩水电站位于云南省丽江市华坪县与四川省攀枝花市交界的金沙江中游河段,电站采用混合坝坝型(碾压混凝土重力坝结合右岸心墙堆石坝),坝顶总长1158m,其中混凝土部分长838m,坝顶高程为1139m,最大坝高为159m。坝区水温和气象条件参照坝区多年气温统计表。在环境温度一定的条件下,坝体混凝土的温度主要由混凝土的水化热、浇筑温度和冷却水管水温共同控制,温度应力主要在施工期和运行期出现。在施工过程中,混凝土采用逐层铺设,坝体不同高程的混凝土浇筑时间可能相差1年甚至几年,因此,坝体温度场计算必须考虑时间变化的因素。

将基于VB和Fortran混合编程的混凝土坝温度应力有限元仿真分析软件应用于观音岩水电站混凝土重力坝中,并进行合理的简化,对左岸1～18号坝段及坝基建立有限元模型进行仿真分析。图1给出了坝段有限元网格与坝体不同混凝土分区,图2给出了第301天坝体沿坐标轴剖面(x=5m)的温度等值线,图3给出了运行期坝体沿坐标轴剖面(x=5 m)的 σy分布,图 4给出了(5 m,-0.5m,1080m)处第一主应力 σ1的时间历程,由 σ1的变化曲线可以看出其值始终满足混凝土强度要求。软件在实际工程中的应用表明,基于VB与Fortran混合编程的混凝土坝温度应力仿真分析软件能提高应用程序的工作效率,提供友好的人机交互界面。

图1 坝段有限元网格与坝体不同混凝土分区

图2 第301天坝体沿坐标轴剖面(x=5m)的温度等值线(单位:℃)

图3 运行期坝体沿坐标轴剖面(x=5m)的σy分布(单位:kPa)

图4 (5m,-0.5m,1080m)处第一主应力σ1的时间历程

5 结 语

采用VB与Fortran混合编程的方法,以动态链接库DLL文件作接口,开发出混凝土坝温度应力仿真分析软件。该软件充分利用了Fortran语言强大的科学计算能力和VB开发界面的强大功能,可实现程序资源的共享和可视化,提高程序的工作效率,

同时具有友好的前后处理界面,能够圆满地实现软件的预期功能。程序中原始数据的建立由VB交互式界面承担,通过调用Fortran动态链接库实现计算功能,再将计算数据传回VB程序并利用VB可视化功能显示和分析。该软件在实际工程的可视化仿真应用表明,基于VB与Fortran混合编程的混凝土坝温度应力仿真分析软件具有2种计算机语言的优点,能提高应用程序的工作效率,提供友好的人机交互界面,在前后处理方面具有明显的优势。

[1]陈国荣.有限单元法原理及应用[M].北京:科学出版社,2009.

[2]牛又奇,孙建国.新编Visual Basic程序设计教程[M].苏州:苏州大学出版社,2002.

[3]朱伯芳.考虑水管冷却效果的混凝土等效热传导方程[J].水利学报,1991,22(3):28-34.

[4]朱伯芳.大体积混凝土温度应力与温度控制[M].北京:中国电力出版社,1999.

[5]陈和群,彭宣茂.有限元法微机程序与图形处理[M].南京:河海大学出版社,1992.

[6]向毅斌,吴诗忄享,刘琦.材料成形模拟中等值线图的生成[J].磨具技术,2001,1(1):4-6.

[7]吴六永.样条曲线(NURBS)原理及其实现技术研究[J].西南交通大学学报,1999,34(6):683-688.

[8]白燕斌,史惠康.OpenGL三维图形库编程指南[M].北京:机械工业出版社,1998.

[9]WRIGHT R S Jr.OpenGL超级宝典[M].4版.北京:人民邮电出版社,2010.

[10]马进荣,王永勇,谢敏.VB与FORTRAN混合编程在河口潮流计算中的应用[J].人民珠江,2005(3):81-82.

[11]徐林春,赵明登,童汉毅.VB与FORTRAN混合编程及其在流动数值模拟可视化技术中的应用[J].武汉大学学报:工学版,2004,37(2):21-24.

[12]彭国伦.Fortran 95程序设计[M].北京:中国电力出版社,2002.

[13]欧阳永忠,王瑞,陆秀平.VC、VB与FORTRAN的混合编程技术及其实现[J].海洋测绘,2004,24(1):54-58.