郭祥偉,陳國榮,劉銀芳

(河海大學(xué)力學(xué)與材料學(xué)院,江蘇南京 210098)

隨著計算機技術(shù)的快速發(fā)展和混凝土壩在我國的建設(shè)進程,有限元軟件在混凝土壩溫度應(yīng)力仿真分析中的應(yīng)用越來越普遍。近年來,以Fortran為核心的有限元仿真程序為解決許多工程技術(shù)問題提供了極大的便利[1]。但有限元前后處理的可視化存在較大困難,因此,應(yīng)用軟件的可視化日益得到工程界的重視。在Windows操作平臺下,VB是用于開發(fā)和創(chuàng)建具有圖形用戶界面應(yīng)用程序的強有力工具之一[2],而且可以方便地調(diào)用OpenGL圖形庫,用戶可以使用OpenGL快速開發(fā)出高質(zhì)量的三維圖形。本文基于Fortran源程序結(jié)合VB仿真平臺,采用混合編程的方法實現(xiàn)混凝土壩溫度應(yīng)力仿真分析的可視化。

1 有限元計算前后期處理

混凝土壩溫度應(yīng)力有限元計算分析[3-4]以Fortran源程序為核心,前后處理采用VB軟件設(shè)計開發(fā),全面實現(xiàn)混凝土壩有限元網(wǎng)格的自動剖分,有限元網(wǎng)格結(jié)點與單元信息、邊界條件與荷載信息的獲取,后期處理包括生成各種計算成果的曲線分布圖等,所有模塊均為子程序定義的程序單元。

1.1 空間有限元網(wǎng)格剖分算法

當(dāng)空間結(jié)構(gòu)的邊界面比較復(fù)雜或結(jié)構(gòu)內(nèi)部區(qū)域的材料不均勻,或者要考慮網(wǎng)格的布置要求時,一般都將空間結(jié)構(gòu)劃分為若干空間子區(qū)域(超單元),并給出各超單元的單元信息、結(jié)點坐標(biāo)以及各個方向的剖分結(jié)點數(shù)、剖分比例因子和材料類型號。超單元的定義一般要根據(jù)結(jié)構(gòu)的幾何特征以及內(nèi)部區(qū)域材料類型、網(wǎng)格的布置要求來確定,并將其變換到相應(yīng)的規(guī)則區(qū)域中,求出規(guī)則區(qū)域中網(wǎng)格結(jié)點的局部坐標(biāo),然后利用參數(shù)變換求出單元內(nèi)各網(wǎng)格結(jié)點的整體坐標(biāo),并相應(yīng)地定義結(jié)點和單元編號以及單元信息,對于邊界上的公共點還必須進行重復(fù)結(jié)點的消除。

1.1.1空間網(wǎng)格結(jié)點的局部坐標(biāo)和整體坐標(biāo)

由于混凝土壩結(jié)構(gòu)邊界的復(fù)雜性和內(nèi)部結(jié)構(gòu)材料的差異性,將其劃分為若干個超單元。下面對八結(jié)點六面體超單元做詳細描述,當(dāng)結(jié)構(gòu)的外形輪廓線呈曲線形狀時,可在邊界上增加結(jié)點,從而采取高次空間超單元進行剖分?？臻g子區(qū)域的幾何形狀由8個基本結(jié)點p1～p8確定,其結(jié)點坐標(biāo)為(xα,yα,zα),α=1,2,…,8,單元信息定義為 M=(p1,p2,…,p8)。將該超單元變換到規(guī)則的正方體(標(biāo)準(zhǔn)單元或母單元)區(qū)域中,此時,超單元的各條棱邊為對 應(yīng)于 正 方體 的,…,假定 l,m,n為對應(yīng)于正方體方向生成的結(jié)點數(shù),現(xiàn)用等值面將正方體在方向上劃分為l-1,m-1和 n-1個單元,從而得到一個規(guī)則的空間立體網(wǎng)格。正方體網(wǎng)格中任意一點的局部坐標(biāo)為[5]

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

而超單元中對應(yīng)網(wǎng)格結(jié)點的整體坐標(biāo)如下:

其中

式中:φα為局部坐標(biāo)插值函數(shù);(xα,yα,zα)為超單元的8個已知結(jié)點坐標(biāo)。

1.1.2空間結(jié)點、單元編號及單元信息

