理論基礎(chǔ):有限元法的較為詳細(xì)的說明

2022-03-27  by:CAE仿真在線  來源:互聯(lián)網(wǎng)

有限元法簡介

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

舉例來說,某函數(shù) u 可能是一個(gè)偏微分方程中的因變量(即溫度、電勢、壓力等)。可以根據(jù)下列表達(dá)式,通過基函數(shù)的線性組合將函數(shù) u 近似為新的函數(shù) uh:

(1)

以及

(2)

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

使用基函數(shù)的線性組合對函數(shù)進(jìn)行逼近的繪圖。

函數(shù) u(藍(lán)色實(shí)線)通過 uh(紅色虛線)進(jìn)行逼近,后者是線性基函數(shù)的線性組合(ψi 用黑色實(shí)線表示)。線性基函數(shù)的系數(shù)由 u0 到 u7 表示。

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

對函數(shù)逼近進(jìn)行離散化。

函數(shù) u(藍(lán)色實(shí)線)通過 uh(紅色虛線)進(jìn)行逼近,后者是線性基函數(shù)的線性組合(ψi 用黑色實(shí)線表示)。線性基函數(shù)的系數(shù)由 u0 到 u7 表示。

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

有限元法的另一個(gè)優(yōu)點(diǎn)是該理論已經(jīng)發(fā)展得較為成熟了,原因就在于偏微分方程問題的數(shù)值表述式和弱表達(dá)式之間的密切關(guān)系(見下面的部分)。例如,當(dāng)數(shù)值模型方程在計(jì)算機(jī)上求解時(shí),該理論在誤差估計(jì)或誤差邊界 估計(jì)方面是較為有效的。

回顧有限元方法的歷史,可知該方法是在 20 世紀(jì) 40 年代初被德裔美國數(shù)學(xué)家 Richard Courant 首次提出的。雖然 Courant 報(bào)道了有限元方法在諸多問題上的應(yīng)用,但幾十年之后該方法才在結(jié)構(gòu)力學(xué)之外的領(lǐng)域獲得了普遍的應(yīng)用,成就了現(xiàn)在的地位。

使用有限元法對輪輞進(jìn)行的結(jié)構(gòu)分析。對輪輞進(jìn)行的結(jié)構(gòu)分析,圖中顯示有限元離散化、應(yīng)力和變形。

代數(shù)方程、常微分方程、偏微分方程和物理定律

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

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

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

(3)

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

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

(4)

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

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

(5)

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

(6)

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

可以通過傳導(dǎo)熱通量的本構(gòu)關(guān)系來描述一個(gè)固體中的熱通量,這也稱為傅里葉定律:

(7)

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

(8)

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

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

(9)

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

類似于上面提到的熱能守恒方程,可以推導(dǎo)出動(dòng)量守恒與質(zhì)量守恒的方程(這兩個(gè)方程構(gòu)成了流體動(dòng)力學(xué)的基礎(chǔ))。此外,亦可以推導(dǎo)出空變與時(shí)變問題中的電磁場和通量方程,從而得到偏微分方程組。

繼續(xù)這一討論,讓我們看看如何從偏微分方程中推導(dǎo)出所謂的弱形式公式

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

假定正在研究的一個(gè)散熱器中的溫度分布由方程(8)給出,但現(xiàn)正處于穩(wěn)定狀態(tài),這就意味著方程(8)中的溫度場的時(shí)間導(dǎo)數(shù)為零。模型域 Ω 的域方程如下:

(10)

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

(11)

(12)

(13)

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

散熱器的數(shù)學(xué)模型。

散熱器數(shù)學(xué)模型的域方程和邊界條件。

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

(14)

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

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

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

(15)

通過要求此等式對希爾伯特空間中的所有 試函數(shù)都成立,可以實(shí)現(xiàn)方程(10)的弱形式公式或稱為變分公式。之所以說是“弱”,是因?yàn)槠浞艑捔?10)的要求,也就是偏微分方程的各項(xiàng)在每一個(gè)點(diǎn)上都必須被明確定義的要求。相反的是,只有在積分時(shí)才要求(14)和(15)是相等的關(guān)系。例如,弱公式化完全允許解的一階導(dǎo)數(shù)不連續(xù),因?yàn)檫@種情況并不妨礙積分。但是,它為二階導(dǎo)數(shù)引入的分布 則并不是普通意義上的函數(shù)。因此,在不連續(xù)點(diǎn)上要求(10)成立是沒有意義的。

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

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

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

