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

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

有限元法簡(jiǎn)介

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

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

(1)

以及

(2)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)

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

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

(4)

如此,該固體中的溫度通過(guò)一個(gè)代數方程(4)來(lái)表示,其中的某個(gè)時(shí)間值 t1 就會(huì )有一個(gè)對應時(shí)間的溫度值 T1。

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

(5)

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

(6)

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

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

(7)

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

(8)

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

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

(9)

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

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

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

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

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

(10)

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

(11)

(12)

(13)

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

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

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

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

(14)

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

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

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

(15)

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

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

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

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

(16)

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

(17)

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

散熱器模型的有限元離散化。來(lái)自之前的散熱器模型圖的有限元離散化。

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

(18)

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

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

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

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

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

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

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

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

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

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

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

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

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

(19)

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

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

(20)

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

(21)

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

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

(22)

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

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

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

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

不同的單元

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

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

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

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

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

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

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

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

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

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

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

構造解方法

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

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

(23)

(24)

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

(25)

(26)

其中,

(27)

這里,

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

(28)

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

(29)

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

目標定向的誤差估計

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

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

網(wǎng)格收斂

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

(30)

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

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

橢圓形膜的結構力學(xué)基準模型。

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

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

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

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

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

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

網(wǎng)格自適應

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

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

在流體流動(dòng)作用下的加熱氣缸周?chē)臏囟葓?chǎng),同時(shí)顯示了未經(jīng)網(wǎng)格細化和經(jīng)過(guò)網(wǎng)格細化的計算結果。在流動(dòng)作用下的受熱圓柱體周?chē)臏囟葓?chǎng)計算結果,上圖未經(jīng)網(wǎng)格細化,下圖經(jīng)過(guò)了網(wǎng)格細化。

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

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

其他有限元公式

在上述例子中,我們?yōu)榛瘮岛驮嚭瘮凳褂昧讼嗤暮瘮导瘉?lái)實(shí)現模型方程的離散化。如果一個(gè)有限元公式可以使試函數不同于基函數,則該公式稱(chēng)為 Petrov-Galerkin 法。這是一種常用的方法;例如,在解決對流-擴散問(wèn)題的過(guò)程中,只會(huì )對流線(xiàn)方向進(jìn)行穩定化處理。其也被稱(chēng)為流線(xiàn)迎風(fēng) /Petrov-Galerkin(SUPG)法。

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

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


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

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

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

全國服務(wù)熱線(xiàn)

1358-032-9919

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