張 濤，黃民水，涂躍亞，周 麟

(1.中交第二公路勘察設(shè)計(jì)研究院有限公司,湖北 武漢 430056；2.武漢工程大學(xué)環(huán)境與城市建設(shè)學(xué)院,湖北 武漢 430074)/

0 引 言

FORTRAN[1]亦譯為福傳，是Formula Translator的縮寫(xiě)，意譯為公式翻譯器，正式發(fā)布于1954年，是世界上最早出現(xiàn)的計(jì)算機(jī)高級(jí)程序語(yǔ)言.它是一種現(xiàn)代的表達(dá)性語(yǔ)言，既適合解決各類(lèi)數(shù)值問(wèn)題，又適合非數(shù)值問(wèn)題，廣泛應(yīng)用于教育部門(mén)和科學(xué)部門(mén)，編寫(xiě)程序簡(jiǎn)潔、高效.ANSYS[2]軟件是一個(gè)強(qiáng)大的通用有限元分析軟件，集結(jié)構(gòu)、熱、電磁、流體、聲學(xué)分析于一體.它不僅兼容性好、計(jì)算能力強(qiáng)、應(yīng)用性廣，而且具有強(qiáng)大的實(shí)體建模功能和多樣的網(wǎng)格劃分手段.FLAC是快速拉格朗日差分分析(Fast Lagrangian Analysis of Continua)的簡(jiǎn)稱(chēng)，F(xiàn)LAC3D[3-4]程序由美國(guó)ITASCA咨詢集團(tuán)公司推出，是目前巖土力學(xué)計(jì)算中的重要數(shù)值方法之一.該程序用于模擬三維土體、巖體及其它材料力學(xué)特性，廣泛應(yīng)用于邊坡、地下洞室、施工過(guò)程(開(kāi)挖、填筑等)、拱壩、隧道和礦山工程等多個(gè)領(lǐng)域的數(shù)值模擬.由于它采用混合離散方法來(lái)模擬材料的屈服或塑性流動(dòng)特性，比有限元方法中采用的降階積分更合理.FLAC3D采用的是顯式方法進(jìn)行求解，它無(wú)需存儲(chǔ)剛度矩陣，采用中等容量的內(nèi)存能求解多單元結(jié)構(gòu)模擬大變形問(wèn)題.即使模擬的系統(tǒng)是靜態(tài)的，仍采用了動(dòng)態(tài)運(yùn)動(dòng)方程，這使得FLAC3D在模擬物理上的不穩(wěn)定過(guò)程不存在數(shù)值上的障礙.然而，F(xiàn)LAC3D在前處理功能較弱，很難建立復(fù)雜三維模型，一般通過(guò)手動(dòng)輸入節(jié)點(diǎn)坐標(biāo)來(lái)建立模型，造成使用者把更多精力浪費(fèi)在了建模上.因此，結(jié)合利用ANSYS強(qiáng)大的建模及網(wǎng)格劃分能力和FORTRAN快速、有效的編程能力，根據(jù)ANSYS與FLAC3D對(duì)應(yīng)的單元和坐標(biāo)關(guān)系，通過(guò)ANSYS建模，然后轉(zhuǎn)換為FLAC3D能識(shí)別的節(jié)點(diǎn)及網(wǎng)格單元文件(*.flac3d)，是一個(gè)能解決FLAC3D建模復(fù)雜問(wèn)題的有效且可行的方法.

1 ANSYS與FLAC3D坐標(biāo)系及模型區(qū)別

ANSYS與FLAC3D單元模型的區(qū)別主要體現(xiàn)在單元節(jié)點(diǎn)編號(hào)[5-6]和坐標(biāo)系的不同，如圖1所示.

圖1 ANSYS與FLAC3D坐標(biāo)系區(qū)別Fig 1 The coordinate system difference between ANSYS and FLAC3D

對(duì)于ANSYS坐標(biāo)系，平面內(nèi)水平向右為x軸正向，豎直向上為y軸正向，垂直平面向外為z軸正向；對(duì)于FLAC3D坐標(biāo)系，平面內(nèi)水平向右為x軸正向，豎直向上為z軸正向，垂直平面向內(nèi)為y軸正向.

ANSYS和FLAC3D的常用單元編號(hào)對(duì)比圖，如圖2所示.FLAC3D的單元節(jié)點(diǎn)編號(hào)排序符合“右手法則”，而ANSYS的單元節(jié)點(diǎn)編號(hào)排序，右下至上螺旋增加.

圖2 單元編號(hào)順序?qū)Ρ葓DFig.2 Comparison chart of element number

2 ANSYS與FLAC3D模型的轉(zhuǎn)換

