復(fù)雜隧道模型的建模技術(shù)

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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)編號坐標(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í)例

武漢軌道交通三號線越江區(qū)間隧道，區(qū)間線路出宗關(guān)站后，沿建一路南側(cè)地塊向西南方向延伸，下穿水廠路中學(xué)、過沿河大道后，區(qū)間隧道下穿漢江進(jìn)入漢陽區(qū)，下穿江漢二橋體育訓(xùn)練基地和龍陽大道與琴臺大道交叉口后，線路轉(zhuǎn)入龍陽大道下繼續(xù)向西南方向延伸，左右線雙繞龍陽大道高架橋樁基，最后到達(dá)漢陽大道與龍陽大道交叉口的王家灣站.項目段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ù)量和計算時間，取前60 m為研究對象.依照前面所述的建模方法，可以方便快捷的建立三維模型及網(wǎng)格劃分，如圖5所示.該網(wǎng)格由4面體和6面體組成，總共47 360個單元，50 840個節(jié)點(diǎn)，從ANSYS導(dǎo)出到利用FORTRAN轉(zhuǎn)換，整個過程僅用了10 s不到.

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

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

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

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

4 結(jié) 語

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

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

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

致 謝

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

參考文獻(xiàn)：

[1] 王保旗. FORTRAN95程序設(shè)計與數(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ù)化編程與命令手冊[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].北京：中國水利水電出版社,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é)報，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] 董鵬, 李長洪. 基于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].中國地質(zhì)災(zāi)害與防治學(xué)報, 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] 劉麗賢, 馬國鷺, 趙登峰. 基于APDL的ANSYS網(wǎng)格劃分應(yīng)用[J]. 重慶科技學(xué)院學(xué)報, 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] 袁國勇. ANSYS網(wǎng)格劃分方法的分析[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)四，魏化中.快開式壓力容器異型密封圈有限元分析[J].武漢工程大學(xué)學(xué)報，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è)計中的應(yīng)用[J].武漢理工大學(xué)學(xué)報，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)