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

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

有限元法簡(jiǎn)介

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

舉例來(lái)說(shuō),某函數(shù) u 可能是一個(gè)偏微分方程中的因變量(即溫度、電勢(shì)、壓力等)?？梢愿鶕?jù)下列表達(dá)式,通過(guò)基函數(shù)的線性組合將函數(shù) u 近似為新的函數(shù) uh:

(1)

以及

(2)

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

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

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

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

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

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

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

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

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

(3)

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

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

(4)

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

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

(5)

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

(6)

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

可以通過(guò)傳導(dǎo)熱通量的本構(gòu)關(guān)系來(lái)描述一個(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ù)來(lái)表示的情況下, 該微分方程就被稱為偏微分方程(PDE),這是因?yàn)槊總€(gè)導(dǎo)數(shù)都可能代表(幾個(gè)可能方向中的)某個(gè)方向上的變化。還需注意的是,常微分方程中的導(dǎo)數(shù)是用 d 來(lái)表示的,而偏微分方程中的導(dǎo)數(shù)則是用更卷曲的 ? 來(lái)表示的。

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

(9)

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

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

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

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

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

(10)

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

(11)

(12)

(13)

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

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

(14)

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

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

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

(15)

通過(guò)要求此等式對(duì)希爾伯特空間中的所有 試函數(shù)都成立,可以實(shí)現(xiàn)方程(10)的弱形式公式或稱為變分公式。之所以說(shuō)是“弱”,是因?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)成立是沒(méi)有意義的。

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

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

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

(16)

由此,對(duì)每個(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(或它們的梯度)的乘積。很難用三維空間來(lái)描述試函數(shù)和基函數(shù)的支集,但其二維的類比卻是能夠被可視化的。

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

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

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

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

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

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

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

(19)

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

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

(20)

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

(21)

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

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

(22)

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

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

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

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

不同的單元

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

對(duì)于二維和三維的線性函數(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)格剖分。在從六面體邊界層單元到四面體單元的過(guò)渡中,錐體單元通常被放置在邊界層單元的頂端。

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

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

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

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

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

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

構(gòu)造解方法

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

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

(23)

(24)

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

(25)

(26)

其中,

(27)

這里,

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

(28)

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

(29)

如此,就可以為不同選擇的離散化程度和  計(jì)算出誤差及其模。如果改動(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)單性和普遍性:既可以用于非線性問(wèn)題和時(shí)變問(wèn)題(瞬態(tài)問(wèn)題),也可以用于任何的數(shù)值方法。

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

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

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

網(wǎng)格收斂

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

(30)

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

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

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

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

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

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

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

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

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

其他有限元公式

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

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

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