ANSYS與FLAC3D模型的轉(zhuǎn)換原理就是坐標(biāo)系的轉(zhuǎn)換和單元對(duì)應(yīng)節(jié)點(diǎn)的轉(zhuǎn)換，達(dá)到模型形狀和方向的一致性.而現(xiàn)有的轉(zhuǎn)換程序，有些只簡(jiǎn)單做了單元的轉(zhuǎn)換，未做坐標(biāo)系的轉(zhuǎn)換，轉(zhuǎn)換后的模型方向不一致，需要手動(dòng)修改，給建模造成困難；有些未考慮單元數(shù)量問(wèn)題，造成模型轉(zhuǎn)換失敗或者單元丟失[7-8].因此，編寫(xiě)轉(zhuǎn)換程序必須解決上述問(wèn)題.

轉(zhuǎn)換的具體步驟如下：

1)ANSYS端采用GUI界面操作結(jié)合命令流的方式，利用內(nèi)置強(qiáng)大布爾操作和網(wǎng)格劃分能力，完成模型建立及網(wǎng)格劃分[9-12]，并賦予材料屬性.以兩個(gè)單元為例，單元如圖3所示，左邊為材料類(lèi)型1，右邊為材料類(lèi)型2.

圖3 ANSYS模型圖Fig.3 ANSYS model fig

2)運(yùn)用ANSYS中內(nèi)置的APDL語(yǔ)言編寫(xiě)轉(zhuǎn)換程序，刪除所有area，僅保留element和node信息 ， 然后把節(jié)點(diǎn)總數(shù)與節(jié)點(diǎn)坐標(biāo)和單元總數(shù)與單元數(shù)據(jù)分別寫(xiě)入“node.txt”和“element.txt”兩個(gè)文本文檔中.節(jié)點(diǎn)的存儲(chǔ)格式是“節(jié)點(diǎn)總數(shù)+依序排列的節(jié)點(diǎn)坐標(biāo)”；單元的存儲(chǔ)格式是“單元總數(shù)+材料總數(shù)(轉(zhuǎn)換到FLAC3D中以group形式表示)+組成單元的節(jié)點(diǎn)號(hào)”.

例子中ANSYS節(jié)點(diǎn)坐標(biāo)文件和單元數(shù)據(jù)格式如表1和表2所示.

3)利用FORTRAN95軟件編寫(xiě)轉(zhuǎn)換程序.

由于編寫(xiě)的轉(zhuǎn)換程序需要是一個(gè)通用轉(zhuǎn)換程序，能運(yùn)用到多個(gè)工程項(xiàng)目中，對(duì)于不同項(xiàng)目節(jié)點(diǎn)及單元數(shù)目可能不同，造成程序中記錄節(jié)點(diǎn)和單元的數(shù)組大小是經(jīng)常變化的，如果采用定義固定數(shù)組大小的方法，就會(huì)導(dǎo)致轉(zhuǎn)換過(guò)程中數(shù)組過(guò)小無(wú)法完成，或者數(shù)組設(shè)置過(guò)大造成內(nèi)存的浪費(fèi)影響速度，遇到類(lèi)似情況需要經(jīng)常修改程序.使用可變大小的數(shù)組就可以在程序執(zhí)行時(shí)根據(jù)需要決定數(shù)組大小，節(jié)省內(nèi)存和時(shí)間.

表1 ANSYS節(jié)點(diǎn)坐標(biāo)文件Table 1 The file of ANSYS node coordinates

表2 ANSYS單元數(shù)據(jù)文件

FORTRAN95可以定義可變數(shù)組，需要經(jīng)過(guò)兩個(gè)步驟：第一步，使用allocatable屬性對(duì)相關(guān)數(shù)組進(jìn)行定義；第二步，運(yùn)用allocate語(yǔ)句為該數(shù)組確切的定義大小，分配內(nèi)存空間，然后再使用該數(shù)組.

轉(zhuǎn)換程序主要思路如下：

1)用allocatable屬性定義需要使用的可變數(shù)組，讀取“node.txt”中第一個(gè)數(shù)據(jù)和“element.txt”文件中的第一、二個(gè)數(shù)據(jù)(記錄節(jié)點(diǎn)和單元數(shù))；

2)根據(jù)讀取的數(shù)據(jù)經(jīng)過(guò)相關(guān)轉(zhuǎn)換，使用allocate語(yǔ)句定義數(shù)組大??；

3)由于ANSYS單元都是由8個(gè)節(jié)點(diǎn)組成的，程序根據(jù)8個(gè)節(jié)點(diǎn)中節(jié)點(diǎn)的重復(fù)程度判斷該單元類(lèi)型，然后將“node.txt”和“element.txt”文件中關(guān)于節(jié)點(diǎn)、單元和材料屬性的數(shù)據(jù)根據(jù)節(jié)點(diǎn)關(guān)系和坐標(biāo)關(guān)系轉(zhuǎn)換好后，寫(xiě)入相應(yīng)數(shù)組中；

