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

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

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

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

密度投影場；泊松方程；神經網絡；有限元；二重積分

0 引言

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

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

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

1 密度投影場的建模

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

圖1 BOS系統光路圖

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

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

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

2 密度投影場求解

2.1 單元剖分

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

圖2 三角單元圖

2.2 偏移角擬合

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

圖3 神經網絡結構圖

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

圖4 三角區域的近似替代

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

2.3 方程求解

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

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

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

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

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

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

合并后，得：

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

即當結點數大于4時，改進算法均比原方法省時。

此外，當用結點預測替代單元預測后，二重積分的被積函數轉為顯式方程，即：

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

3 實驗及結果分析

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

3.1 插值法與神經網絡法對比

為了驗證相比于傳統插值算法，神經網絡對源項的預測更精確，本文針對泊松方程：

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

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

分別用插值和神經網絡的方法對源項進行預測，得到如表1所示的誤差結果。

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

3.2 仿真實驗

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

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

（1）誤差分析

已知，在均勻網格下，結點z上的誤差為：

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

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

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

表2 誤差對比

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

（2）時間分析

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

表3 時間性能

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

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

3.3 真實實驗

圖5 真實實驗二維結果圖

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

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

4 結語

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

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

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

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

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

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