有限元法(FEM) 附有限元仿真實踐原理下載
有限元法簡介
空間和時間相關問題的物理定律通常用偏微分方程(PDE)來描述。對于絕大多數(shù)的幾何結構和所面對的問題來說,可能無法求出這些偏微分方程的解析解。不過,在通常的情況下,可以根據(jù)不同的離散化 類型來構造出近似的方程,得出與這些偏微分方程近似的數(shù)值模型方程,并可以用數(shù)值方法求解。如此,這些數(shù)值模型方程的解就是相應的偏微分方程真實解的近似解。有限元法(FEM)就是用來計算出這些近似解的。
在此,ψi 代表這些基函數(shù),而 ui 則代表用來對 u 進行近似的 uh 函數(shù)中的系數(shù)。下圖用一個一維問題闡明這一原理。例如,u 可以表示某一均勻加熱的桿在特定長度(x)處的溫度。此圖中的線性基函數(shù)的值,在各自的節(jié)點處為 1,在其他節(jié)點處為 0。在這個例子中,函數(shù) u 的定義域所在的 x-軸部分(即這根桿的長度)共有七個單元。
函數(shù) u(藍色實線)通過 uh(紅色虛線)進行逼近,后者是線性基函數(shù)的線性組合(ψi 用黑色實線表示)。線性基函數(shù)的系數(shù)由 u0 到 u7 表示。
使用有限元法的好處之一就是該方法在離散度的選擇方面提供了極大的自由(同時包括用于離散空間和離散基函數(shù)的單元的離散度選擇)。比如,在上圖中,這些單元均勻地分布在 x-軸上(雖然并不總會是這種情況)。在函數(shù) u 的一個梯度較大的區(qū)域中,也可以使用較小的單元,如下所示。
函數(shù) u(藍色實線)通過 uh(紅色虛線)進行逼近,后者是線性基函數(shù)的線性組合(ψi 用黑色實線表示)。線性基函數(shù)的系數(shù)由 u0 到 u7 表示。
有限元法的另一個優(yōu)點是該理論已經(jīng)發(fā)展得較為成熟了,原因就在于偏微分方程問題的數(shù)值表述式和弱表達式之間的密切關系(見下面的部分)。例如,當數(shù)值模型方程在計算機上求解時,該理論在誤差估計或誤差邊界 估計方面是較為有效的。
回顧有限元方法的歷史,可知該方法是在 20 世紀 40 年代初被德裔美國數(shù)學家 Richard Courant 首次提出的。雖然 Courant 報道了有限元方法在諸多問題上的應用,但幾十年之后該方法才在結構力學之外的領域獲得了普遍的應用,成就了現(xiàn)在的地位。
對輪輞進行的結構分析,圖中顯示有限元離散化、應力和變形。
代數(shù)方程、常微分方程、偏微分方程和物理定律
物理定律通常使用數(shù)學語言來表達。例如,各類守恒定律(如能量守恒定律、質量守恒定律和動量守恒定律等)都可以用偏微分方程(PDE)來表達。這些定律也可以用相關變量(包括溫度、密度、速度、電勢以及其他因變量)的本構關系來表達。
微分方程包含有相應的表達式,可以在自變量(x, y, z, t)發(fā)生變化時確定因變量的小幅變化。這一小幅變化也被稱為因變量對應于自變量的導數(shù)。
在此, 表示密度,而 Cp 則代表熱容量。溫度 T 是因變量,時間 t 是自變量。函數(shù)
可以描述隨溫度和時間而變化的一個熱源。方程 (3) 表明,如果溫度在隨著時間而變化,則它必然會由熱源
所平衡(或所引起)。此方程是用一個自變量(t)的導數(shù)所表示的一個微分方程。這種微分方程被稱為常微分方程(ODE)。
在某些情況下,當某一時間的溫度 t0 為已知時(稱為初始條件),即可得到方程 (3) 的一個解析解,表達式如下:
(4)
如此,該固體中的溫度通過一個代數(shù)方程(4)來表示,其中的某個時間值 t1 就會有一個對應時間的溫度值 T1。
因此,方程(5)表明,在所有方向上都有了改變時,如果凈通量發(fā)生了變化,以至于 q 的發(fā)散(變化的總和)不為零,則必須有一個熱源以及/或者隨時間變化的溫度變化來進行平衡(或引發(fā))。
在此,導數(shù)是以 t、x、y 和 z 表示的。在某個微分方程是用一個以上的自變量的導數(shù)來表示的情況下, 該微分方程就被稱為偏微分方程(PDE),這是因為每個導數(shù)都可能代表(幾個可能方向中的)某個方向上的變化。還需注意的是,常微分方程中的導數(shù)是用 d 來表示的,而偏微分方程中的導數(shù)則是用更卷曲的 ? 來表示的。
在不用解析法求解偏微分方程的前提下,另一種方案就是通過尋找近似的數(shù)值解 來求解數(shù)值模型方程。有限元法正是這種類型的方法——一種求解偏微分方程的數(shù)值方法。
類似于上面提到的熱能守恒方程,可以推導出動量守恒與質量守恒的方程(這兩個方程構成了流體動力學的基礎)。此外,亦可以推導出空變與時變問題中的電磁場和通量方程,從而得到偏微分方程組。
繼續(xù)這一討論,讓我們看看如何從偏微分方程中推導出所謂的弱形式公式。
源自于弱公式的有限元法:基函數(shù)和試函數(shù)
其中,h 表示傳熱系數(shù),Tamb 表示環(huán)境溫度。邊界表面上向外的單位法向矢量由 n 表示。方程(10)至(13)描述了這一散熱器的數(shù)學模型,如下圖所示。
下一步是將方程(10)的兩邊都乘以一個試函數(shù) φ,并在域 Ω 上積分:
(14)
試函數(shù) φ 與方程的解 T 被假定屬于希爾伯特空間(Hilbert space)。希爾伯特空間是一個具有無限維度的函數(shù)空間,并帶有具備特定屬性的函數(shù)。它可以被看作是具有一定屬性的函數(shù)的集合;這樣一來,這些函數(shù)可以同向量空間中的普通向量一樣被方便地操作。例如,可以在該集合中生成函數(shù)的線性組合(這些函數(shù)有明確的長度,稱為模 ),并且可以像歐幾里德矢量一樣測量函數(shù)之間的角度。
實際上,可以通過有限元方法簡單地將這些函數(shù)轉換為普通的矢量。有限元法是一種系統(tǒng)性的方法,將無限維函數(shù)空間中的函數(shù)轉換為有限維函數(shù)空間中的一類函數(shù),最后再轉換為可以用數(shù)值方法處理的普通矢量(在某一矢量空間中)。
通過要求此等式對希爾伯特空間中的所有 試函數(shù)都成立,可以實現(xiàn)方程(10)的弱形式公式或稱為變分公式。之所以說是“弱”,是因為其放寬了(10)的要求,也就是偏微分方程的各項在每一個點上都必須被明確定義的要求。相反的是,只有在積分時才要求(14)和(15)是相等的關系。例如,弱公式化完全允許解的一階導數(shù)不連續(xù),因為這種情況并不妨礙積分。但是,它為二階導數(shù)引入的分布 則并不是普通意義上的函數(shù)。因此,在不連續(xù)點上要求(10)成立是沒有意義的。
有時可以對某個分布進行積分,以使(14)被明確定義。可以證明的是,弱公式化以及通過(13)得到的邊界條件(11)都是與通過逐點公式化求出的解直接相關的。此外,對于解可微分 的情況(即二階導數(shù)明確定義),這些解是相同的。這些公式化是等效的,因為從(10)推導(15)的過程依賴于格林第一恒等式,而其只有在 T 有連續(xù)的二階導數(shù)的情況下才成立。
這是有限元公式化的第一步。利用弱公式化,就有可能對數(shù)學模型方程進行離散化,從而得到數(shù)值模型方程。可以利用伽遼金法——許多可能的有限元法公式化中的一種——來進行離散化。
這里的未知數(shù),就是函數(shù) T(x) 的近似解中的系數(shù) Ti。隨后,方程(17)就變成了一個方程組,該方程組與有限維函數(shù)空間擁有相同的維度。如果使用了試函數(shù) ψj 中的數(shù)字 n,使 j 從 1 一直變到 n,那么就可以根據(jù)(17)得到一個方程數(shù)量為 n 的方程組。方程(16)中也有 n 個未知的系數(shù)(Ti)。