在空間超單元變換剖分法中,網(wǎng)格的結(jié)點編號、單元編號和單元信息的定義可在規(guī)則區(qū)域網(wǎng)格中進行,對于不同的空間單元類型,其結(jié)點、單元的自動編號及單元信息的定義有所不同。當(dāng)生成的單元類型為空間八結(jié)點六面體單元時,假定結(jié)點的編號累計順序是先從 ζ=1到 ζ=-1,再從 η=1到 η=-1,最后從ξ=1到 ξ=-1,那么網(wǎng)格內(nèi)任一點的一維編號為[5]如果假定單元的自動編號順序和結(jié)點編號順序一致,則網(wǎng)格內(nèi)任一個單元的編號為而對應(yīng)的單元信息可定義為

1.2 等值線生成算法

等值線法是對大量離散的具有一定規(guī)律的物理量值進行處理,并用圖形表達這些量值的分布規(guī)律的方法,其特性主要表現(xiàn)為:等值線上的點對應(yīng)的物理量值相等;等值線是連續(xù)曲線,可能是封閉曲線,也可能由于目標(biāo)區(qū)域的限定,從區(qū)域的邊界開始到區(qū)域的另一邊界結(jié)束;也可能某一標(biāo)稱值對應(yīng)有多條等值線。

本文結(jié)合混凝土壩溫度應(yīng)力有限元模擬的特點,提出一種簡單易行的等值線生成算法,該算法的思想是將整個混凝土壩的等值線劃分轉(zhuǎn)化為對單元進行等值線劃分;體單元可以通過剖面或體單元表面轉(zhuǎn)化為平面單元,對每個單元進行判斷,根據(jù)判斷結(jié)果決定是否需要細分[6];然后在單元內(nèi)部采用直接網(wǎng)格法進行等值點的計算,再按一定的方位判別法連接各等值點得到等值線。由于單元內(nèi)部等值點成對出現(xiàn),對于四邊形單元和三角形單元,遍歷所有單元,當(dāng)單元內(nèi)等值點數(shù)目不超過2個時,直接將各單元的等值點連接起來。這樣生成的等值線要么是閉合的,要么終止于邊界上,且不相互交錯。該算法的實現(xiàn)步驟如下:計算單元的物理量數(shù)值范圍,求出最大和最小物理量;計算最大、最小物理量對應(yīng)的等值線范圍區(qū)間;對單元進行細分,對每一次細分得到的小單元重新進行上述步驟,對所有單元進行循環(huán)計算。

按上述方法生成等值線,直接計算生成一系列直線段,繪出所有線段便組成了各條等值線,但是這樣生成的等值線表現(xiàn)為折線圖,當(dāng)網(wǎng)格比較稀疏而且剖切面不垂直于坐標(biāo)軸時,不光滑程度較高,甚至表現(xiàn)為鋸齒狀,因此應(yīng)使用樣條曲線擬合方法[7]使其光滑。

2 OpenGL在VB環(huán)境中的應(yīng)用

2.1 Vbog1.tlb的使用

在VB環(huán)境中應(yīng)用OpenGL進行三維模型的設(shè)計操作大多通過第三方函數(shù)庫VB OpenGL type library(Vbog1.tlb)來進行,可省去大量的底層編程工作,在應(yīng)用程序設(shè)計中可以起到事半功倍的效果,TLB是一種OLE或ActiveX定義文件,它包括常數(shù)、類、接口等的定義。將Vbog1.tlb應(yīng)用于Microsoft平臺的具體操作方法為:①將Vbog1.tlb安裝在適當(dāng)?shù)墓ぷ髂夸?②注冊Vbog1.tlb,可以使用regsvr32.exe進行注冊,也可以在VB集成環(huán)境的“工程”菜單下通過“引用”子菜單將TLB文件加入項目;③在“對象瀏覽器”中查看該函數(shù)庫;④在代碼窗口中調(diào)用OpenGL函數(shù)。

2.2 VB環(huán)境下OpenGL程序框架的建立

基于VB程序框架,在VB集成環(huán)境的“工程”菜單下通過“引用”子菜單將VB OpenGL API加入到項目中,在代碼窗口中調(diào)用OpenGL函數(shù),運用網(wǎng)格剖分、等值線和云圖等算法,可以在VB環(huán)境中應(yīng)用OpenGL進行三維模型的設(shè)計操作[8-9]。

a.編寫窗體的Load事件過程。此過程主要設(shè)置像素格式,建立渲染描述表,設(shè)置投影矩陣的模式,利用顯示列表建立三維圖形。具體代碼如下:

Private Sub Form_load()

Dim hGLRC As long

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

pfd.nSize=Len(pfd) '設(shè)置像素格式大小

pfd.nVersion=1 '版本號

PixelFormat:ChoosePixelFormat hDC.pfd '匹配設(shè)備上下文的像素與指定的像素

SetPixelFormat hDC,PixelFormat,pfd '設(shè)置設(shè)備上下文的像素格式

hGLRC=wglCreateContext hDC '創(chuàng)建渲染描述表

wglMakeCurrent hDC,hGLRC '設(shè)置為當(dāng)前渲染描述表

glViewport…… '設(shè)置視區(qū)

Form.Paint '窗體重繪

End Sub

b.編寫窗體的Paint事件過程。當(dāng)窗體重繪時發(fā)生Paint事件,在此過程中設(shè)置視點,調(diào)用顯示列表并顯示所繪圖形。具體代碼如下:

Private Sub Form_Paint()

glLoadIdentity '重置當(dāng)前指定的矩陣為單位矩陣,即初始化矩陣

gluLookAt…… '設(shè)置視點

glCallList…… '調(diào)用顯示列表

SwapBuffers hDC '設(shè)置雙緩沖

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混合編程的實現(xiàn)方法

實現(xiàn)VB和Fortran2種語言的混合編程,首先要解決相關(guān)接口問題,包括2種語言之間的調(diào)用約定,數(shù)據(jù)的傳遞和2種語言的命名約定等。在正確處理和充分考慮相關(guān)接口后,采用call()函數(shù)實現(xiàn)VB調(diào)用Fortran生成的動態(tài)鏈接庫DLL文件的目的。動態(tài)鏈接庫DLL文件是一個庫函數(shù),可以是一個函數(shù)也可以是多個函數(shù),且這些函數(shù)獨立與主程序進行編譯、鏈接和儲存,可為程序提供良好的設(shè)計平臺。Fortran動態(tài)鏈接庫建立的具體步驟參見文獻[10-11]。結(jié)合混凝土壩溫度應(yīng)力數(shù)值仿真Fortran源程序,VB和Fortran混合編程的具體步驟包括:

步驟1創(chuàng)建通用的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屬性表明該函數(shù)或子過程能被其他程序或DLL調(diào)用,WENDU為該程序的過程名。語句2的作用是打開(建立)一個名為'data_direct1.txt'的文件,文件號為20,文件中每個記錄長度為158B,采用直接(隨機)有格式的存取方式。語句3表示調(diào)用溫度應(yīng)力計算模塊。語句4表示在20號文件中按照格式輸出每一個記錄的數(shù)據(jù)。RETURN語句把Fortran計算結(jié)果返回VB主程序。

