郭祥偉,陳國(guó)榮,劉銀芳

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

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

1 有限元計(jì)算前后期處理

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

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

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

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

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

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

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

其中

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

1.1.2空間結(jié)點(diǎn)、單元編號(hào)及單元信息

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

1.2 等值線生成算法

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

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

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

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

2.1 Vbog1.tlb的使用

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

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

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

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)窗體重繪時(shí)發(fā)生Paint事件,在此過程中設(shè)置視點(diǎn),調(diào)用顯示列表并顯示所繪圖形。具體代碼如下:

Private Sub Form_Paint()

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

gluLookAt…… '設(shè)置視點(diǎn)

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

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

End Sub

c.編寫窗體的Unload事件過程。在此過程中主要?jiǎng)h除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混合編程的實(shí)現(xiàn)方法

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

步驟1創(chuàng)建通用的Fortran計(jì)算源程序[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接口的引入點(diǎn),ATTRIBUTES是Fortran的元命令,用于聲明微軟擴(kuò)展屬性,DLLEXPORT屬性表明該函數(shù)或子過程能被其他程序或DLL調(diào)用,WENDU為該程序的過程名。語句2的作用是打開(建立)一個(gè)名為'data_direct1.txt'的文件,文件號(hào)為20,文件中每個(gè)記錄長(zhǎng)度為158B,采用直接(隨機(jī))有格式的存取方式。語句3表示調(diào)用溫度應(yīng)力計(jì)算模塊。語句4表示在20號(hào)文件中按照格式輸出每一個(gè)記錄的數(shù)據(jù)。RETURN語句把Fortran計(jì)算結(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動(dòng)態(tài)鏈接庫(kù)的聲明語句,Alias屬性表示子過程的標(biāo)識(shí)符,以防子過程名中出現(xiàn)VB不可識(shí)別的字符,括號(hào)內(nèi)為定義的參數(shù)。語句6表示將從打開文件對(duì)話框中選定的結(jié)果數(shù)據(jù)文件及其路徑賦給fname變量。語句7體現(xiàn)數(shù)據(jù)的傳遞,可以看出只要將數(shù)組的第1個(gè)元素傳遞到DLL過程中,就可以訪問數(shù)組的所有元素。語句8表示打開結(jié)果數(shù)據(jù)文件。語句9為VB讀取結(jié)果數(shù)據(jù)文件。

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

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

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

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

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

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

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

5 結(jié) 語

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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