基于ANSYS的波動(dòng)有限元程序前處理

基于ANSYS的波動(dòng)有限元程序前處理研究摘要：波動(dòng)有限元能夠真實(shí)反映地震動(dòng)傳播特性，因而逐漸成為地震工程和巖土工程重要的數(shù)值模擬方法之一，但現(xiàn)有波動(dòng)有限元程序建模難度大，缺乏方便而快捷的前處理程序，使其在科研工作和實(shí)際工程中的應(yīng)用受到了很大的限制。為此，提出了一套快速易學(xué)的前處理思路和方法，即以ANSYS有限元程序完成的復(fù)雜地質(zhì)體網(wǎng)格劃分為基礎(chǔ)，采用VisualBasic語言編寫了波動(dòng)有限元程序的前處理數(shù)據(jù)轉(zhuǎn)換程序和波動(dòng)有限元程序控制性參數(shù)的可視化輸入界面，使波動(dòng)有限元程序前處理變得輕松自如，簡單易用。經(jīng)實(shí)例驗(yàn)證，該方法快速高效，對(duì)其它程序的前處理具有重要的借鑒意義。關(guān)鍵詞：波動(dòng)有限元；復(fù)雜模型；前處理程序；自動(dòng)生成引言我國學(xué)者廖振鵬教授將時(shí)域集中質(zhì)量有限單元法和人工透射邊界結(jié)合起來給出一種波動(dòng)問題的時(shí)域有限元解法[1，2]，其方法解波動(dòng)過程最為直接，可以很好地模擬地震動(dòng)的傳播過程，因而逐漸受到了科研人員和工程人員的高度重視和廣泛應(yīng)用。同時(shí)該方法在實(shí)現(xiàn)的過程中采用了顯式解法，降低了對(duì)計(jì)算機(jī)硬件的要求，因而受到大家的歡迎。楊柏坡主持編寫的“二維顯式有限元波動(dòng)程序包LAND”于1992年通過部級(jí)鑒定，為科研和工程應(yīng)用奠定了堅(jiān)實(shí)的基礎(chǔ)[3]，周正華[4]、李山有[5]等也對(duì)該程序的功能進(jìn)行了擴(kuò)展和改進(jìn)。需要指出的是，原有的計(jì)算程序均采用FORTRAN77語言，數(shù)組采用固定數(shù)組，為了能夠充分利用計(jì)算機(jī)內(nèi)存，不得不結(jié)合解題規(guī)模進(jìn)行編譯，使用起來比較麻煩。為了能夠使波動(dòng)有限元程序自動(dòng)結(jié)合解題規(guī)模充分利用內(nèi)存，作者采用具有先進(jìn)的FORTRAN95對(duì)上述程序進(jìn)行了改造升級(jí)，采用了動(dòng)態(tài)數(shù)組，使其能夠高效利用計(jì)算機(jī)內(nèi)存，避免了多次編譯，同時(shí)將程序重新命名為WPA2D(2-DimensionWavePropagationAnalysis)，并重新編寫兩側(cè)邊界不等高的邊界處理程序、增加了彈粘塑性分析、粘彈性非線性等功能，作者將另文介紹。盡管主分析程序的計(jì)算穩(wěn)定性和功能已經(jīng)達(dá)到了實(shí)用化程序，但現(xiàn)有程序在模型建立及數(shù)據(jù)準(zhǔn)備等前處理問題上卻存在以下不足，使其在科研工作和實(shí)際工程中的應(yīng)用受到了很大的限制：（1）模型的建立只能靠數(shù)據(jù)文件來實(shí)現(xiàn)，不是很直觀，不能像ANSYS或ALGOR等有限元軟件，可以直接進(jìn)行圖形的處理；（2）對(duì)于比較復(fù)雜的工程模型，在建模時(shí)需要各控制點(diǎn)詳細(xì)的數(shù)據(jù)，容易出錯(cuò)，檢查起來也不是很容易；（3）建模工作量大，花費(fèi)時(shí)間長，直接造成了模擬計(jì)算的周期長、難度大。為解決顯式波動(dòng)有限元程序建模的不足，楊柏坡等用Fortran語言編寫了顯式波動(dòng)有限元的前處理程序[3]，對(duì)于水平成層的場地實(shí)現(xiàn)了快速、便捷的建模。吳兆營編寫了基于ANSYS專門針對(duì)土石壩的前處理程序[6]，該程序只能針對(duì)兩側(cè)邊界等高的情況，但是對(duì)于具有復(fù)雜的地質(zhì)結(jié)構(gòu)等地質(zhì)體，其顯式波動(dòng)有限元模型的建立及網(wǎng)格劃分仍然非常不便。然而，許多有限元程序?qū)τ趶?fù)雜工程地質(zhì)體數(shù)值模擬的建立有明顯優(yōu)勢。由于研究對(duì)象的單元節(jié)點(diǎn)信息，與大型有限元軟件完全一致，即將研究對(duì)象劃分為具有一定形狀且有一定編制規(guī)則的節(jié)點(diǎn)組成的單元面。因此，作者利用通用的建模與網(wǎng)格劃分功能強(qiáng)大的ANSYS軟件對(duì)復(fù)雜工程地質(zhì)體建立相應(yīng)數(shù)值模型（包括網(wǎng)格劃分），再通過數(shù)據(jù)轉(zhuǎn)換實(shí)現(xiàn)WPA2D前處理的自動(dòng)生成。WPA2D程序的前處理目前，WPA2D只包含平面四節(jié)點(diǎn)四邊形，對(duì)應(yīng)于ANSYS有限元程序中的PLANE42單元類型，剖分時(shí)，直接選取該單元類型。根據(jù)波動(dòng)分析要求，計(jì)算出單元所需的尺寸界限，以保證波動(dòng)分析的有效性和穩(wěn)定性，并完成復(fù)雜地質(zhì)體的建模、網(wǎng)格劃分。由于WPA2D的計(jì)算數(shù)據(jù)與ANSYS前處理所得到的數(shù)據(jù)存在一定的差異，因此，利用ANSYS所得的數(shù)據(jù)轉(zhuǎn)換成WPA2D模型的思路和接口程序是本文建議建模方法的關(guān)鍵，實(shí)施原則是最大程度地發(fā)揮ANSYS數(shù)據(jù)的功效，盡量減少人工輸入。其實(shí)施的具體步驟如下：ANSYS模型的建立、ANSYS和WPA2D的數(shù)據(jù)轉(zhuǎn)換和波動(dòng)控制性參數(shù)的輸入。1.1ANSYS模型的建立ANSYS作為通用的有限元軟件，充分綜合了CAD，CAE，CAM等圖象處理工具[7]，是建立復(fù)雜計(jì)算模型有效而又方便快捷的平臺(tái)。ANSYS可以自上而下直接建立實(shí)體模型，還可通過自下而上依次生成點(diǎn)、線、面和體，從而建立實(shí)體模型。強(qiáng)大的布爾運(yùn)算工具可以實(shí)現(xiàn)實(shí)體之間加、減、分類、搭接、粘接和分割等復(fù)雜運(yùn)算，大大提高了建立復(fù)雜地質(zhì)體模型的效率。對(duì)于實(shí)體模型的網(wǎng)格剖分，ANSYS提供了功能強(qiáng)大的控制工具，如單元大小和形狀的控制、網(wǎng)格的劃分類型（自由和映射）以及網(wǎng)格的清除和細(xì)化；此外，還可對(duì)實(shí)體模型圖直接劃分網(wǎng)格；最后，輸出各單元節(jié)點(diǎn)坐標(biāo)及單元信息文件NLIST.LIS和ELIST.LIS。1.2WPA2D與ANSYS數(shù)據(jù)的關(guān)系要將ANSYS前處理所得的單元信息和節(jié)點(diǎn)信息準(zhǔn)確轉(zhuǎn)換為WPA2D所需的數(shù)據(jù)，必須分析兩者之間的相同與不同之處。ANSYS前處理所得的節(jié)點(diǎn)信息包括節(jié)點(diǎn)編碼和節(jié)點(diǎn)坐標(biāo)，單元信息包括單元材料序號(hào)、單元類型、實(shí)常數(shù)、單元系統(tǒng)的坐標(biāo)指針、剖面代號(hào)和單元節(jié)點(diǎn)編碼[6]。除了單元類型、實(shí)常數(shù)、單元系統(tǒng)的坐標(biāo)指針、剖面代號(hào)外，其它信息也是WPA2D所要求的數(shù)據(jù)。與ANSYS數(shù)據(jù)的差別還在于WPA2D將節(jié)點(diǎn)分為內(nèi)節(jié)點(diǎn)和邊界節(jié)點(diǎn)，且邊界節(jié)點(diǎn)進(jìn)一步細(xì)分為左邊界節(jié)點(diǎn)、右邊界節(jié)點(diǎn)和底邊界節(jié)點(diǎn)，并且需要提供邊界區(qū)節(jié)點(diǎn)編碼及邊界節(jié)點(diǎn)所處的層數(shù)（自底邊界算起，底邊界為1），邊界節(jié)點(diǎn)處的人工波速，還有左右邊界節(jié)點(diǎn)縱坐標(biāo)信息（從小到大）及對(duì)應(yīng)節(jié)點(diǎn)處的材料編碼，供邊界區(qū)自由場計(jì)算使用。仔細(xì)總結(jié)分析WPA2D所需數(shù)據(jù)的特點(diǎn)，發(fā)現(xiàn)有以下規(guī)律可以利用：①左、右邊界節(jié)點(diǎn)橫坐標(biāo)為系統(tǒng)橫坐標(biāo)的最小值和最大值，底邊界的縱坐標(biāo)為系統(tǒng)縱坐標(biāo)的最小值，可以充分利用該特點(diǎn)來判斷節(jié)點(diǎn)類型和邊界節(jié)點(diǎn)類型；②左右邊界節(jié)點(diǎn)所在的層數(shù)與其縱坐標(biāo)從小到大所對(duì)應(yīng)的排列序號(hào)相同，可以根據(jù)縱坐標(biāo)利用排序程序?qū)ψ笥疫吔绻?jié)點(diǎn)來進(jìn)行重新排序，從而獲得側(cè)邊界節(jié)點(diǎn)所在的層數(shù)；③與普通有限元的區(qū)別在于，由于采用了透射邊界，要求左右邊界區(qū)至少有3到5個(gè)單元為規(guī)則網(wǎng)格，且單元的水平長度應(yīng)該相等，底邊界至少有3到5個(gè)單元為規(guī)則網(wǎng)格，且單元的豎向長度應(yīng)該相等，根據(jù)以上兩個(gè)特點(diǎn)，首先對(duì)底邊界節(jié)點(diǎn)按橫坐標(biāo)從小到大進(jìn)行重新排序，利用重新排序后的底邊界相臨兩點(diǎn)的橫坐標(biāo)之差的絕對(duì)值給出左右側(cè)邊界區(qū)節(jié)點(diǎn)的水平間隔，同時(shí)利用左右邊界相臨兩點(diǎn)縱坐標(biāo)之差的絕對(duì)值給出底邊界區(qū)的豎向間隔，利用該特點(diǎn)可以找出邊界透射層的節(jié)點(diǎn)編碼。關(guān)于左右邊界節(jié)點(diǎn)處的材料編碼的實(shí)現(xiàn)，首先規(guī)定左右邊界節(jié)點(diǎn)處的材料編碼（自底向上）代表該節(jié)點(diǎn)至上一節(jié)點(diǎn)間的材料號(hào)，除了底部、頂部腳點(diǎn)只包含于一個(gè)單元內(nèi)（可判斷最底部的節(jié)點(diǎn)處的材料號(hào)等于其所在單元的材料號(hào)），其它節(jié)點(diǎn)均包含于兩個(gè)單元，因此節(jié)點(diǎn)材料數(shù)組必須要開一個(gè)臨時(shí)的二維數(shù)組，為了暫時(shí)保存一個(gè)節(jié)點(diǎn)所在的兩個(gè)單元的材料編碼；根據(jù)已排序的邊界節(jié)點(diǎn)，按單元循環(huán)，并將單元材料信息先賦予臨時(shí)的二維數(shù)組，并判斷每個(gè)節(jié)點(diǎn)臨時(shí)兩維數(shù)組與下部相臨節(jié)點(diǎn)的材料號(hào)的異同，如果兩者都相同，即可代表該節(jié)點(diǎn)出的材料號(hào)，否則不同者為節(jié)點(diǎn)處材料號(hào)。具體實(shí)現(xiàn)思路如下：從ANSYS導(dǎo)出輸入單元節(jié)點(diǎn)信息和節(jié)點(diǎn)坐標(biāo)信息，找出，，判斷節(jié)點(diǎn)類型找出邊界節(jié)點(diǎn)編碼和邊界節(jié)點(diǎn)類型，對(duì)左右邊界節(jié)點(diǎn)按縱坐標(biāo)從小到大進(jìn)行排序和底邊界節(jié)點(diǎn)按橫坐標(biāo)進(jìn)行從小到大排序，重新對(duì)邊界節(jié)點(diǎn)進(jìn)行排序，并給出相臨邊界層編碼及層數(shù)，同時(shí)找出左右邊界節(jié)點(diǎn)處的材料號(hào)，輸入材料強(qiáng)度參數(shù)和動(dòng)力計(jì)算控制參數(shù)，依次求出邊界節(jié)點(diǎn)的人工波速（如果節(jié)點(diǎn)上下兩段的材料號(hào)不同時(shí)，取其剪切波速的平均值）。1.3波動(dòng)有限元控制參數(shù)輸入除了利用ANSYS轉(zhuǎn)換的數(shù)據(jù)外，有些波動(dòng)控制參數(shù)尚尚需人工輸入，本文在設(shè)計(jì)程序界面時(shí)，進(jìn)行提供交互式輸入，減少人工勞動(dòng)，其輸入的參數(shù)見圖1和圖2。在上述基礎(chǔ)上，采用VisualBasic語言編寫了波動(dòng)有限元程序的前處理數(shù)據(jù)轉(zhuǎn)換和波動(dòng)有限元程序控制性參數(shù)的可視化輸入界面程序-PreWPA2D（圖1和圖2），并按照程序WPA2D的輸入數(shù)據(jù)格式保存到相應(yīng)的數(shù)據(jù)文件，完成波動(dòng)有限元程序WPA2D前處理工作。圖1波動(dòng)控制性參數(shù)輸入和轉(zhuǎn)換程序界面應(yīng)用實(shí)例本文選取一個(gè)比較復(fù)雜的且具有代表性的含軟弱夾層的巖質(zhì)邊坡作為示例，其模型尺寸如圖2所示。對(duì)于這樣一個(gè)地表形態(tài)和地質(zhì)結(jié)構(gòu)復(fù)雜的邊坡，采用本文提出的前處理方法很理想地建立其計(jì)算模型，具體操作如下：圖2某邊坡地質(zhì)模型示意圖（單位：m）（1）ANSYS網(wǎng)格剖分及數(shù)據(jù)輸出。首先，找出邊坡的幾何控制點(diǎn)作為ANSYS文件的關(guān)鍵點(diǎn)輸入，將其連線，生成實(shí)體模型，定義材料類型，選擇PLANE42單元類型；然后根據(jù)波動(dòng)計(jì)算要求控制單元尺寸，用MESHTOOL命令完成網(wǎng)格剖分，并用COMPRESSNUMBER命令使節(jié)點(diǎn)和單元編號(hào)達(dá)到最優(yōu)，避免冗余，模型劃分3578個(gè)四邊形單元，3726個(gè)節(jié)點(diǎn)；最后用LIST命令欄中的NODES…和ELEMENTS命令輸出節(jié)點(diǎn)信息文件NLIST.LIS和單元節(jié)點(diǎn)信息文件ELIST.LIS。（2）運(yùn)行PreWPA2D。根據(jù)PreWPA2D提示信息，依次完成所有的選擇和輸入，并點(diǎn)擊命令“開始轉(zhuǎn)換”，程序自動(dòng)完成前處理工作，完成后相應(yīng)的命令欄會(huì)出現(xiàn)提示性信息“完成轉(zhuǎn)換”，關(guān)閉該程序，即可運(yùn)行程序WPA2D進(jìn)行計(jì)算，生成的數(shù)值模型如圖3所示。圖3某邊坡數(shù)值模型（單位：m）整個(gè)建模從數(shù)據(jù)采集到模型開始計(jì)算，共用時(shí)不到1個(gè)小時(shí)，目前這是其它方法無可比擬的。3結(jié)論本文提出了基于ANSYS平臺(tái)復(fù)雜地質(zhì)體波動(dòng)有限元程序WPA2D數(shù)值模型的自動(dòng)生成方法。該方法不僅克服了建立復(fù)雜計(jì)算模型的困難，而且實(shí)現(xiàn)了建模的自動(dòng)化。其流程的每一個(gè)操作均簡單、易行，所要求的基本數(shù)據(jù)僅為幾何圖形數(shù)據(jù)點(diǎn)，大大減少了建模所需的時(shí)間、精力，提高了數(shù)值模擬的頻度，有利于充分發(fā)揮WPA2D解決科研和實(shí)際工程問題的強(qiáng)大功能，也有助于WPA2D在地震工程和巖土工程領(lǐng)域的普及和應(yīng)用。該方法的主要思路著眼于發(fā)揮和結(jié)合WPA2D與ANSYS軟件的優(yōu)點(diǎn)，并開發(fā)了WPA-ANSYS接口程序PreWPA2D，解決了復(fù)雜工程地質(zhì)體的建模難題。對(duì)于復(fù)雜工程地質(zhì)問題的數(shù)值模擬，本文提出的方法值得借鑒。參考文獻(xiàn)：[1]廖振鵬.工程波動(dòng)理論導(dǎo)論[M].北京:科學(xué)出版社,2002.[2]楊柏坡,陳慶彬.顯式有限元法在地震工程中的應(yīng)用[J].世界地震工程,1992,8(4):31-40.[3]董瑩.“二維顯式有限元波動(dòng)程序包”科技成果在哈爾濱通過部級(jí)鑒定[J].世界地震工程,1993,9(1):43.[4]周正華.地基土介質(zhì)阻