步驟2建立VB標(biāo)準(zhǔn)工程[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調(diào)用Fortran動態(tài)鏈接庫的聲明語句,Alias屬性表示子過程的標(biāo)識符,以防子過程名中出現(xiàn)VB不可識別的字符,括號內(nèi)為定義的參數(shù)。語句6表示將從打開文件對話框中選定的結(jié)果數(shù)據(jù)文件及其路徑賦給fname變量。語句7體現(xiàn)數(shù)據(jù)的傳遞,可以看出只要將數(shù)組的第1個元素傳遞到DLL過程中,就可以訪問數(shù)組的所有元素。語句8表示打開結(jié)果數(shù)據(jù)文件。語句9為VB讀取結(jié)果數(shù)據(jù)文件。

4 觀音巖水電站混凝土壩溫度應(yīng)力仿真分析

觀音巖水電站位于云南省麗江市華坪縣與四川省攀枝花市交界的金沙江中游河段,電站采用混合壩壩型(碾壓混凝土重力壩結(jié)合右岸心墻堆石壩),壩頂總長1158m,其中混凝土部分長838m,壩頂高程為1139m,最大壩高為159m。壩區(qū)水溫和氣象條件參照壩區(qū)多年氣溫統(tǒng)計表。在環(huán)境溫度一定的條件下,壩體混凝土的溫度主要由混凝土的水化熱、澆筑溫度和冷卻水管水溫共同控制,溫度應(yīng)力主要在施工期和運行期出現(xiàn)。在施工過程中,混凝土采用逐層鋪設(shè),壩體不同高程的混凝土澆筑時間可能相差1年甚至幾年,因此,壩體溫度場計算必須考慮時間變化的因素。

將基于VB和Fortran混合編程的混凝土壩溫度應(yīng)力有限元仿真分析軟件應(yīng)用于觀音巖水電站混凝土重力壩中,并進行合理的簡化,對左岸1～18號壩段及壩基建立有限元模型進行仿真分析。圖1給出了壩段有限元網(wǎng)格與壩體不同混凝土分區(qū),圖2給出了第301天壩體沿坐標(biāo)軸剖面(x=5m)的溫度等值線,圖3給出了運行期壩體沿坐標(biāo)軸剖面(x=5 m)的 σy分布,圖 4給出了(5 m,-0.5m,1080m)處第一主應(yīng)力 σ1的時間歷程,由 σ1的變化曲線可以看出其值始終滿足混凝土強度要求。軟件在實際工程中的應(yīng)用表明,基于VB與Fortran混合編程的混凝土壩溫度應(yīng)力仿真分析軟件能提高應(yīng)用程序的工作效率,提供友好的人機交互界面。

圖1 壩段有限元網(wǎng)格與壩體不同混凝土分區(qū)

圖2 第301天壩體沿坐標(biāo)軸剖面(x=5m)的溫度等值線(單位:℃)

圖3 運行期壩體沿坐標(biāo)軸剖面(x=5m)的σy分布(單位:kPa)

圖4 (5m,-0.5m,1080m)處第一主應(yīng)力σ1的時間歷程

5 結(jié) 語

采用VB與Fortran混合編程的方法,以動態(tài)鏈接庫DLL文件作接口,開發(fā)出混凝土壩溫度應(yīng)力仿真分析軟件。該軟件充分利用了Fortran語言強大的科學(xué)計算能力和VB開發(fā)界面的強大功能,可實現(xiàn)程序資源的共享和可視化,提高程序的工作效率,

同時具有友好的前后處理界面,能夠圓滿地實現(xiàn)軟件的預(yù)期功能。程序中原始數(shù)據(jù)的建立由VB交互式界面承擔(dān),通過調(diào)用Fortran動態(tài)鏈接庫實現(xiàn)計算功能,再將計算數(shù)據(jù)傳回VB程序并利用VB可視化功能顯示和分析。該軟件在實際工程的可視化仿真應(yīng)用表明,基于VB與Fortran混合編程的混凝土壩溫度應(yīng)力仿真分析軟件具有2種計算機語言的優(yōu)點,能提高應(yīng)用程序的工作效率,提供友好的人機交互界面,在前后處理方面具有明顯的優(yōu)勢。

[1]陳國榮.有限單元法原理及應(yīng)用[M].北京:科學(xué)出版社,2009.

[2]牛又奇,孫建國.新編Visual Basic程序設(shè)計教程[M].蘇州:蘇州大學(xué)出版社,2002.

[3]朱伯芳.考慮水管冷卻效果的混凝土等效熱傳導(dǎo)方程[J].水利學(xué)報,1991,22(3):28-34.

[4]朱伯芳.大體積混凝土溫度應(yīng)力與溫度控制[M].北京:中國電力出版社,1999.

[5]陳和群,彭宣茂.有限元法微機程序與圖形處理[M].南京:河海大學(xué)出版社,1992.

[6]向毅斌,吳詩忄享,劉琦.材料成形模擬中等值線圖的生成[J].磨具技術(shù),2001,1(1):4-6.

[7]吳六永.樣條曲線(NURBS)原理及其實現(xiàn)技術(shù)研究[J].西南交通大學(xué)學(xué)報,1999,34(6):683-688.

[8]白燕斌,史惠康.OpenGL三維圖形庫編程指南[M].北京:機械工業(yè)出版社,1998.

[9]WRIGHT R S Jr.OpenGL超級寶典[M].4版.北京:人民郵電出版社,2010.

[10]馬進榮,王永勇,謝敏.VB與FORTRAN混合編程在河口潮流計算中的應(yīng)用[J].人民珠江,2005(3):81-82.

[11]徐林春,趙明登,童漢毅.VB與FORTRAN混合編程及其在流動數(shù)值模擬可視化技術(shù)中的應(yīng)用[J].武漢大學(xué)學(xué)報:工學(xué)版,2004,37(2):21-24.

[12]彭國倫.Fortran 95程序設(shè)計[M].北京:中國電力出版社,2002.

[13]歐陽永忠,王瑞,陸秀平.VC、VB與FORTRAN的混合編程技術(shù)及其實現(xiàn)[J].海洋測繪,2004,24(1):54-58.