理論基礎:有限元法的較為詳細的說明

2022-03-27  by:CAE仿真在線  來源:互聯網

有限元法簡介

空間和時間相關問題的物理定律通常用偏微分方程(PDE)來描述。對于絕大多數的幾何結構和所面對的問題來說,可能無法求出這些偏微分方程的解析解。不過,在通常的情況下,可以根據不同的離散化 類型來構造出近似的方程,得出與這些偏微分方程近似的數值模型方程,并可以用數值方法求解。如此,這些數值模型方程的解就是相應的偏微分方程真實解的近似解。有限元法(FEM)就是用來計算出這些近似解的。

舉例來說,某函數 u 可能是一個偏微分方程中的因變量(即溫度、電勢、壓力等)?？梢愿鶕铝斜磉_式,通過基函數的線性組合將函數 u 近似為新的函數 uh:

(1)

以及

(2)

在此,ψi 代表這些基函數,而 ui 則代表用來對 u 進行近似的 uh 函數中的系數。下圖用一個一維問題闡明這一原理。例如,u 可以表示某一均勻加熱的桿在特定長度(x)處的溫度。此圖中的線性基函數的值,在各自的節點處為 1,在其他節點處為 0。在這個例子中,函數 u 的定義域所在的 x-軸部分(即這根桿的長度)共有七個單元。

使用有限元法的好處之一就是該方法在離散度的選擇方面提供了極大的自由(同時包括用于離散空間和離散基函數的單元的離散度選擇)。比如,在上圖中,這些單元均勻地分布在 x-軸上(雖然并不總會是這種情況)。在函數 u 的一個梯度較大的區域中,也可以使用較小的單元,如下所示。

這兩幅圖都表明,選定的線性基函數在 x-軸方向上獲得的支持(僅有一個較為狹窄的非零區間)和重疊非常有限。根據手頭的問題,可以選擇其他類型的函數而不是線性函數。

有限元法的另一個優點是該理論已經發展得較為成熟了,原因就在于偏微分方程問題的數值表述式和弱表達式之間的密切關系(見下面的部分)。例如,當數值模型方程在計算機上求解時,該理論在誤差估計或誤差邊界 估計方面是較為有效的。

回顧有限元方法的歷史,可知該方法是在 20 世紀 40 年代初被德裔美國數學家 Richard Courant 首次提出的。雖然 Courant 報道了有限元方法在諸多問題上的應用,但幾十年之后該方法才在結構力學之外的領域獲得了普遍的應用,成就了現在的地位。

代數方程、常微分方程、偏微分方程和物理定律

物理定律通常使用數學語言來表達。例如,各類守恒定律(如能量守恒定律、質量守恒定律和動量守恒定律等)都可以用偏微分方程(PDE)來表達。這些定律也可以用相關變量(包括溫度、密度、速度、電勢以及其他因變量)的本構關系來表達。

微分方程包含有相應的表達式,可以在自變量(x, y, z, t)發生變化時確定因變量的小幅變化。這一小幅變化也被稱為因變量對應于自變量的導數。

假設一個固體具有時變溫度,但在空間上的變化忽略不計。在這種情況下,通過內能(熱)守恒方程,就可以導出在熱源 g 的作用下,隨著時間的小幅變化而發生的溫度變化的方程式:

(3)

在此, 表示密度,而 Cp 則代表熱容量。溫度 T 是因變量,時間 t 是自變量。函數  可以描述隨溫度和時間而變化的一個熱源。方程 (3) 表明,如果溫度在隨著時間而變化,則它必然會由熱源  所平衡(或所引起)。此方程是用一個自變量(t)的導數所表示的一個微分方程。這種微分方程被稱為常微分方程(ODE)。

在某些情況下,當某一時間的溫度 t0 為已知時(稱為初始條件),即可得到方程 (3) 的一個解析解,表達式如下:

(4)

如此,該固體中的溫度通過一個代數方程(4)來表示,其中的某個時間值 t1 就會有一個對應時間的溫度值 T1。

