基于有限元的密度投影場快速求解研究

陳沁梅，張征宇，周潤，嚴華，祝福順，李果陽

（1.四川大學電子信息學院，成都610065；2.中國空氣動力研究與發(fā)展中心，綿陽621000）

某流場中，密度投影場與擾動引起的光線偏移角之間滿足泊松方程的關(guān)系，且該泊松方程的源項是一系列離散的偏移角，無法利用現(xiàn)有方法直接求解。在有限元法的基礎(chǔ)上，利用神經(jīng)網(wǎng)絡(luò)擬合源項中偏折角與采樣坐標的關(guān)系，得到剖分網(wǎng)格點上的源項；同時，針對神經(jīng)網(wǎng)絡(luò)的加入使得有限元法中單元載荷向量的求解過于耗時的問題，在求解源項時，論文將三角單元整體預測的二重積分表達式近似替換為三角單元結(jié)點預測的常數(shù)表達式。仿真實驗表明，引入神經(jīng)網(wǎng)絡(luò)擬合偏折角，相對于傳統(tǒng)的插值方法，可以取得更高精度的結(jié)果；而且在相同誤差下，提出的算法大大提升了運算速度。將該算法應用到真實流場中，得到的該密度投影場的特性與搭建的真實環(huán)境的結(jié)果相似，進一步說明所提算法的有效性。

密度投影場；泊松方程；神經(jīng)網(wǎng)絡(luò)；有限元；二重積分

0 引言

流場中物理模型的密度投影場信息是飛行器氣動外形優(yōu)化的重要數(shù)據(jù)，因此流動密度投影場的定量分析一直都是流場的一個重要研究方向[1]。目前有效的方法是將背景紋影技術(shù)（Background Oriented Schlieren，BOS）光路和視頻測量技術(shù)融合起來研究試驗中的密度投影場，通過測量標記點光束穿過擾動流場到相機攝影中心的光偏折數(shù)據(jù)，以解決紋影和陰影技術(shù)定量化不足、干涉測量技術(shù)抗干擾能力差和波面?zhèn)鞲衅骺臻g分辨率不夠等問題[2-3]。

密度投影場在數(shù)學上是一個帶有源項的泊松方程，一般使用有限元法進行求解[4]。有限元法是從變分法和有限差分的基礎(chǔ)上發(fā)展而來的，它不但保留了有限差分的離散處理，還借鑒里茲法的試函數(shù)（即下文提到的形函數(shù)）思想[5-6]。早在上世紀四十年代就出現(xiàn)了有限元的計算思想，但隨著高速電子計算機的出現(xiàn)，有限元法才得以真正用于解決工程問題。1956年Turner等人發(fā)表了第一篇有限元法在結(jié)構(gòu)力學中應用的論文。1966年，Zienkiewiez等人第一次提出將有限元法應用于求解位勢流問題[7]。這之后，有限元法被迅速應用到流體力學中[8]?？梢哉J為，有限元法是流體力學中理論研究、實驗分析及解決工程問題的強有力的工具。

考慮到本文所求解方程的特殊性，即常見泊松方程的源項是解析式，而本文需要求解的泊松方程源項是由一系列離散點表示，無法直接利用有限元進行求解。因此，需先擬合源項中偏移角與采樣坐標x,y的關(guān)系，將其化為一般的泊松方程。傳統(tǒng)的方法有曲面擬合和插值，但精度都比較低。由于神經(jīng)網(wǎng)絡(luò)能夠以任意精度逼近連續(xù)函數(shù)，泛化能力強，顯著地優(yōu)于二次曲面擬合結(jié)果[9]。因此，本文采用神經(jīng)網(wǎng)絡(luò)來實現(xiàn)擬合。

1 密度投影場的建模