(16)

由此,對每個(gè)試函數(shù) ψj 而言,方程(15)的離散化形式即變?yōu)?

(17)

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

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

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

(18)

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

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

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

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

共享兩個(gè)三角形單元的基函數(shù)在二維幾何域中互相重疊。帳篷形狀的線性基函數(shù),在相應(yīng)節(jié)點(diǎn)上的值為 1,在所有其他節(jié)點(diǎn)上的值為 0。兩個(gè)基函數(shù)共享一個(gè)單元時(shí)會發(fā)生重疊。

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

比如說,假設(shè)現(xiàn)在這兩個(gè)基函數(shù)更進(jìn)一步地分開了。這兩個(gè)函數(shù)不共享單元,但它們有一個(gè)共同的單元頂點(diǎn)。如下圖所示,它們不重疊。

只含一個(gè)共同單元頂點(diǎn)的基函數(shù)在二維幾何域中不重疊。共享一個(gè)單元頂點(diǎn)的兩個(gè)基函數(shù)在二維域中不發(fā)生重疊。

當(dāng)這兩個(gè)基函數(shù)重疊時(shí),方程(17)具有非零值,且對系統(tǒng)矩陣的貢獻(xiàn)也是非零的。當(dāng)沒有重疊時(shí),積分為零,因此對系統(tǒng)矩陣的貢獻(xiàn)也為零。

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

使用有限元法在散熱器模型中生成的溫度場近似圖。對散熱器中的溫度場進(jìn)行的有限元近似。

瞬態(tài)問題(時(shí)變問題)

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

(19)

在此,系數(shù) Ti 是時(shí)變函數(shù),而基函數(shù)和試函數(shù)則僅依賴于空間坐標(biāo)。再者,在時(shí)間域上的時(shí)間導(dǎo)數(shù)不是離散的。

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

(20)

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

(21)

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

第二個(gè)公式則基于 t 處的解:

(22)

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

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

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

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

不同的單元

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

對于二維和三維的線性函數(shù),最常見的元素如下圖所示。此圖上圖 給出的是線性基函數(shù)(被定義在三角形網(wǎng)格中,形成了三角形的線性單元)?;瘮?shù)被表示為節(jié)點(diǎn)位置(二維時(shí):x 和 y;三維時(shí):x、y 和 z)的函數(shù)。

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

該示意圖顯示二維和三維線性單元的節(jié)點(diǎn)的幾何與位置。二維和三維線性單元的節(jié)點(diǎn)位置與幾何形狀。

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

與線性單元對應(yīng)的二階單元。二階單元。如果移除灰色節(jié)點(diǎn),便可得到相應(yīng)的巧湊邊點(diǎn)單元。黑色、白色和灰色節(jié)點(diǎn)都存在于拉格朗日單元中。

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

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

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

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

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

構(gòu)造解方法

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

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

(23)

(24)

此方法可用于對改動(dòng)后的問題進(jìn)行求解

(25)

(26)

其中,

(27)

這里,

是可以被自由選擇的一個(gè)解析表達(dá)式。另外,如果

(28)

則  是改動(dòng)后的問題的精確解,此時(shí)可以直接計(jì)算出誤差大小:

(29)

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

目標(biāo)定向的誤差估計(jì)

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

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

網(wǎng)格收斂

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

(30)

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

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

橢圓形膜的結(jié)構(gòu)力學(xué)基準(zhǔn)模型。

橢圓薄膜的基準(zhǔn)模型,其中假設(shè)沿 x 和 y 軸(滾動(dòng)支座)的邊呈對稱分布,并在外部邊上施加載荷。

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

該圖顯示用于二次基函數(shù)的拉格朗日矩形單元。用于二次基函數(shù)的矩形單元。

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