4)按照FLAC3D模型文件“*.flac3d”的特點(diǎn)將數(shù)組中的數(shù)據(jù)寫(xiě)入文件里.

5)將“*.flac3d”導(dǎo)入FLAC3D后，根據(jù)group名賦予材料屬性，轉(zhuǎn)換后的圖如圖4所示，左邊為group 1，右邊為group 2.例子中關(guān)于“*.flac3d”文件的數(shù)據(jù)如表3和表4所示.

圖4 FLAC3D模型圖Fig.4 FLAC3D model fig

節(jié)點(diǎn)編號(hào)坐標(biāo)XYZG10.1000E+010.0000E+000.1000E+01G20.0000E+000.0000E+000.1000E+01G30.0000E+000.0000E+000.0000E+00G40.1000E+010.0000E+000.0000E+00G50.2000E+010.0000E+000.1000E+01G60.2000E+010.0000E+000.0000E+00G70.1000E+010.1000E+010.1000E+01G80.0000E+000.1000E+010.1000E+01G90.0000E+000.1000E+010.0000E+00G100.1000E+010.1000E+010.0000E+00G110.2000E+010.1000E+010.1000E+01G120.2000E+010.1000E+010.0000E+00

表4 FLAC3D單元數(shù)據(jù)文件Table 4 The file of FLAC3D element datas

3 應(yīng)用實(shí)例

武漢軌道交通三號(hào)線越江區(qū)間隧道，區(qū)間線路出宗關(guān)站后，沿建一路南側(cè)地塊向西南方向延伸，下穿水廠路中學(xué)、過(guò)沿河大道后，區(qū)間隧道下穿漢江進(jìn)入漢陽(yáng)區(qū)，下穿江漢二橋體育訓(xùn)練基地和龍陽(yáng)大道與琴臺(tái)大道交叉口后，線路轉(zhuǎn)入龍陽(yáng)大道下繼續(xù)向西南方向延伸，左右線雙繞龍陽(yáng)大道高架橋樁基，最后到達(dá)漢陽(yáng)大道與龍陽(yáng)大道交叉口的王家灣站.項(xiàng)目段CK11+253.876～CK11+577.50隧道底板埋深30.0～38.5 m左右，標(biāo)高+4.0～-19.5 m左右，基底依次位于(18-2)中風(fēng)化灰?guī)r、中風(fēng)化石英砂巖、中風(fēng)化砂巖、泥質(zhì)砂巖.考慮到單元數(shù)量和計(jì)算時(shí)間，取前60 m為研究對(duì)象.依照前面所述的建模方法，可以方便快捷的建立三維模型及網(wǎng)格劃分，如圖5所示.該網(wǎng)格由4面體和6面體組成，總共47 360個(gè)單元，50 840個(gè)節(jié)點(diǎn)，從ANSYS導(dǎo)出到利用FORTRAN轉(zhuǎn)換，整個(gè)過(guò)程僅用了10 s不到.

圖5 三維網(wǎng)絡(luò)模型Fig.5 Three-dimensional network model

為檢驗(yàn)網(wǎng)格劃分和單元幾何形狀是否滿足要求，進(jìn)行自重應(yīng)力下初試地應(yīng)力場(chǎng)計(jì)算，豎向應(yīng)力云圖如圖6所示，分布均勻且連續(xù)，說(shuō)明網(wǎng)格劃分及模型良好.圖7為隧道開(kāi)挖完成，并完成錨噴支護(hù)后，中間段面(縱向30 m處)的位移圖，拱底因開(kāi)挖卸荷，以隆起為主，最大回彈量為6.9 mm.拱頂因開(kāi)挖發(fā)生沉降，最終穩(wěn)定在10.2 mm.

圖6 豎向應(yīng)力云圖Fig.6 Vertical stress

圖7 豎向位移圖Fig.7 Vertical displacement

4 結(jié) 語(yǔ)

a.與一般轉(zhuǎn)換方法相比較，少了另外修改坐標(biāo)那一步，模型所建即所得.

b.借助FORTRAN高效的計(jì)算能力和獨(dú)特的可變數(shù)組定義功能，可以在節(jié)約內(nèi)存的前提下提高ANSYS模型轉(zhuǎn)換FLAC3D模型的效率.47 360個(gè)單元、50 840個(gè)節(jié)點(diǎn)的模型轉(zhuǎn)換只花了不到10 s的時(shí)間，相比其它方法更為高效.

c.通過(guò)工程實(shí)例驗(yàn)證了通過(guò)ANSYS建模轉(zhuǎn)換FLAC3D計(jì)算的可行性，具有良好的工程應(yīng)用前景.