本文通過將BOS技術(shù)和視頻測量技術(shù)結(jié)合的方式對密度投影場進行建模。在BOS系統(tǒng)中，相機的觀測對象設(shè)置為具有隨機點陣分布的背景，當二者之間存在流場擾動時，相機上背景點陣的像會相較于無擾動時有一定的偏移[10]。由于該偏移量與擾動流場的折射率分布直接相關(guān)，便可定量得到擾動流場的密度分布及通過擾動流場的氣動光學特性信息[11-12]。

圖1 BOS系統(tǒng)光路圖

已知，BOS將待測擾動流場的密度變化與光線的偏折角聯(lián)系在一起，且偏折角與氣體密度存在一定關(guān)系[13]。通過搭建光路，如圖1所示，可獲得光線的偏折角。

式中：空氣中標準狀態(tài)下，KGD的取值為0.226，單位為cm3g，ρ0為流場試驗的總壓。可見，經(jīng)過擾動流場后，光線的偏折角與密度的梯度以及光路長度成正比，兩式微分后相加，可得到如下的含有源項的泊松方程[14]。

它表示密度投影場與x,y方向上偏折角的關(guān)系。而本文的目標，就是求解這一泊松方程。

2 密度投影場求解

2.1 單元剖分

為了更加準確地求解密度投影場的泊松方程，采用三角剖分將待求解區(qū)域離散化，如圖2所示。剖分過程包括三個步驟：首先將待求解區(qū)域( x,y)按照一定間隔k劃分成矩形網(wǎng)絡(luò)，其次將上一步剖分好的矩形網(wǎng)格沿著對角線劃分成一系列三角形單元；最后為每一個三角形單元和其頂點編號，為之后的循環(huán)單元求解做準備[15]。

圖2 三角單元圖

2.2 偏移角擬合

密度投影場的泊松方程源項中偏移角與采樣坐標x,y的關(guān)系的擬合是本文算法的重點。它的精確度直接關(guān)系到最終方程求解的精度。本文利用神經(jīng)網(wǎng)絡(luò)來實現(xiàn)擬合，得到偏移角與采樣坐標x,y的關(guān)系，進而利用訓練得到的網(wǎng)絡(luò)模型預測剖分點對應的偏移角，最后由偏移角求出源項。

圖3 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖

出于性能和效率的考慮，本文采用一個具有三層全連接層的神經(jīng)網(wǎng)絡(luò)來擬合偏移角和采樣坐標x,y的關(guān)系。網(wǎng)絡(luò)結(jié)構(gòu)圖如上圖3所示。首先，將數(shù)據(jù)處理得到的偏移角ε和坐標值( )x,y，作為輸入，利用均方差損失約束網(wǎng)絡(luò)，經(jīng)過多次迭代，可以訓練得到一個偏折角ε在x,y方向上對應于背景板坐標值的擬合網(wǎng)絡(luò)。其中，激活函數(shù)為tanh，優(yōu)化函數(shù)為Adam。網(wǎng)絡(luò)訓練完畢后，將三角網(wǎng)格剖分得到的坐標點輸入網(wǎng)絡(luò)，就可以預測其對應的偏移角。

圖4 三角區(qū)域的近似替代

式中：h為求解步長，其值為1 k。

2.3 方程求解

單元載荷向量的求解與泊松方程的源項有關(guān)。即對于本文而言，要求解單元載荷向量，需先預測三角單元的源項，再對其進行坐標變換，再求解如下所示二重積分：

式中：fi表示一個單元上第i個形函數(shù)與源項的作用，x( λ1,λ2),y( λ1,λ2)表示單元上的點經(jīng)坐標變換后的坐標，|J|表示單元變換的雅克比行列式。

顯然，由于神經(jīng)網(wǎng)絡(luò)預測的加入，單元載荷向量的求解時間效率會有所下降。因此，參考定積分的幾何定義，對其進行優(yōu)化。

由于三角形區(qū)域較小，本文采用結(jié)點坐標確定的平面代替整個三角形單元，如圖4。

即可定義一個三維平面，并用其代替三角單元預測，再對其進行二重積分得到單元載荷向量，如式（6）。