物理屬性常常會隨著時間和空間發生變化。例如,該固體中靠近熱源處的溫度可能比其他位置略高。這種溫度差異進而引起該固體內部不同部分之間產生熱通量。在這種情況下,根據能量守恒定律就可以導出一個傳熱方程,該方程同時具有時間變量和空間變量(x),如:

(5)

同之前一樣,T 是因變量,而 x(x = (x, y, z))和 t 則是自變量。該固體中的熱通量矢量由 q = (qx, y, qz) 表示,而 q 的發散 則描述了熱通量沿著空間坐標的變化。在笛卡爾坐標系中,q 的發散被定義為:

(6)

因此,方程(5)表明,在所有方向上都有了改變時,如果凈通量發生了變化,以至于 q 的發散(變化的總和)不為零,則必須有一個熱源以及/或者隨時間變化的溫度變化來進行平衡(或引發)。

可以通過傳導熱通量的本構關系來描述一個固體中的熱通量,這也稱為傅里葉定律:

(7)

在上述方程式中,k 表示導熱系數。方程(7)表明,在導熱系數為比例常數的情況下,熱通量與溫度梯度成正比。方程(7)((5) 中)給出了以下的微分方程:

(8)

在此,導數是以 t、x、y 和 z 表示的。在某個微分方程是用一個以上的自變量的導數來表示的情況下, 該微分方程就被稱為偏微分方程(PDE),這是因為每個導數都可能代表(幾個可能方向中的)某個方向上的變化。還需注意的是,常微分方程中的導數是用 d 來表示的,而偏微分方程中的導數則是用更卷曲的 ? 來表示的。

除了方程(8),還可以知道的就是某個時間 t0 上的溫度或者某個位置 x0 上的熱通量。此類知識可應用于方程(8)的初始條件和邊界條件。在許多情況下,偏微分方程都無法通過解析方法來求解(即得出不同時間和位置下的因變量的值)。有時,要得到一個如下的解析表達式,可能非常困難,甚至幾乎是不可能的,例如方程(8)中的:

(9)

在不用解析法求解偏微分方程的前提下,另一種方案就是通過尋找近似的數值解 來求解數值模型方程。有限元法正是這種類型的方法——一種求解偏微分方程的數值方法。

類似于上面提到的熱能守恒方程,可以推導出動量守恒與質量守恒的方程(這兩個方程構成了流體動力學的基礎)。此外,亦可以推導出空變與時變問題中的電磁場和通量方程,從而得到偏微分方程組。

繼續這一討論,讓我們看看如何從偏微分方程中推導出所謂的弱形式公式。

源自于弱公式的有限元法:基函數和試函數

假定正在研究的一個散熱器中的溫度分布由方程(8)給出,但現正處于穩定狀態,這就意味著方程(8)中的溫度場的時間導數為零。模型域 Ω 的域方程如下:

(10)

此外,假定沿邊界(?Ω1)的溫度已知,同時垂直于其他一些邊界(?Ω2)的熱通量的表達式也已知。在其余的邊界上,熱通量在向外的方向(?Ω3)上為零。這些邊界上的邊界條件就成為:

(11)

(12)

(13)

其中,h 表示傳熱系數,Tamb 表示環境溫度。邊界表面上向外的單位法向矢量由 n 表示。 方程(10)至(13)描述了這一散熱器的數學模型,如下圖所示。

下一步是將方程(10)的兩邊都乘以一個試函數 φ,并在域 Ω 上積分:

(14)

試函數 φ 與方程的解 T 被假定屬于希爾伯特空間(Hilbert space)。希爾伯特空間是一個具有無限維度的函數空間,并帶有具備特定屬性的函數。它可以被看作是具有一定屬性的函數的集合;這樣一來,這些函數可以同向量空間中的普通向量一樣被方便地操作。例如,可以在該集合中生成函數的線性組合(這些函數有明確的長度,稱為模 ),并且可以像歐幾里德矢量一樣測量函數之間的角度。

實際上,可以通過有限元方法簡單地將這些函數轉換為普通的矢量。有限元法是一種系統性的方法,將無限維函數空間中的函數轉換為有限維函數空間中的一類函數,最后再轉換為可以用數值方法處理的普通矢量(在某一矢量空間中)。