致 謝

感謝中國(guó)交通第二公路勘察設(shè)計(jì)院為本研究提供支持和幫助！

參考文獻(xiàn)：

[1] 王保旗. FORTRAN95程序設(shè)計(jì)與數(shù)據(jù)結(jié)構(gòu)基礎(chǔ)教程[M]. 天津：天津大學(xué)出版社, 2007.

WANG Baoqi. Foundation Course about Program Design and Data Structure of FORTRAN95[M]. Tianjin: Tianjin University Press, 2007. (in Chinese)

[2] 龔曙光, 謝桂蘭, 黃云清. ANSYS參數(shù)化編程與命令手冊(cè)[M]. 北京：機(jī)械工業(yè)出版社, 2009.

GONG Shu-guang, XIE Gui-lan, HUANG Yun-qing. Parametric Programming and Manual Command of ANSYS[M]. Beijing: China Machine Press, 2009.(in Chinese)

[3] 陳育民, 徐鼎平. FLAC/FLAC3D基礎(chǔ)與工程實(shí)例[M].北京：中國(guó)水利水電出版社,2009.

CHEN Yu-min XU Ding-ping. FLAC/FLAC3D foundation and practical engineering [M]. Beijing: China water conservancy and hydropower press, 2009. (in Chinese)

[4] 歐湘萍，白楷，朱云升.基于FLAC-3D的強(qiáng)度折減法邊坡穩(wěn)定性分析[J].武漢理工大學(xué)學(xué)報(bào)，2008，31(9)：59-61.

OU Xiang-ping,BAI Kai,ZHU Yun-sheng.The Strength Reduction Method for Stability Analysis of Slope Based on FLAC-3D[J].Journal of Wuhan University of Technology,2008,31(9):59-61.(in Chinese)

[5] 彭文斌. FLAC3D實(shí)用教程[M]. 北京：機(jī)械工業(yè)出版社, 2007.

PENG Wen-bin. FLAC3D practical tutorial [M]. Beijing: China Machine Press, 2007.(in Chinese)

[6] 王新敏. ANSYS工程結(jié)構(gòu)數(shù)值分析[M]. 北京：人民交通出版社, 2007.

WANG Xin-min. ANSYS engineering structure numerical analysis [M]. Beijing: China Communications Press, 2007. (in Chinese)

[7] 董鵬, 李長(zhǎng)洪. 基于ANSYS軟件的FLAC3D實(shí)體建模[J]. 山西建筑, 2007, 33(3): 350-351.

DONG Peng,LI Chang-hong. Modeling in FLAC3D software based on ANSYS software entity[J].Shanxi Architecture, 2007, 33(3): 350-351.(in Chinese)

[8] 劉秀軍. 基于GOCAD的復(fù)雜地質(zhì)體FLAC3D模型生成技術(shù)[J].中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào), 2011, 22(4): 41-45.

LIU Xiu-jun. FLAC3D modeling for complex geologic body based on GOCAD [J]. The Chinese Journal of Geological Hazard and Control,2011, 22 (4) : 41-45. (in Chinese)

[9] 劉麗賢, 馬國(guó)鷺, 趙登峰. 基于APDL的ANSYS網(wǎng)格劃分應(yīng)用[J]. 重慶科技學(xué)院學(xué)報(bào), 2008, 10(5): 122 -123.

LIU Li-xian, MA Guo-lu, ZHAO Deng-feng. Grid Division of ANSYS and Application Based on APDL [J]. Journal of Chongqing University of Science and Technology, 2008, 10(5): 122 -123. (in Chinese)

[10] 袁國(guó)勇. ANSYS網(wǎng)格劃分方法的分析[J]. 計(jì)算機(jī)應(yīng)用, 2009(6): 59-60.

YUAN Guo-yong. Analysis of Meshing Division Method [J]. Modern Machinery,2009(6): 59-60. (in Chinese)

[11] 丁克勤，劉關(guān)四，魏化中.快開(kāi)式壓力容器異型密封圈有限元分析[J].武漢工程大學(xué)學(xué)報(bào)，2013，35(8)：52-56.

DING Ke-qin,LIU Guan-si,WEI Hua-zhong.Finite element analysis of special-shaped sealing ring on quick actuating pressuke vessel[J].Journal of Wuhan Institute of technology,2013,35(8):52-56.

[12] 吳作為.ANSYS軟件在盾構(gòu)管片設(shè)計(jì)中的應(yīng)用[J].武漢理工大學(xué)學(xué)報(bào)，2010,37(3)：655-657.

WU Zou-wei.Application of ANSYS Software in shield Segments Pesign[J].Journal of Wuhan University of Technology,2010,37(3):655-657.(in Chinese)