時間上，對于整個單元的源項預測，等價于預測單元上的所有點的源項，且由于一個單元上有3個形函數(shù)，需進行三次單元預測；而對于結(jié)點的源項預測，只需預測所有單元結(jié)點的源項。設(shè)每個單元需預測M個點（該點數(shù)僅與二重積分預置函數(shù)integral的運算規(guī)則，計算機性能有關(guān)，與三角單元大小無關(guān)），剖分點數(shù)為n2，則按照本文的剖分原則，剖分單元有2( n-1)2，求解使得改進前后時間性能相同時的臨界剖分點數(shù)，即求解方程：

合并后，得：

由求根公式求得其正根為：

即當結(jié)點數(shù)大于4時，改進算法均比原方法省時。

此外，當用結(jié)點預測替代單元預測后，二重積分的被積函數(shù)轉(zhuǎn)為顯式方程，即：

積分區(qū)域為直角三角形單元，即{(λ1,λ2)|0≤λ1≤1-λ2,0≤λ2≤1}，可利用積分公式將二重積分化為常數(shù)表達式，在保證誤差不變的情況下，也可提升速度。

3 實驗及結(jié)果分析

為了驗證所提泊松方程解法的精度和時間性能，本文首先對比了在不同疏密的網(wǎng)絡(luò)剖分下，神經(jīng)網(wǎng)絡(luò)與傳統(tǒng)插值方法對同一方程源項的預測效果；再設(shè)計了兩組仿真實驗，利用一個已知真解的泊松方程，在不同的網(wǎng)格剖分精度下，分別測試了原有限元法和本文改進的有限元法的求解精度和求解時間。實驗結(jié)果如表1-表3。進一步，本文也在真實的密度投影場場景下，利用BOS系統(tǒng)得到的數(shù)據(jù)測試了所提方法的有效性，結(jié)果如圖5所示。本文實現(xiàn)算法的環(huán)境是MATLABR2017b以及Python 3.6.5。

3.1 插值法與神經(jīng)網(wǎng)絡(luò)法對比

為了驗證相比于傳統(tǒng)插值算法，神經(jīng)網(wǎng)絡(luò)對源項的預測更精確，本文針對泊松方程：

式中：Ω=( 0 ,10)2,g1( x,y)=-6x-6y，真解為u( x,y)=x3+y3。進行了不同疏密的網(wǎng)絡(luò)剖分，從源項g1( x,y)中產(chǎn)生數(shù)據(jù)對( xn,yn,g1n( xn,yn))，并分出10%作為測試集，其余作為訓練集。

采用L2范數(shù)進行誤差估計[16]，其公式為：

分別用插值和神經(jīng)網(wǎng)絡(luò)的方法對源項進行預測，得到如表1所示的誤差結(jié)果。

從表1中可見，在取樣規(guī)模較大時，隨著網(wǎng)格剖分的精度增加，神經(jīng)網(wǎng)絡(luò)和插值方法的誤差都會減小，但是在同一規(guī)模下神經(jīng)網(wǎng)絡(luò)的誤差遠小于插值方法的誤差說明神經(jīng)網(wǎng)絡(luò)擬合精度遠高于插值方法。

3.2 仿真實驗

由于實際BOS系統(tǒng)測得的流場密度所滿足的泊松方程的源項為一系列離散的點，而不是確切的解析表達式，因此無法求出泊松方程的解析解。為了驗證所提解法的精度和時間性能，本文采用一個已知真解、源項為已知解析式的泊松方程

式中：Ω=( 0 ,1)2,g2( x,y)=2π2sinπx sinπy，真解為u( x,y)=sinπx sinπy。使用與實際BOS系統(tǒng)中相同的數(shù)據(jù)采樣方式從源項g2( x,y)中產(chǎn)生相同數(shù)目的數(shù)據(jù)對( xn,yn,g2n( xn,yn))作為神經(jīng)網(wǎng)絡(luò)的輸入，來擬合得到源項g2( x,y)與坐標點x,y之間的關(guān)系，繼而利用該擬合得到的源項通過原有限元法和本文的方法求解公式（13）所示的泊松方程。