如果要求(14)對試函數空間中的所有試函數都成立,而不是方程(10)對 Ω 中的所有點都成立,則可以得到弱形式公式。因此,基于方程(10)的問題公式有時也稱為逐點公式。在我們所說的伽遼金法 中,假設解 T 同測試函數屬于相同的希爾伯特空間。這通常寫為 φ ? h 和 T ? h,其中 H 表示希爾伯特空間。 使用格林第一恒等式(實質上是進行分部積分), 就可以推出以下方程(14):

(15)

通過要求此等式對希爾伯特空間中的所有 試函數都成立,可以實現方程(10)的弱形式公式或稱為變分公式。之所以說是“弱”,是因為其放寬了(10)的要求,也就是偏微分方程的各項在每一個點上都必須被明確定義的要求。相反的是,只有在積分時才要求(14)和(15)是相等的關系。例如,弱公式化完全允許解的一階導數不連續,因為這種情況并不妨礙積分。但是,它為二階導數引入的分布 則并不是普通意義上的函數。因此,在不連續點上要求(10)成立是沒有意義的。

有時可以對某個分布進行積分,以使(14)被明確定義?？梢宰C明的是,弱公式化以及通過(13)得到的邊界條件(11)都是與通過逐點公式化求出的解直接相關的。此外,對于解可微分 的情況(即二階導數明確定義),這些解是相同的。這些公式化是等效的,因為從(10)推導(15)的過程依賴于格林第一恒等式,而其只有在 T 有連續的二階導數的情況下才成立。

這是有限元公式化的第一步。利用弱公式化,就有可能對數學模型方程進行離散化,從而得到數值模型方程?？梢岳觅み|金法——許多可能的有限元法公式化中的一種——來進行離散化。

首先,要實現離散化,就意味著要在希爾伯特空間 H 的有限維子空間中尋找方程(15)的近似解;如此,T ≈ Th。這就是說,近似解被表示為一組屬于子空間的基函數 ψi 的線性組合:

(16)

由此,對每個試函數 ψj 而言,方程(15)的離散化形式即變為:

(17)

這里的未知數,就是函數 T(x) 的近似解中的系數 Ti。隨后,方程(17)就變成了一個方程組,該方程組與有限維函數空間擁有相同的維度。如果使用了試函數 ψj 中的數字 n,使 j 從 1 一直變到 n,那么就可以根據(17)得到一個方程數量為 n 的方程組。方程(16)中也有 n 個未知的系數(Ti)。

一旦體系被離散化并被施加了邊界條件后, 根據以下表達式就可以得到一個方程組:

(18)

其中,T 是未知矢量,且 T h = {T1, .., Ti, …, Tn};A 則是一個 nxn 的矩陣,其元素 Aji 中的每個方程 j 都含有系數 Ti。右邊是維度從 1 到 n 的矢量。A 是系統矩陣,通常稱為(消除)的剛度矩陣 ——這是有限元方法的首次應用,也是其在結構力學中的用途。

如果源函數在溫度方面是非線性的,或者傳熱系數取決于溫度,那么該方程組也是非線性的,矢量 b 就成為了未知系數 Ti 的一個非線性函數。

有限元方法的優點之一是它能夠選擇試函數和基函數。在非常小的幾何區域的支集之上,是有可能選擇試函數和基函數的。這意味著,方程(17)在任意一處都為零——除非是在函數 ψj 和 ψi 重疊的非常有限的區域上,因為上面所有的積分都包括了函數 i 和 j(或它們的梯度)的乘積。很難用三維空間來描述試函數和基函數的支集,但其二維的類比卻是能夠被可視化的。

假設有一個二維的幾何域,并且選用了 x 和 y 的線性函數,每個函數在點 i 上的值為 1,但在其他點 k 上的值為零。下一步是使用三角形對這一二維域進行離散化,并為某一三角形網格中的兩個相鄰節點 i 和 j 給出基函數(試函數或形函數)。