來自之前的散熱器模型圖的有限元離散化。

其中,T 是未知矢量,且 T h = {T1, .., Ti, …, Tn};A 則是一個 nxn 的矩陣,其元素 Aji 中的每個方程 j 都含有系數(shù) Ti。右邊是維度從 1 到 n 的矢量。A 是系統(tǒng)矩陣,通常稱為(消除)的剛度矩陣 ——這是有限元方法的首次應用,也是其在結構力學中的用途。
如果源函數(shù)在溫度方面是非線性的,或者傳熱系數(shù)取決于溫度,那么該方程組也是非線性的,矢量 b 就成為了未知系數(shù) Ti 的一個非線性函數(shù)。
有限元方法的優(yōu)點之一是它能夠選擇試函數(shù)和基函數(shù)。在非常小的幾何區(qū)域的支集之上,是有可能選擇試函數(shù)和基函數(shù)的。這意味著,方程(17)在任意一處都為零——除非是在函數(shù) ψj 和 ψi 重疊的非常有限的區(qū)域上,因為上面所有的積分都包括了函數(shù) i 和 j(或它們的梯度)的乘積。很難用三維空間來描述試函數(shù)和基函數(shù)的支集,但其二維的類比卻是能夠被可視化的。
假設有一個二維的幾何域,并且選用了 x 和 y 的線性函數(shù),每個函數(shù)在點 i 上的值為 1,但在其他點 k 上的值為零。下一步是使用三角形對這一二維域進行離散化,并為某一三角形網(wǎng)格中的兩個相鄰節(jié)點 i 和 j 給出基函數(shù)(試函數(shù)或形函數(shù))。