繪圖顯示不同單元和單元尺寸的相對誤差。顯示不同單元和單元尺寸(單元尺寸 = h)的前圖中計(jì)算點(diǎn)處的 σx 的相對誤差。四邊形是指矩形單元,它們可以是線性的,也可以具有二次基函數(shù)。

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

網(wǎng)格自適應(yīng)

在計(jì)算出了這些數(shù)值方程的解 uh 之后,就可以用后驗(yàn) 局部誤差估計(jì)值來創(chuàng)建一個(gè)密度更大的網(wǎng)格,該網(wǎng)格具有較大的誤差。然后可以使用細(xì)化的網(wǎng)格來計(jì)算出第二個(gè)近似解。

下圖描述了一個(gè)被加熱的圓柱體在受到穩(wěn)態(tài)流體流動(dòng)作用下的溫度場。對這一穩(wěn)態(tài)問題進(jìn)行了兩次求解:一次是用基礎(chǔ)網(wǎng)格,另一次是用一個(gè)細(xì)化網(wǎng)格(被基本網(wǎng)格計(jì)算出的誤差估計(jì)所控制)。該細(xì)化網(wǎng)格在溫度和熱通量方面的計(jì)算精度更高,而這一點(diǎn)可能正是該實(shí)例所需要的。

在流體流動(dòng)作用下的加熱氣缸周圍的溫度場,同時(shí)顯示了未經(jīng)網(wǎng)格細(xì)化和經(jīng)過網(wǎng)格細(xì)化的計(jì)算結(jié)果。在流動(dòng)作用下的受熱圓柱體周圍的溫度場計(jì)算結(jié)果,上圖未經(jīng)網(wǎng)格細(xì)化,下圖經(jīng)過了網(wǎng)格細(xì)化。

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

在噴墨打印機(jī)模型仿真中應(yīng)用網(wǎng)格細(xì)化。在一個(gè)瞬態(tài)兩相流問題中,對噴墨打印機(jī)中的一串墨滴進(jìn)行網(wǎng)格細(xì)化。

其他有限元公式

在上述例子中,我們?yōu)榛瘮?shù)和試函數(shù)使用了相同的函數(shù)集來實(shí)現(xiàn)模型方程的離散化。如果一個(gè)有限元公式可以使試函數(shù)不同于基函數(shù),則該公式稱為 Petrov-Galerkin 法。這是一種常用的方法;例如,在解決對流-擴(kuò)散問題的過程中,只會對流線方向進(jìn)行穩(wěn)定化處理。其也被稱為流線迎風(fēng) /Petrov-Galerkin(SUPG)法。

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

COMSOL Multiphysics 仿真軟件中混合有限元方法的設(shè)置。COMSOL Multiphysics 軟件中用于流體流動(dòng)分析的混合單元法的設(shè)置,其中二次形函數(shù)(基函數(shù))用于計(jì)算速度,線性形函數(shù)用于計(jì)算壓力。


開放分享:優(yōu)質(zhì)有限元技術(shù)文章,助你自學(xué)成才

相關(guān)標(biāo)簽搜索:理論基礎(chǔ):有限元法的較為詳細(xì)的說明 Ansys有限元培訓(xùn) Ansys workbench培訓(xùn) ansys視頻教程 ansys workbench教程 ansys APDL經(jīng)典教程 ansys資料下載 ansys技術(shù)咨詢 ansys基礎(chǔ)知識 ansys代做 Fluent、CFX流體分析 HFSS電磁分析 Abaqus培訓(xùn) 

編輯
在線報(bào)名:
  • 客服在線請直接聯(lián)系我們的客服,您也可以通過下面的方式進(jìn)行在線報(bào)名,我們會及時(shí)給您回復(fù)電話,謝謝!
驗(yàn)證碼

全國服務(wù)熱線

1358-032-9919

廣州公司:
廣州市環(huán)市中路306號金鷹大廈3800
電話:13580329919
          135-8032-9919
培訓(xùn)QQ咨詢:點(diǎn)擊咨詢 點(diǎn)擊咨詢
項(xiàng)目QQ咨詢:點(diǎn)擊咨詢
email:kf@www.br5w05v.cn