現階段結構力學分析的江湖幾乎已被有限單元法一統江山。但是實際上,有限單元法只是數值方法求解微分方程的眾多手法之一,在此文中,我將介紹現階段幾位一樣能夠參與結構分析華山論劍的絕頂高手。
一、 武林盟主——有限單元法
如今的江湖,各路豪傑並起,熱鬧非凡,上至通用分析軟體中鼎鼎大名的ABAQUS、ANSYS、SAP2000等,下至結構工程師的老朋友PKPM、YJK、Midas,若要追根溯源,有限單元法是他們共同的基礎。在結構設計界,據我所知,有限單元法已經一統江湖,號令天下,莫敢不從,儼然武林盟主之態。那麼這位盟主出身如何呢,且看以下分解。
盟主的身世——有限單元法的發展
有限單元法思想最早出現在1943年,Courant採用定義在一系列三角形區域的分片連續函數結合最小勢能原理求解St.Venant扭轉問題。由於當時的計算工具的限制,這種方法並沒有在接下來的時間裡得到較好的發展。直到1960年以後,電子計算機的廣泛應用和發展,現在意義上的有限單元法才真正發展起來。
現代意義上的有限單元法第一個成功的嘗試出現在1956年,這是由Tuner和Clough等人在分析飛機結構時提出的。他們使用有限單元法三角形單元得出了平面應力問題的正確解答。1960年Clough將其進一步推廣到平面彈性問題,並第一次提出「有限單元法」這一名稱。而在之後的1965年O. C. Zienkiewicz和Y. K. Cheung發現存在變分形式的所有場問題即滿足變分原理的所有場問題,都可以用與固體力學有限元法相同的步驟求解。
但是並不是所有問題都能找到對應的變分原理,在不能找到變分原理的問題中應用有限單元法就必須找到能夠給出有限元格式的方法。在1969年,Oden J T從所謂的能量平衡法出發,得出熱彈性問題的有限元分析方程組。在1969年B. A. Szabo和G. C. Lee將加權餘量法,特別是Galerkin法引入有限單元法導出標準的有限元過程來求解非結構問題。加權餘量法引入有限單元法後使得有限單元法的應用範圍得到了極大拓展。而且有限單元法同門中還有兩位骨骼驚奇的奇俠,一般的有限單元法,形函數都是一階線性的,對於複雜的幾何邊界條件比較吃力,所以有了等參元法和Wilson非協調元,這兩位形函數中含二次項的選手。
20世紀70年代,有限單元法的發展已達到高峰,其基本理論和方法已趨成熟。
盟主的招式——有限單元法的基本步驟
盟主有限單元法在解決敵手時總有個固定的套路。
起手式,將求解域離散為一系列單元。就是現在各路軟體中劃分網格的過程。
緊接一招隔空猜物,就是在各個單元內用基於單元各結點值的單元插值函數來擬合此單元上的未知場函數,就我們常見的力學求解問題,一般都是先求解整體結構的位移函數。
然後一招萬法歸宗,根據力學原理求出各單元的剛度矩陣。再由單元剛度矩陣凝聚為總剛矩陣,配合邊界條件解出數值解。
在一般的結構力學問題中,這幾招下來基本無堅不摧。
盟主的罩門——有限單元法的瑕疵
盟主雖然已然一統江湖,但是其餘高手仍未被除盡,可見盟主武功中也存在弱點。
首先就是單元的劃分,雖然許多有限元程序有了自動劃分單元的功能,但是至今還無法保證劃分的網格一定合理,由於網格劃分不合理導致計算精度差甚至不收斂也不能完全避免。
其次,現有的有限單元法大都採用一階線性形函數,在擬合複雜問題時效率不高,而等參元法擬合的二次形函數並不完備,Wilson非協調元由於不符合變形協調條件,不能保證收斂。
再者,有限單元法在涉及大變形問題時,由於計算中網格會出現嚴重的扭曲,這樣不但需要網格的重構,而且也會導致計算精度的下降。在處理裂紋發展問題時,由於裂紋發展的不確定性,在計算中隨著裂紋的發展不得不不斷地重新劃分單元來模擬整個裂紋發展的動態過程。
二、隱俠—— 無網格法
隱俠的由來——無網格法的發展
針對盟主的罩門,另一位高手修煉了另一門絕學——無網格法,其計算與網格無關。這種無網格法的出現,在大變形和裂縫發展問題中,不像有限單元法一樣需要重新劃分網格,因此無網格法在處理這類問題時不僅能夠提高計算精度,並且還能減小計算難度。
無網格法的出現最早可以追溯到20世紀70年代。當時人們開始研究對非規則網格有限元差分法的研究,這是無網格法的前身。1977年,Lucy提出了「Smoothed Particle Hydrodynamics」方法,此方法即被簡稱為SPH方法。在1995年Liu等人在其基礎上提出了再生核粒子方法——「Reproducing Kernel Particle Method」簡稱為RKPM方法。RKPM方法和SPH方法一樣都是基於核近似的無網格方法,它在SPH方法基礎上引入修正函數施加再生條件並且採用高斯積分。RKPM方法不僅解決了SPH方法在邊界條件上不一致性的問題,而且完全消除了SPH方法的張力不穩定。但是由於RKPM方法不滿足Kronecker delta 函數性質,這使得這種方法施加邊界條件比較困難。
除了基於核近似的無網格法以外,還有採用移動最小二乘法近似的移動最小二乘法——「Moving Least Square」簡稱為MLS方法。在1992年Nayroles等人將Galerkin法引入移動最小二乘法,從而提出了彌散單元法——「Diffuse Element Method」簡稱DEM方法。此後的1994年Belytschko等人對彌散單元法進行了改進,在形函數的導數中保留了被彌散單元法省去的項,並且引入Lagrange乘子法應用邊界條件,這就是無網格Galerkin法——「The element-free Galerkin Method」簡稱為EFG方法。其後又有諸多如有限點法——「The Finite Point Method」簡稱FPM方法、無網格局部Petrov-Galerkin法——Meshless Local Petrov- Garlerkin Method簡稱MLPG方法、再生核階譜單位分解方法等方法。
隱俠的招式——無網格法的基本思想
鄙人對無網格法涉獵不深,從基本思想上來講無網格法是在求解域上布置有限數量的結點,在每個結點上定義覆蓋結點周圍一定區域的結點權函數或稱結點核函數來擬合改結點附近某一區域內的物理量。每個結點所影響的區域稱為其支撐域(支撐域可部分互相重合),由所有結點的支撐域來覆蓋整個求解域。這樣通過各個結點的權函數的疊加就可以擬合待求解函數。
為何歸隱——無網格法的缺點
隱俠無網格法一身奇絕武藝,甚至擺脫了有限單元法網格的桎枯,可以無視網格劃分是否合理達到無招勝有招境界,可是為何現實中卻無一款商業軟體採用此法呢?歸根結底是由於隱俠招法過於奇詭,難以駕馭。就以我們熟悉的結構分析為例,我們已經習慣於在模型的一些的節點上愉快的施加邊界條件。給一個節點剛接約束節點就不能有位移的暢快體驗在無網格法中幾乎成為泡影。因為無網格法中節點的支撐域重合,通常不滿足Kronecker函數特性,即一個節點上的位移並不僅由此節點的形函數確定,周圍的節點都可能會對其有影響,所以就算要施加一個位移為0的約束也是千難萬難,而且無網格法的積分計算較有限元法複雜,所以這位高手暫時只能在特定領域內棲身,江湖中往往僅僅只留下他的傳說。
三、盟主的覬覦者—— 廣義有限單元法
武林中從來不缺驚才絕豔的新人。廣義有限單元法就是這麼一位可能承繼盟主之位的新人。廣義有限單元法最早是在20世紀末由Babuska I等人提出。將傳統有限單元法中的插值函數作為單位分解函數,從中我們能夠看出廣義有限單元法實際上是傳統有限單元法和無網格法中單位分解法結合產生的。廣義有限單元法對比與傳統有限單元法最大的區別就是將能夠反映待求解問題特性的函數引入到插值函數中,從而能避免網格重構,提高計算精度與計算效率。
在大致相同時期,我國也有學者提出了類似的思想。梁國平教授在1995年借鑑石根華的流形思想提出了廣義有限元方法。此後,欒茂田等人引入廣義結點的概念,也就是在傳統有限元方法的結點中增加新的自由度,從而得到不增加結點就可以模擬高階變形的廣義結點有限元方法。辛婭雲在程堯舜教授的指導下改進了欒茂田等人提出的廣義結點有限元方法。具體做法是將各個結點的局部坐標原點設置為改該結點本身,給出了新的插值函數。在2007年,S. Rajendran和B.R. Zhang綜合了最小二乘法、有限單元法和點插值法創造了一種新的四結點四邊形單元,並應用該單元解決了平面應力問題。之後,又將這種單元應用於自由振動分析和非線性柯西材料的平面應力分析。此後,李剛在程堯舜教授的指導下又在辛婭雲的廣義結點有限元方法上進一步改良,引入鄰近點概念結合點插值法提出了半彌散單元法。半彌散單元法是在原有的有限單元法形函數上加入新的自由度使之成為二次完備的形函數,在擬合複雜問題時較傳統有限單元法收斂的更快,精度更高。
初試身手——半彌散單元法的應用
本部分是半彌散單元法展示實力的部分。
本節中的符號約定:
SDEM——表示本文中半彌散單元法;
T3——表示傳統三結點三角形單元的有限單元法。
本節以數值算例比較半彌散單元法與傳統三結點三角形單元在穩態熱傳導問題計算過程中的優劣。
四邊邊界溫度已知的矩形溫度場問題
如下圖所示一無限長矩形柱體的橫斷面,邊長分別為L1與L2,材料性質均勻,且柱體內部無熱源。由於柱體無限長,而邊界上熱交換條件與縱向坐標無關。因此問題簡化為二維穩態熱傳導問題。
該問題的數學模型可以用以下微分方程表示:
邊界條件為
區域內的溫度分布解析解為:
為了方便計算,本文中將L1 、L2都取為1,邊界條件中TW也取作1;同時,材料的熱傳導係數Kx、Ky均取為1;選取計算區域中心位置P(0.5,0.5)為考察點。
根據問題的解析解式(3)可以求出P位置處溫度值應為
半彌散單元法和傳統有限單元法在計算中採用同樣的網格劃分。網格劃分如下圖所示。
採用以上四種網格劃分方法,分別採用半彌散單元法和傳統有限單元法求解該問題,求得的考察點P溫度值與用式(3)求得的考察點P溫度值解析解列於下表。
四種網格劃分方式下相對誤差對比可以下圖表示。
如以自由度數作為比較標準。在表2中列出了半彌散單元法和傳統有限單元法自由度數對比。
從表2可知,半彌散單元法的自由度數與傳統有限單元法的自由度數相同,這是由於將原廣義結點有限元中增加的自由度可用結點信息表示,在形成剛度矩陣時凝聚掉。
下面進行基於自由度數目的誤差分析,由於自由度數目在網格劃分加密的過程中並不是線性增加的,在分析中我們將對自由度數目取自然對數作為衡量指標。
可以看出半彌散單元法無論是在相同網格還是在相同自由度數目下,在計算的收斂速度和精度方面都遠遠優於傳統有限單元法。
四、亂點英雄譜
觀今時結構分析手段之江湖,豪傑並起,百舸爭流。而有限單元法憑藉其成熟的理論,嚴整的招法,為各大軟體開發商所喜,因之一統江湖,一時之間風光無兩。但其嚴整的招法既有好處也帶來了缺陷,對於一些大變形問題,由於對手招式過於詭奇,應付起來難免措手不及。而無網格法招式輕靈多變,應付有限單元法所不能解決的對手時長袖善舞,應付自如,但也因自身套路難以琢磨而難以為眾人所用。廣義有限單元法則有些集二者之長的意味,但作為新晉高手,至今仍未顯大名於世。
本人學識粗鄙,天下英雄實不能盡識,如有錯漏,還望各位看官海涵。
參考文獻
[1] 王勖成.有限單元法.北京:清華大學出版社,2003
[2] 陳曉珞. 二次廣義有限元性質及其邊界條件的處理方法:[碩士學位論文]. 上海:同濟大學,2011
[3] 李剛. 半彌散單元法:[碩士學位論文]. 上海:同濟大學,2007
來源:非解構(ID:non-structure),本文已獲授權,對原作者表示感謝。