兩個相鄰的基函數共享兩個三角形的單元。因此,兩個基函數之間有一些重疊,如上所示。此外,請注意,如果 i = j,則函數之間會完全重疊。這些貢獻形成了未知矢量 T 的系數,這一未知矢量與系統矩陣的對角線分量 Ajj 相對應。

比如說,假設現在這兩個基函數更進一步地分開了。這兩個函數不共享單元,但它們有一個共同的單元頂點。如下圖所示,它們不重疊。

當這兩個基函數重疊時,方程(17)具有非零值,且對系統矩陣的貢獻也是非零的。當沒有重疊時,積分為零,因此對系統矩陣的貢獻也為零。

這意味著,在從 1 到 n 的節點上,對(17)的方程組中的每個方程來說,它們都只能從共享同一個單元的相鄰節點中得到若干個非零的項。方程(18)中的系統矩陣 A 變得稀疏,而對應于重疊 ij:s 的矩陣分量才有非零項。這一代數方程組的解可以作為該偏微分方程的近似解。網格越稠密,近似解就越接近真實解。

瞬態問題(時變問題)

可以在瞬態(時變)的情況下進一步定義該散熱器中的熱能平衡。根據伽遼金方法,每個試函數 ψj 的離散弱公式化可以寫作:

(19)

在此,系數 Ti 是時變函數,而基函數和試函數則僅依賴于空間坐標。再者,在時間域上的時間導數不是離散的。

一種方法是對時間域也使用有限元法,但這種做法可能會耗費大量的計算資源。經常采取的另一種方案則是通過直線法來對時間域進行獨立的離散化。比如可以使用有限差分法。其最簡單的形式可以用下面的差分近似法來表示:

(20)

給出的是方程(19)中的兩個可能有限差分逼近。當未知的系數 Ti,t 以 t + Δt 的形式表示時,就可以得到第一個式子:

(21)

在面對線性問題時,在每一個時間步長上都需要求解一個線性方程組。如果是非線性的問題,則必須在每個時間步長內求解相應的非線性方程組。由于在 t + Δt 處的解是被方程(21)隱含地給出的,所以這種時間推進方案被稱為隱式法。

第二個公式則基于 t 處的解:

(22)

該式表明,一旦在某一給定時間上的解(Ti,t)已知,那么方程(22)就能顯式地給出在 t + Δt (Ti, t+Δt) 處的解。換言之,對于一個顯式的時間推進方案,不需要在每個時間步長上都求解一個方程組。顯式時間推進方案的缺點是它們有一個穩定性方面的時間步進限制。對于熱問題來說(如此處所強調的情況),顯式方法需要非常短的時間步長。隱式方案允許更大的時間步長,減少了如(22)這樣的方程所需的計算資源(在每一個時間步長上都要對這些方程進行求解)。

在實踐中,現代化的時間步進算法會根據具體問題自動在顯式和隱式步進法之間切換。此外,方程(20)中的差分方程被替換為一個多項式,其階次和步長可以發生變化,具體取決于所要解決的問題和求解所需的時間?，F代化的時間推進方案會根據數值解的時間演化來自動地控制多項式的階次以及步長。

下面有幾個例子,對最常用的幾種方法加以說明:

• 向后微分公式(BDF)法
• 廣義 α 法
• 不同的 Runge-Kutta 法

不同的單元

如上所述,伽遼金法采用了與基函數和試函數相同的函數集。然而,即使是這種方法,也可以通過很多種方式(理論上是無窮多的)來定義基函數(即伽遼金有限元公式中的單元)。讓我們來回顧一下最常用的幾種單元。

對于二維和三維的線性函數,最常見的元素如下圖所示。此圖和上圖 給出的是線性基函數(被定義在三角形網格中,形成了三角形的線性單元)?；瘮当槐硎緸楣濣c位置(二維時:x 和 y;三維時:x、y 和 z)的函數。

在二維面上,矩形單元常常被用于結構力學分析。它們還可用于計算流體動力學(CFD)和傳熱建模中的邊界層網格剖分。它們的三維類比就是所謂的六面體單元,后者也常被應用于結構力學和邊界層網格剖分。在從六面體邊界層單元到四面體單元的過渡中,錐體單元通常被放置在邊界層單元的頂端。

