基于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.