兩個相鄰的基函數(shù)共享兩個三角形的單元。因此,兩個基函數(shù)之間有一些重疊,如上所示。此外,請注意,如果 i = j,則函數(shù)之間會完全重疊。這些貢獻形成了未知矢量 T 的系數(shù),這一未知矢量與系統(tǒng)矩陣的對角線分量 Ajj 相對應。
比如說,假設現(xiàn)在這兩個基函數(shù)更進一步地分開了。這兩個函數(shù)不共享單元,但它們有一個共同的單元頂點。如下圖所示,它們不重疊。
共享一個單元頂點的兩個基函數(shù)在二維域中不發(fā)生重疊。
當這兩個基函數(shù)重疊時,方程(17)具有非零值,且對系統(tǒng)矩陣的貢獻也是非零的。當沒有重疊時,積分為零,因此對系統(tǒng)矩陣的貢獻也為零。
這意味著,在從 1 到 n 的節(jié)點上,對(17)的方程組中的每個方程來說,它們都只能從共享同一個單元的相鄰節(jié)點中得到若干個非零的項。方程(18)中的系統(tǒng)矩陣 A 變得稀疏,而對應于重疊 ij:s 的矩陣分量才有非零項。這一代數(shù)方程組的解可以作為該偏微分方程的近似解。網(wǎng)格越稠密,近似解就越接近真實解。
對散熱器中的溫度場進行的有限元近似。
瞬態(tài)問題(時變問題)
在此,系數(shù) Ti 是時變函數(shù),而基函數(shù)和試函數(shù)則僅依賴于空間坐標。再者,在時間域上的時間導數(shù)不是離散的。
一種方法是對時間域也使用有限元法,但這種做法可能會耗費大量的計算資源。經(jīng)常采取的另一種方案則是通過直線法來對時間域進行獨立的離散化。比如可以使用有限差分法。其最簡單的形式可以用下面的差分近似法來表示:
(20)
在面對線性問題時,在每一個時間步長上都需要求解一個線性方程組。如果是非線性的問題,則必須在每個時間步長內求解相應的非線性方程組。由于在 t + Δt 處的解是被方程(21)隱含地給出的,所以這種時間推進方案被稱為隱式法。
該式表明,一旦在某一給定時間上的解(Ti,t)已知,那么方程(22)就能顯式地給出在 t + Δt (Ti, t+Δt) 處的解。換言之,對于一個顯式的時間推進方案,不需要在每個時間步長上都求解一個方程組。顯式時間推進方案的缺點是它們有一個穩(wěn)定性方面的時間步進限制。對于熱問題來說(如此處所強調的情況),顯式方法需要非常短的時間步長。隱式方案允許更大的時間步長,減少了如(22)這樣的方程所需的計算資源(在每一個時間步長上都要對這些方程進行求解)。
在實踐中,現(xiàn)代化的時間步進算法會根據(jù)具體問題自動在顯式和隱式步進法之間切換。此外,方程(20)中的差分方程被替換為一個多項式,其階次和步長可以發(fā)生變化,具體取決于所要解決的問題和求解所需的時間。現(xiàn)代化的時間推進方案會根據(jù)數(shù)值解的時間演化來自動地控制多項式的階次以及步長。
下面有幾個例子,對最常用的幾種方法加以說明:
向后微分公式(BDF)法
廣義 α 法
不同的 Runge-Kutta 法
不同的單元
如上所述,伽遼金法采用了與基函數(shù)和試函數(shù)相同的函數(shù)集。然而,即使是這種方法,也可以通過很多種方式(理論上是無窮多的)來定義基函數(shù)(即伽遼金有限元公式中的單元)。讓我們來回顧一下最常用的幾種單元。
對于二維和三維的線性函數(shù),最常見的元素如下圖所示。此圖和上圖 給出的是線性基函數(shù)(被定義在三角形網(wǎng)格中,形成了三角形的線性單元)。基函數(shù)被表示為節(jié)點位置(二維時:x 和 y;三維時:x、y 和 z)的函數(shù)。
在二維面上,矩形單元常常被用于結構力學分析。它們還可用于計算流體動力學(CFD)和傳熱建模中的邊界層網(wǎng)格剖分。它們的三維類比就是所謂的六面體單元,后者也常被應用于結構力學和邊界層網(wǎng)格剖分。在從六面體邊界層單元到四面體單元的過渡中,錐體單元通常被放置在邊界層單元的頂端。
二維和三維線性單元的節(jié)點位置與幾何形狀。
下圖顯示的是相應的二階單元(二次單元)。在此,面對一個域邊界的邊和面通常是彎曲的,而面對該域內部的邊和面則是直線或平面。但是請注意,也可以將所有的邊和曲都定義為是彎曲的。拉格朗日單元和巧湊邊點元是二維和三維建模中最常用的單元類型。拉格朗日單元使用下面所有的節(jié)點(黑色、白色和灰色),而巧湊邊點元則不使用灰色的節(jié)點。
二階單元。如果移除灰色節(jié)點,便可得到相應的巧湊邊點單元。黑色、白色和灰色節(jié)點都存在于拉格朗日單元中。
博客“在多物理場模型中追蹤單元階次”中給出了二階(二次)拉格朗日元的二維圖形,非常漂亮。在上述單元的內部,很難用三維的形式描述這些二次基函數(shù)的基,但是可以用色塊來表示單元表面的函數(shù)數(shù)值。
在討論有限元法時,需要考慮的一個重要因素就是誤差估計。原因在于,當達到估計出的誤差寬容度時,就會發(fā)生收斂。請注意,這里的討論具有更一般的性質,而不是局限于特定的有限元方法。
有限元法給出的是數(shù)學模型方程的一個近似解。數(shù)值方程的解與數(shù)學模型方程的精確解之間的差值就是誤差:e = u - uh。
在許多情況下,可以在得出數(shù)值方程的解之前就估計出誤差的大小(即先驗 誤差估計)。先驗 估計通常僅用于預測所用有限元方法的收斂階數(shù)。例如,如果某個問題是適定的,并且相應的數(shù)值方法可以收斂,那么根據(jù) O(hα)(其中 α 表示收斂階數(shù)),隨著通常的單元尺寸 h 的減小,誤差的模也會減小。由此可見,隨著網(wǎng)格密度的增加,誤差的模也會快速地降低。
不過,只有簡單的問題才能進行先驗 估計。此外,估計出的結果往往會包含不同的未知常數(shù),從而不可能給出定量的預測。后驗估計使用的則是近似解,并結合了相關問題的其他近似,以估計出誤差的模。
構造解方法
一個非常簡單但卻通用的誤差估計方法(用于數(shù)值方法和偏微分問題),就是對問題進行略微改動——如這一篇博客文章 所述—— 使預定義的解析表達式成為改動后的問題的真實解。這種方法的優(yōu)點是未對數(shù)值方法或其背后的數(shù)學問題進行過假設。此外,由于解是已知的,所以可以很容易地計算出誤差的大小。通過謹慎地選擇分析表達式,就可以對方法和問題的不同方面進行研究。
如此,就可以為不同選擇的離散化程度和 計算出誤差及其模。如果改動后問題的解與未改動問題的解具有相同的特性,那么改動后問題的誤差就可以用作未改動問題的近似誤差。在實踐中,可能很難知道是否是這種情況——這是此方法的缺點。這種方法的優(yōu)點在于它的簡單性和普遍性:既可以用于非線性問題和時變問題(瞬態(tài)問題),也可以用于任何的數(shù)值方法。
目標定向的誤差估計
如果可以從近似解中選擇出一個函數(shù)(或泛函數(shù)),并將其作為一個重點物理量來進行誤差估計,那么就可以通過解析方法精準地估算出此物理量的計算誤差(或界限)。此類估計依賴于對偏微分方程殘差的后驗 計算,以及對所謂的對偶問題 進行的近似求解。對偶問題與所選擇的函數(shù)是直接相關的(并由其定義)。
這種方法的缺點在于其依賴于“對偶問題”的精確計算,而且只給出了所選函數(shù)的誤差估計(而沒有涉及其他物理量)。這種方法的優(yōu)勢在于其較高的普適性和較合理的資源消耗(用于誤差計算)。
網(wǎng)格收斂
在實踐中,要對非常精細的網(wǎng)格剖分方案(比實際所需精細得多的方案)進行近似求解,其實是較困難的。因此,在習慣上會使用最精細網(wǎng)格的近似來達到此目的。對每個網(wǎng)格細化來說,也可以從所得解的變化情況中估算出收斂性。如果近似解位于一個收斂的區(qū)域之中,那么這個解的變化幅度會隨著網(wǎng)格的細化而收窄;如此,所得的近似解也就會越來越接近于真實解。
下圖顯示了一個橢圓形膜的結構力學基準模型;得益于對稱性,只需要對該膜的四分之一部分進行計算就可以了。載荷作用在該幾何體的外邊緣上,而沿 x 和 y 軸的邊界被認為是對稱的。
橢圓薄膜的基準模型,其中假設沿 x 和 y 軸(滾動支座)的邊呈對稱分布,并在外部邊上施加載荷。
對不同網(wǎng)格類型和單元尺寸的數(shù)值模型方程進行求解。例如,下圖描述了用于二次基函數(shù)的矩形拉格朗日單元,這些基函數(shù)是根據(jù)上圖得出的。
根據(jù)更早給出的這幅圖,對該點上的應力和應變進行了計算。下面的圖表顯示的是此點上的 σx 所得的相對值。此值應為零,因此與零值的任何差異都是一種誤差。為了得到一個相對誤差,將計算出的 σx 除以計算出的 σy,以便為相對誤差的估計給出正確的數(shù)量級。
顯示不同單元和單元尺寸(單元尺寸 = h)的前圖中計算點處的 σx 的相對誤差。四邊形是指矩形單元,它們可以是線性的,也可以具有二次基函數(shù)。
上圖表明,隨著各單元的單元尺寸(h)的減小,相對誤差也相應減小。在這種情況下,隨著基函數(shù)階次(單元階次)的升高,收斂曲線也變得更為陡峭。不過,需要注意的是,在單元尺寸一定的情況下,隨著階次的升高,數(shù)值模型中的未知項的數(shù)量也會增加。這就意味著,當我們增加單元的階次時,我們也要為更高的精確度而付出代價,這種代價就是計算耗時的增加。如果不使用更高階的單元,還可以采取的另一種方法就是為較低階次的單元選用更細化的網(wǎng)格。
網(wǎng)格自適應
在計算出了這些數(shù)值方程的解 uh 之后,就可以用后驗 局部誤差估計值來創(chuàng)建一個密度更大的網(wǎng)格,該網(wǎng)格具有較大的誤差。然后可以使用細化的網(wǎng)格來計算出第二個近似解。
下圖描述了一個被加熱的圓柱體在受到穩(wěn)態(tài)流體流動作用下的溫度場。對這一穩(wěn)態(tài)問題進行了兩次求解:一次是用基礎網(wǎng)格,另一次是用一個細化網(wǎng)格(被基本網(wǎng)格計算出的誤差估計所控制)。該細化網(wǎng)格在溫度和熱通量方面的計算精度更高,而這一點可能正是該實例所需要的。
在流動作用下的受熱圓柱體周圍的溫度場計算結果,上圖未經(jīng)網(wǎng)格細化,下圖經(jīng)過了網(wǎng)格細化。
對時變(瞬態(tài))的對流問題來說,也可以通過前序時步的解來實現(xiàn)對流網(wǎng)格的細化。在下圖給出的例子中,相場被用來計算噴墨打印機中的墨水液滴與空氣之間的界面。該界面是由相場函數(shù)的等值面所給出的,其值等于 0.5。在這個界面上,相場函數(shù)的值迅速地從 1 變?yōu)?0。在此相場函數(shù)的這些陡峭梯度的周圍,我們可以使用誤差估計來自動完成網(wǎng)格細化的工作,而流場則可以用來對流網(wǎng)格細化,以便僅在相場等值面的面前才使用更細的網(wǎng)格。
在一個瞬態(tài)兩相流問題中,對噴墨打印機中的一串墨滴進行網(wǎng)格細化。
其他有限元公式
在上述例子中,我們?yōu)榛瘮?shù)和試函數(shù)使用了相同的函數(shù)集來實現(xiàn)模型方程的離散化。如果一個有限元公式可以使試函數(shù)不同于基函數(shù),則該公式稱為 Petrov-Galerkin 法。這是一種常用的方法;例如,在解決對流-擴散問題的過程中,只會對流線方向進行穩(wěn)定化處理。其也被稱為流線迎風 /Petrov-Galerkin(SUPG)法。
在耦合方程組的求解過程中,不同的因變量可能會用到不同的基函數(shù)。一個典型的例子是納維-斯托克斯方程的求解,其中的壓力往往比速度更平滑、更易進行近似。在某類方法中,如果一個耦合方程組中不同的因變量的基函數(shù)(以及試函數(shù))屬于不同的函數(shù)空間,那么這類方法便稱為混合有限元法。
COMSOL Multiphysics 軟件中用于流體流動分析的混合單元法的設置,其中二次形函數(shù)(基函數(shù))用于計算速度,線性形函數(shù)用于計算壓力。
下載地址:有限元仿真實踐原理
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP




