下圖顯示的是相應的二階單元(二次單元)。在此,面對一個域邊界的邊和面通常是彎曲的,而面對該域內部的邊和面則是直線或平面。但是請注意,也可以將所有的邊和曲都定義為是彎曲的。拉格朗日單元和巧湊邊點元是二維和三維建模中最常用的單元類型。拉格朗日單元使用下面所有的節點(黑色、白色和灰色),而巧湊邊點元則不使用灰色的節點。

博客“在多物理場模型中追蹤單元階次”中給出了二階(二次)拉格朗日元的二維圖形,非常漂亮。在上述單元的內部,很難用三維的形式描述這些二次基函數的基,但是可以用色塊來表示單元表面的函數數值。

在討論有限元法時,需要考慮的一個重要因素就是誤差估計。原因在于,當達到估計出的誤差寬容度時,就會發生收斂。請注意,這里的討論具有更一般的性質,而不是局限于特定的有限元方法。

有限元法給出的是數學模型方程的一個近似解。數值方程的解與數學模型方程的精確解之間的差值就是誤差:e = u - uh。

在許多情況下,可以在得出數值方程的解之前就估計出誤差的大小(即先驗 誤差估計)。先驗 估計通常僅用于預測所用有限元方法的收斂階數。例如,如果某個問題是適定的,并且相應的數值方法可以收斂,那么根據 O(hα)(其中 α 表示收斂階數),隨著通常的單元尺寸 h 的減小,誤差的模也會減小。由此可見,隨著網格密度的增加,誤差的模也會快速地降低。

不過,只有簡單的問題才能進行先驗 估計。此外,估計出的結果往往會包含不同的未知常數,從而不可能給出定量的預測。后驗 估計使用的則是近似解,并結合了相關問題的其他近似,以估計出誤差的模。

構造解方法

一個非常簡單但卻通用的誤差估計方法(用于數值方法和偏微分問題),就是對問題進行略微改動——如這一篇博客文章 所述—— 使預定義的解析表達式成為改動后的問題的真實解。這種方法的優點是未對數值方法或其背后的數學問題進行過假設。此外,由于解是已知的,所以可以很容易地計算出誤差的大小。通過謹慎地選擇分析表達式,就可以對方法和問題的不同方面進行研究。

讓我們來看一個例子,對這一點進行說明。假設有一種數值方法可以對一個單位正方形(Ω)上的泊松方程進行求解,且該正方形具有齊次邊界條件

(23)

(24)

此方法可用于對改動后的問題進行求解

(25)

(26)

其中,

(27)

這里,

是可以被自由選擇的一個解析表達式。另外,如果

(28)

則  是改動后的問題的精確解,此時可以直接計算出誤差大小:

(29)

如此,就可以為不同選擇的離散化程度和  計算出誤差及其模。如果改動后問題的解與未改動問題的解具有相同的特性,那么改動后問題的誤差就可以用作未改動問題的近似誤差。在實踐中,可能很難知道是否是這種情況——這是此方法的缺點。這種方法的優點在于它的簡單性和普遍性:既可以用于非線性問題和時變問題(瞬態問題),也可以用于任何的數值方法。

目標定向的誤差估計

如果可以從近似解中選擇出一個函數(或泛函數),并將其作為一個重點物理量來進行誤差估計,那么就可以通過解析方法精準地估算出此物理量的計算誤差(或界限)。此類估計依賴于對偏微分方程殘差的后驗 計算,以及對所謂的對偶問題 進行的近似求解。對偶問題與所選擇的函數是直接相關的(并由其定義)。

這種方法的缺點在于其依賴于“對偶問題”的精確計算,而且只給出了所選函數的誤差估計(而沒有涉及其他物理量)。這種方法的優勢在于其較高的普適性和較合理的資源消耗(用于誤差計算)。

網格收斂

網格收斂是一種簡單的方法,該方法比較了不同的網格剖分方案所得到的近似解。在理想情況下,一套非常精細的網格剖分方案所得出的近似解就可以作為真實解的近似了。較粗化的網格剖分方案所得出的近似解的誤差,可以由下式直接估算出:

(30)

在實踐中,要對非常精細的網格剖分方案(比實際所需精細得多的方案)進行近似求解,其實是較困難的。因此,在習慣上會使用最精細網格的近似來達到此目的。對每個網格細化來說,也可以從所得解的變化情況中估算出收斂性。如果近似解位于一個收斂的區域之中,那么這個解的變化幅度會隨著網格的細化而收窄;如此,所得的近似解也就會越來越接近于真實解。

下圖顯示了一個橢圓形膜的結構力學基準模型;得益于對稱性,只需要對該膜的四分之一部分進行計算就可以了。載荷作用在該幾何體的外邊緣上,而沿 x 和 y 軸的邊界被認為是對稱的。

對不同網格類型和單元尺寸的數值模型方程進行求解。例如,下圖描述了用于二次基函數的矩形拉格朗日單元,這些基函數是根據上圖得出的。

根據更早給出的這幅圖,對該點上的應力和應變進行了計算。下面的圖表顯示的是此點上的 σx 所得的相對值。此值應為零,因此與零值的任何差異都是一種誤差。為了得到一個相對誤差,將計算出的 σx 除以計算出的 σy,以便為相對誤差的估計給出正確的數量級。

上圖表明,隨著各單元的單元尺寸(h)的減小,相對誤差也相應減小。在這種情況下,隨著基函數階次(單元階次)的升高,收斂曲線也變得更為陡峭。不過,需要注意的是,在單元尺寸一定的情況下,隨著階次的升高,數值模型中的未知項的數量也會增加。這就意味著,當我們增加單元的階次時,我們也要為更高的精確度而付出代價,這種代價就是計算耗時的增加。如果不使用更高階的單元,還可以采取的另一種方法就是為較低階次的單元選用更細化的網格。

網格自適應

在計算出了這些數值方程的解 uh 之后,就可以用后驗 局部誤差估計值來創建一個密度更大的網格,該網格具有較大的誤差。然后可以使用細化的網格來計算出第二個近似解。

下圖描述了一個被加熱的圓柱體在受到穩態流體流動作用下的溫度場。對這一穩態問題進行了兩次求解:一次是用基礎網格,另一次是用一個細化網格(被基本網格計算出的誤差估計所控制)。該細化網格在溫度和熱通量方面的計算精度更高,而這一點可能正是該實例所需要的。

對時變(瞬態)的對流問題來說,也可以通過前序時步的解來實現對流網格的細化。在下圖給出的例子中,相場被用來計算噴墨打印機中的墨水液滴與空氣之間的界面。該界面是由相場函數的等值面所給出的,其值等于 0.5。在這個界面上,相場函數的值迅速地從 1 變為 0。在此相場函數的這些陡峭梯度的周圍,我們可以使用誤差估計來自動完成網格細化的工作,而流場則可以用來對流網格細化,以便僅在相場等值面的面前才使用更細的網格。

其他有限元公式

在上述例子中,我們為基函數和試函數使用了相同的函數集來實現模型方程的離散化。如果一個有限元公式可以使試函數不同于基函數,則該公式稱為 Petrov-Galerkin 法。這是一種常用的方法;例如,在解決對流-擴散問題的過程中,只會對流線方向進行穩定化處理。其也被稱為流線迎風 /Petrov-Galerkin(SUPG)法。

在耦合方程組的求解過程中,不同的因變量可能會用到不同的基函數。一個典型的例子是納維-斯托克斯方程的求解,其中的壓力往往比速度更平滑、更易進行近似。在某類方法中,如果一個耦合方程組中不同的因變量的基函數(以及試函數)屬于不同的函數空間,那么這類方法便稱為混合有限元法。

相關標簽搜索：理論基礎:有限元法的較為詳細的說明 Ansys有限元培訓 Ansys workbench培訓 ansys視頻教程 ansys workbench教程 ansys APDL經典教程 ansys資料下載 ansys技術咨詢 ansys基礎知識 ansys代做 Fluent、CFX流體分析 HFSS電磁分析 Abaqus培訓