有限單元法綜述

有限單元法發(fā)展與趨勢綜述有限單元法發(fā)展與趨勢綜述摘要：本文簡述了有限單元法的提出過程，發(fā)展歷史及基本思想，重點論述了有限單元法的基本解題步驟及有限單元法的發(fā)展趨勢。文章最后介紹了一個基于ANSYS軟件運用有限單元法分析某隧道的實例關鍵詞：有限單元法；發(fā)展趨勢；ANSYS；隧道ReviewoftheFiniteElementMethod’sDevelopmentsandTrendsAbstract：Thispaperbrieflydiscussestheproposedprocessofthefiniteelement(FEM),thedevelopmentoftheFEMandthebasicideaoftheFEM,withemphasisonthebasicproblem-solvingstepsoftheFEMandthedevelopmenttrendoftheFEM.Finally,thepaperanalyzesatunnelproblembyusingtheFEM,whichisbasedontheANSYSsoftwareplatform.Keywords：thefiniteelement；developmenttrends；ANSYS；tunnel引言隨著現代工業(yè)、生產技術的發(fā)展，不斷要求設計高質量、高水平的大型、復雜和精密的機械及工程結構。為此目的，人們必須預先通過有效的計算手段，確切的預測即將誕生的機械和工程結構，在未來工作時所發(fā)生的應力、應變和位移。但是傳統(tǒng)的一些方法往往難以完成對工程實際問題的有效分析。彈性力學的經典理論，由于求解偏微分方程邊值問題的困難，只能解決結構形狀和承受載荷較簡單的問題，對于幾何形狀復雜、不規(guī)則邊界、有裂縫或厚度突變，以及幾何非線性、材料非線性等問題往往遇到很多麻煩，試圖按經典的彈性力學方法獲得解析解是十分困難的，甚至是不可能的。因此，需要尋求一種簡單而又精確的數值分析方法。有限單元法正是適應這種要求而產生和發(fā)展起來的一種十分有效的數值計算方法。有限單元法概述有限單元法早在40年代初期就已有人提出，但當時由于沒有計算工具而擱置，一直到50年代中期，高速數字電子計算機的出現和發(fā)展為有限單元法的應用提供了重要的物質條件，才使有限單元法得以迅速發(fā)展。有限單元法在西方起源于飛機和導彈的結構設計，發(fā)表這方面文章最早而且最有影響的是西德的J.H.Argyris教授，于1954–1955年間，他在《Aircraftengineering》上發(fā)表了許多有關這方面的論文，并在此基礎上寫成了《能量原理與結構分析》，此書成為有限單元法的理論基礎。美國的M.T.Turner，R.W.Clough,H.C.Martin和L.J.Topp等人于1956年發(fā)表了一篇題為《復雜結構的剛度和撓度分析》一文，此文提出了計算復雜結構剛度影響系數的方法，說明了如何利用計算機進行分析。美國教授R.W.Clough于1960年在一篇介紹平面應力分析的論文中，首次提出了有限單元法的名字。1965年英國的Zienliewice教授及其合作者解決了將有限元應用于所有場的問題，使有限單元法的應用范圍更加廣泛。有限單元法的優(yōu)點很多，其中最突出的優(yōu)點是應用范圍廣。發(fā)展至今，不僅能解決靜態(tài)的、平面的、最簡單的桿系結構，而且還可以解決空間問題、板殼問題、結構的穩(wěn)定性問題、動力學問題、彈塑性問題和粘彈性問題、疲勞和脆性斷裂問題以及結構的優(yōu)化設計問題。而且不論物體的結構形式和邊界條件如何復雜，也不論材料的性質和外載荷的情況如何，原則上都能應用。有限單元法的基本思想有限單元法的基本思想，是在力學模型上將一個原來連續(xù)的物體離散成為有限個具有一定大小的的單元，這些單元僅在有限個節(jié)點上相連接，并在節(jié)點上引進等效力以代替實際作用于單元上的外力。對于每個單元，根據分塊近似的思想，選擇一種簡單的函數來表示單元內位移的分布規(guī)律，并按彈性理論中的能量原理（或用變分原理）建立單元節(jié)點力和節(jié)點位移之間的關系。最后，把所有的單元的這種關系式集合起來，就得到一組以節(jié)點位移為未知量的代數方程組，解這些方程組就可以求出物體上有限個離散節(jié)點上的位移。有限單元法的基本步驟有限單元法分析計算的基本步驟可歸納為以下五點：1.結構的離散化結構的離散化是有限單元法分析的第一步，它是有限元法的基礎。將某個機械結構劃分為由各種單元組成的計算模型，這一步稱作單元劃分。離散后單元與單元之間利用單元的節(jié)點相互連接起來，將求解區(qū)域變成為用點、線或面劃分的有限數目的單元組合成的集合體。單元的形狀原則上是任意的。例如，在平面問題中通常采用三角形單元，有時也采用矩形或任意四邊形單元。在空間問題中，可以采用四面體、長方體或任意六面體單元?？梢?，不管單元取什么樣的形狀，在一般情況下，單元的便捷總不可能與求解區(qū)域的真實邊界完全吻合，這就帶來了有限元法的一個基本近似性——幾何近似。在一個具體的結構中，確定單元的類型和數目以及哪些部位的單元可以取得大一些，哪些部位單元應該取得小一些，需要由經驗來做出判斷。單元劃分越細則描述變形情況越精確，即越接近實際變形，但計算量越大。所以有限元法中分析的結構已不是原有的物體或結構物，而是同樣材料的眾多單元以一定方式連接成的離散物體。這樣，用有限元分析計算所得的結果只是近似的。如果劃分單元數目非常多而又合理，則所獲得的計算結果就越逼近實際情況。1)網格劃分基礎與劃分原則復雜結構的離散是有限元分析的基礎，也決定著計算結果的精確度。一個復雜的結構總可以離散為一維、二維、三維的小單元。當然對二維和三維單元，其離散后的形狀可以為任意的，但是為了計算的方便性和精確性的結合，二維單元一般采用三角形和四邊形，而三維單元則采用四面體和六面體。簡單的說，復雜結構的離散就是網格的劃分。有限元網格[3]的劃分有很多原則，一是網格數量，網格數量直接影響計算精度和計算時耗,網格數量增加會提高計算精度,但同時計算時耗也會增加。當網格數量較少時增加網格計算精度可明顯提高,但計算時耗不會有明顯增加;當網格數量增加到一定程度后,再繼續(xù)增加網格時精度提高就很小,而計算時耗卻大幅度增加。所以在確定網格數量時應權衡這兩個因素綜合考慮。二是網格密度，為了適應應力等計算數據的分布特點,在結構不同部位需要采用大小不同的網格。如在孔的附近有集中應力，因此網格需要加密，周邊應力梯度相對較小，網格劃分較稀。該網格反映了疏密不同的網格劃分原則：在計算數據變化梯度較大的部位。為了較好地反映數據變化規(guī)律，需要采用比較密集的網格；而在計算數據變化梯度較小的部位，為減小模型規(guī)模，網格則應相對稀疏。三是單元階次，單元階次與有限元的計算精度有著密切的關聯，單元一般具有線性、二次和三次等形式，其中二次和三次形式的單元稱為高階單元。高階單元的曲線或曲面邊界能夠更好地逼近結構的曲線和曲面邊界，且高次插值函數可更高精度地逼近復雜場函數，所以增加單元階次可提高計算精度。但增加單元階次的同時網格的節(jié)點數也會隨之增加，在網格數量相同的情況下由高階單元組成的模型規(guī)模相對較大,，因此在使用時應權衡考慮計算精度和時耗。四是網格形狀，網格單元形狀的好壞對計算精度有著很大的影響，單元形狀太差的網格甚至會中止計算。在網格劃分時應保證合理的單元形狀，即使只有一個單元形狀很差或畸形時,也可能給計算結果帶來很大的誤差，甚至使得計算無法進行下去。2）對網格的評價單元形狀評價一般有以下幾個指標：(1)單元的邊長比、面積比或體積比以正三角形、正四面體、正六面體為參考基準，理想單元的邊長比為一,線性單元可接受的邊長比小于三,二次單元小于十。(2)扭曲度:單元面內的扭轉和面外的翹曲程度。(3)節(jié)點編號:節(jié)點編號對于求解過程中總剛矩陣的帶寬和波前因數有較大的影響,從而影響計算時耗和存儲容量的大小。因此合理的節(jié)點編號有利于剛度矩陣對稱、帶狀分布等求解效率,從而提高計算速度。3）不同維數模型劃分介紹我們對各維模型的單元劃分做簡要的介紹。一維單元可分為兩種。一類是單元的節(jié)點參數中只包含場函數的節(jié)點值C0型，另一類是單元的節(jié)點參數中，除場函數的結點值外，還包含場函數導數的節(jié)點值的C1型單元。這分別是拉格朗日單元和Hermite單元。也就是說拉格朗日是一次插值單元，而后者是二次插值，這樣就能保證導數的連續(xù)性，也就是能保證在連接處除了位移連續(xù)，連接的交點也是光滑的。對二維單元，可以采用三角形和四邊形單元。對三角形單元，如同一維單元的情形，可以利用總體笛卡爾坐標，也可以利用無量綱的局部自然坐標以構造三角形單元的插值函數。利用總體笛卡爾坐標構造三結點三角形單元的差值函數較復雜，更普遍采用的是局部自然坐標來直接構造一般三角行單元的差值函數，這時運算比較簡單。三角形單元的插值一般采用面積坐標，把一個三角形用線段分成等分塊，由插值函數的性質等可以推導出差值函數。通常情況下，采用矩形單元比三角形單元更為方便而有效。其差值函數的推導和一維情況也很相似，也可以構造二維的拉格朗日矩形單元和Hermite矩形單元。此時后者的精度同樣比拉格朗日單元的精度要高。由于有時四邊形單元的節(jié)點在矩形內部，所以一個偶然的發(fā)現，Serendipity四邊形單元被發(fā)現，這個單元有很多優(yōu)勢，一方面由于在實際用用中優(yōu)勢希望統(tǒng)一單元的不同邊界有不同數目的節(jié)點，這樣可以實現不同階次單元之間的過渡，從而可能在求解的不同區(qū)域采用不同精度的單元，另一方面通過它闡述構造單元插值函數的一般方法。三維單元可能有的幾何形狀要比二維單元多得多，在應用中只討論幾種常用的形狀，又因為構造其插值函數的方法只是二維的推廣，所以其形式是很容易構造出來的。其四面體單元也可以用體積坐標，同時也存在Serendipity單元。2.單元分析1）選擇位移模式位移模式是表示單元內任意點的位移隨位置變化的函數式，由于所采用的函數是一種近似的試函數，一般不能精確地反映單元中真實的位移分布，這就帶來了有限元法的另一種基本近似性。采用位移法時，物體或結構物離散化之后，就可把單元中的一些物理量如位移、應變和應力等由節(jié)點位移來表示。這時可以對單元中位移的分布采用一些能逼近原函數的近似函數予以描述。通常，有限單元法中我們將位移表示為坐標變量的簡單函數，這種函數稱為位移模式或位移函數，如式中，是待定系數；是與坐標有關的某種函數。2）建立單元剛度矩陣選定單元的類型和位移模式以后，就可按虛功原理或最小勢能原理建立單元剛度方程，它實際上是單元各個節(jié)點的平衡方程，其系數矩陣稱為單元剛度矩陣式中，為單元編號；為單元的節(jié)點位移向量；為單元的節(jié)點力向量；為單元剛度矩陣，它的每一個元素都反映了一定的剛度特性。根據單元的材料性質、形狀、尺寸、節(jié)點數目、位置及其含義等，找出單元節(jié)點力和節(jié)點位移的關系式，這是單元分析中的關鍵一步。此時需要應用彈性力學中的幾何方程和物理方程來建立力和位移的方程式，從而導出單元剛度矩陣，這是有限元法的基本步驟之一。3）計算等效節(jié)點力物體離散化后，假定力是通過節(jié)點從一個單元傳遞到另一個單元。但是，對于實際的連續(xù)體，力是從單元的公共邊界傳遞到另一個單元中去的。因而，這種作用在單元邊界的表面力、體積力或集中力都需要等效的移到節(jié)點上去，也就是用等效的節(jié)點力來代替所有作用在單元上的力。3.整體分析有限元法的分析過程就是先分后合，即先進行單元分析，在建立了單元剛度方程以后，再進行整體分析，把這些方程集成起來，形成求解區(qū)域的剛度方程，稱為有限元位移法基本方程。集成所遵循的原則是各相鄰單元在共同節(jié)點處具有相同的位移。利用結構力的平衡條件和邊界條件把各個單元按原來的結構重新連接起來；形成整體的有限元方程式中，K為整體結構的剛度矩陣；為整體節(jié)點位移向量；F為整體載荷向量。4.求解方程，得出節(jié)點位移解有限元方程式得出位移、這里，可以根據方程組的具體特點來選擇合適的計算方法。有限元方程是一組方程組，方程在經歷平衡問題中就是以節(jié)點位移為基本未知量的系統(tǒng)結點平衡方程。有限元求解的效率及計算結果的精確很大成都上取決于線性代數方程組的解法。特別是隨著研究對象的更加復雜，有限元分析需要采用更多單元的離散模型來近似實際結構或力學問題的幾何構形時，線性代數方程租的階數就愈來愈高。因而，線性方程組采用何種有限的方法求解，以保證求解的效率和精度就稱為更加重要的問題。不僅在線性靜力分析中，求解代數方程組的時間在整個解題時間中占有很大比重，而且在動力分析和非線性分析中這部分比重也是相當大的。若不采用適當的求解方法，不僅計算費用大量增加，嚴重時可能導致求解過程的不穩(wěn)定和求解的失敗。線性代數方程組的解法可以分為兩大類，即直接解法和迭代解法。直接解法的特點是，選定某種形式的直接解法以后，對于一個給定的線性代數方程組，事先可以按規(guī)定的算法步驟計算出它所需要的算數運算操作數，直接給出最后的結果。迭代解法的特點是，對于一個給定的線性代數方程租，首先假設一個初始解，然后按一定的算法共識進行迭代。在每次迭代過程中對解的誤差進行檢查，并通過增加迭代次數不斷降低解的誤差，直至滿足解的精度要求，并輸出最后的解答。迭代解法的優(yōu)點之一是，它不要求保存洗漱矩陣中高度輪廓線以下的零元素，并且不對它們進行運算，即它們保持為零不變。這樣一來，計算機只需存儲洗漱矩陣的非零元素以及記錄它們位置的輔助數組。這不僅可以最大限度地節(jié)約了存儲空間，而且提高了計算效率。另一方面，迭代解法在計算過程中可以對解的誤差進行檢查，并通過增加迭代次數來降低誤差，直至滿足解的精度要求。其不足之處是，每一種迭代算法可能只適合某一類問題，常缺乏通用的有效性，如使用不當，可能會出現迭代收斂很慢，甚至不收斂的情況。5.由節(jié)點位移計算單元的應變與應力解出節(jié)點位移以后，根據需要，可由彈性力學的幾何方程和物理方程來計算應變和應力。通過上述分析，可以看出，有限單元法的基本思想是“一分一合”，化整為零，集零為整，把復雜的結構看成由有限個單元組成的整體。有限單元法的發(fā)展趨勢有限單元法已經成為現代力學領域分析問題的一個最重要的途徑，為了方便用戶的使用和適應問題復雜性的要求，目前有限單元法發(fā)展方向主要集中在以下幾個方面：更為強大的網格處理能力有限元法求解問題的基本過程主要包括：分析對象的離散化、有限元求解、計算結果的后處理三部分。由于結構離散后的網格質量直接影響到求解時間及求解結果的正確性與否，近年來各軟件開發(fā)商都加大了其在網格處理方面的投入，使網格生成的質量和效率都有了很大的提高，但在有些方面卻一直沒有得到改進，如對三維實體模型進行自動六面體網格劃分和根據求解結果對模型進行自適應網格劃分，除了個別商業(yè)軟件做得較好外，大多數分析軟件仍然沒有此功能。自動六面體網格劃分是指對三維實體模型程序能自動的劃分出六面體網格單元，現在大多數軟件都能采用映射、拖拉、掃略等功能生成六面體單元，但這些功能都只能對簡單規(guī)則模型適用，對于復雜的三維模型則只能采用自動四面體網格劃分技術生成四面體單元。對于四面體單元，如果不使用中間節(jié)點，在很多問題中將會產生不正確的結果，如果使用中間節(jié)點將會引起求解時間、收斂速度等方面的一系列問題，因此人們迫切的希望自動六面體網格功能的出現。自適應性網格劃分是指在現有網格基礎上，根據有限元計算結果估計計算誤差、重新劃分網格和再計算的一個循環(huán)過程。對于許多工程實際問題，在整個求解過程中，模型的某些區(qū)域將會產生很大的應變，引起單元畸變，從而導致求解不能進行下去或求解結果不正確，因此必須進行網格自動重劃分。自適應網格往往是許多工程問題如裂紋擴展、薄板成形等大應變分析的必要條件。由單一結構場求解發(fā)展到耦合場問題的求解有限元分析方法最早應用于航空航天領域，主要用來求解線性結構問題[5]，實踐證明這是一種非常有效的數值分析方法。而且從理論上也已經證明，只要用于離散求解對象的單元足夠小，所得的解就可足夠逼近于精確值。現在用于求解結構線性問題的有限元方法和軟件已經比較成熟，發(fā)展方向是結構非線性、流體動力學和耦合場問題的求解。例如由于摩擦接觸而產生的熱問題，金屬成形時由于塑性功而產生的熱問題,需要結構場和溫度場的有限元分析結果交叉迭代求解，即"熱力耦合"的問題。當流體在彎管中流動時，流體壓力會使彎管產生變形，而管的變形又反過來影響到流體的流動，這就需要對結構場和流場的有限元分析結果交叉迭代求解，即所謂"流固耦合"的問題。由于有限元的應用越來越深入，人們關注的問題越來越復雜，耦合場的求解必定成為有限元的發(fā)展方向。由求解線性問題發(fā)展到求解非線性問題隨著科學技術的發(fā)展，線性理論已經遠遠不能滿足設計的要求，許多工程問題如材料的破壞與失效、裂紋擴展等僅靠線性理論根本不能解決，必須進行非線性分析求解，例如薄板成形就要求同時考慮結構的大位移、大應變（幾何非線性）和塑性（材料非線性）；而對塑料、橡膠、陶瓷、混凝土及巖土等材料進行分析或需考慮材料的塑性、蠕變效應時則必須考慮材料非線性。眾所周知，非線性問題的求解是很復雜的，它不僅涉及到很多專門的數學問題，還必須掌握一定的理論知識和求解技巧，學習起來也較為困難。為此國外一些公司花費了大量的人力和物力開發(fā)非線性求解分析軟件，如ADINA、ABAQUS等。它們的共同特點是具有高效的非線性求解器、豐富而實用的非線性材料庫，ADINA還同時具有隱式和顯式兩種時間積分方法。程序面向用戶的開放性隨著商業(yè)化的提高，各軟件開發(fā)商為了擴大自己的市場份額，滿足用戶的需求，在軟件的功能、易用性等方面花費了大量的投資，但由于用戶的要求千差萬別，不管他們怎樣努力也不可能滿足所有用戶的要求，因此必須給用戶一個開放的環(huán)境，允許用戶根據自己的實際情況對軟件進行擴充，包括用戶自定義單元特性、用戶自定義材料本構、用戶自定義邊界條件、用戶自定義結構斷裂判據和裂紋擴展規(guī)律等等。與圖形軟件如CAD軟件的無縫集成當今有限元分析軟件的一個發(fā)展趨勢是與通用CAD軟件的集成使用，即在用CAD軟件完成部件和零件的造型設計后，能直接將模型傳送到CAE軟件中進行有限元網格劃分并進行分析計算，如果分析的結果不滿足設計要求則重新進行設計和分析，直到滿意為止，從而極大地提高了設計水平和效率。為了滿足工程師快捷地解決復雜工程問題的要求，許多商業(yè)化有限元分析軟件都開發(fā)了和著名的CAD軟件的接口?；贏NSYS的計算實例分析題目:某雙線鐵路隧道，長700米，除進出口洞段（約2100m）巖體為碎石狀松散結構體，巖體級別為V級外，其余500m洞段主要為塊碎石狀鑲嵌結構體，巖體級別為IV級，埋深50m~150m不等?，F根據《鐵路隧道設計規(guī)范》(TB10003—2005)規(guī)定，IV級巖體段隧道襯砌采用復合襯砌，其中的初期支護及二次襯砌采用工程類比法進行設計，并通過數值方法進行驗算。初襯和二襯設計參數見表1，圍巖及支護結構物理力學參數見表2。表1鐵路雙線隧道IV級圍巖復合式襯砌設計參數圍巖級別初期支護二次襯砌(C30)噴射混凝土(C20)厚度(cm)錨桿(直徑25mm)鋼筋網鋼架拱、墻(cm)仰拱(cm)拱、墻仰拱位置長度(m)間距(m)IV2020拱、墻3.01.0拱、墻、仰拱@2525—4045表2圍巖及支護結構基本物理力學參數項目密度/kg/m3彈性(變形)模量/GPa泊松比內摩擦角/粘聚力/MPaIV級圍巖22003.20.32350.5C20噴砼2500210.2——C30混凝土2500300.2——錨桿78002000.3——圖1隧道開挖輪廓為兩段圓弧(見圖1)，頂拱及邊墻段圓弧圓心在隧道軸線上，半徑5.5m，圓心角約223；仰拱段圓弧由頂拱及邊墻段圓弧的兩個角點及仰拱頂點組成，仰拱頂點位于頂拱及邊墻段圓弧圓心以下3.8m處。施工采用CRD工法。計算當隧道埋深分別為50m、100m和150m（至隧道軸線）時，初襯內力和圍巖的應力分布特征，并進行比較。計算方法:對于雙線鐵路隧道，采用CRD工法施工，其襯砌包括初期支護和二次襯砌。本次計算中采用開挖數值模擬（地層結構模型）來分別計算隧道埋深為50m，100m，150m時初襯內力和圍巖的應力分布。地層結構模型需要用到ansys軟件中平面單元（Plane42）模擬圍巖和桿單元（Link）模擬錨桿，采用Beam3單元模擬噴射混凝土和鋼筋網。錨桿布置于邊墻軸線以上的部分，噴混凝土為封閉結構，不考慮鋼筋網對噴砼的加強作用。計算過程:荷載計算：隧道埋深分別為50m，100m，150m時，模擬面均取為50m埋深。100m，150m時模擬面以上的巖體以荷載形式施加，頂面垂直荷載為上部巖體（未模擬部分）自重，側邊水平荷載，按彈性力學方法計算，等于垂直荷載乘側壓力系數μ/(1-μ)，本題中側壓力系數為0.4706。由于該鐵路隧道地下水不豐富，按水土合算側壓力。模擬的地層范圍，橫向兩端延伸至4倍洞徑即44m，高度為埋深50m，隧道底下3倍洞高即31.7m。整個模擬地層為99m×81.7m。（1）．前處理主要包括材料、實常數和單元類型定義，建立幾何模型和單元網格劃分。具體見命令流。生成的關鍵點，生成的線，生成的面，劃分單元大小，單元網格劃分如下列圖所示：關鍵點圖創(chuàng)建的線圖創(chuàng)建的面圖單元大小圖單元網格圖（2）.加載與自重應力場求解Ⅰ.埋深為50米時:輸入相應的命令流加載，可看到