（1）誤差分析

已知，在均勻網(wǎng)格下，結(jié)點z上的誤差為：

式中：u為真實解，uh是步長為h時所求近似解，w(z) h2是精度為O(h4)時的誤差。若剖分間隔k擴大一倍，即h縮小一倍，則：

由公式（14）和公式（15）可得到h 2步長時的可計算誤差為：

因此，針對不同的剖分步長h，我們對比了本文所提方法和原有限元法的可計算誤差，誤差對比結(jié)果如表2所示。

表2 誤差對比

其中，兩種方法不同點在于單元載荷向量求解方法不同。有限元法為對整個三角單元預測進行二重積分；而本文方法為將三角單元結(jié)點預測輸入到二重積分化簡的解析式中。由表2可知，剖分間隔越大，網(wǎng)格剖分得越密，兩種方法的求解誤差均越小，且隨著采樣數(shù)的增大，誤差出現(xiàn)震蕩，最后可認為改善剖分間隔，已對其誤差的減小作用不大。

（2）時間分析

為了驗證所提解法能夠有效地減少離散源項的泊松方程的求解時間，本文在不同的網(wǎng)格剖分間隔k下測試了原有限元法和本文方法的求解時間如表3所示。

表3 時間性能

從表3可知，在實驗中測試的所有的網(wǎng)格剖分精度下，本文所提的解法時間性能大大優(yōu)于原有限元法，而且網(wǎng)格剖分得越細，本文所提解法的時間效率優(yōu)勢越明顯。

結(jié)合表2和表3可看出，本文所提的解法在犧牲了很小的求解精度下可以大大提高求解的時間效率。

3.3 真實實驗

圖5 真實實驗二維結(jié)果圖

真實實驗數(shù)據(jù)是某真實流場提供的。根據(jù)密度投影場模型，如圖1所示，我們可以得到觀測點A坐標( )x,y及有擾動時的點A＇相對無擾動時A的偏移量對應的像素值( )?x,?y。欲求出方程的源項，還需得到背景板上的偏移值，再將坐標值投影到擾動場中。再由式（17）的偏折角公式，就可以得到泊松方程源項中的偏移角。

觀察兩個結(jié)果圖發(fā)現(xiàn)，整個試驗區(qū)域內(nèi)，大部分均與自然來流密度相同，擾動造成的密度變化主要出現(xiàn)在區(qū)域內(nèi)，其中存在兩個峰值（均位于兩張圖的右下角），一個約為自然來流密度的0.8倍，一個約是自然來流密度的1.15倍。說明本論文求解結(jié)果與實驗結(jié)果基本一致，證明了其有效性。

4 結(jié)語

考慮到本文所要求解密度投影場的特殊性，即源項不是解析式，而是一系列離散的數(shù)值時，我們在有限元的基礎(chǔ)上，做了如下改進：

（1）為了利用有限元法直接求解，首先通過神經(jīng)網(wǎng)絡(luò)擬合源項中的偏移角與采樣區(qū)域的x,y坐標的關(guān)系，將神經(jīng)網(wǎng)絡(luò)擬合的結(jié)果作為偏移角的表達式，計算剖分網(wǎng)格點上的偏移角的值，從而得到源項；

（2）神經(jīng)網(wǎng)絡(luò)的加入導致有限元法中單元載荷向量的求解時間性能降低，因此用單元節(jié)點預測近似代替整個單元預測，并將二重積分化為解析式的形式。

為了驗證本文方法的有效性，通過三組實驗可以得到如下結(jié)論：

（1）仿真實驗驗證，近似優(yōu)化使求解在達到相同誤差時，明顯縮減了時間消耗；

（2）真實實驗表明，利用本論文所提方法求解實際流場密度泊松方程的結(jié)果，跟真實環(huán)境下測得結(jié)果相符，且密度投影場的變化反映了擾動場的變化，說明了本文算法的有效性。