喻 敏 王志鴻 許文輝 章 軒

(武漢理工大學(xué)交通學(xué)院1) 武漢 430063) (中國艦船研究設(shè)計(jì)中心2) 武漢 430063)

0 引 言

流固耦合(fluid-structure interaction,FSI)是描述流體動(dòng)力學(xué)定律和結(jié)構(gòu)力學(xué)定律之間的多物理場(chǎng)耦合.流固耦合的關(guān)鍵是實(shí)現(xiàn)精確的耦合分析，其中分區(qū)耦合較適用于數(shù)值計(jì)算，在CFD數(shù)值計(jì)算中應(yīng)用廣泛.分區(qū)耦合涉及流體到結(jié)構(gòu)的網(wǎng)格節(jié)點(diǎn)插值、荷載和變形傳輸問題，即邊界數(shù)據(jù)傳遞.利用徑向基函數(shù)(radial basis function,RBF)插值方法進(jìn)行數(shù)據(jù)傳遞[1]，可有效解決邊界網(wǎng)格不匹配而導(dǎo)致節(jié)點(diǎn)數(shù)據(jù)難以直接傳遞的問題[2]，特別是對(duì)復(fù)雜非線性流固耦合問題具有很好的適用性，具有良好的工程應(yīng)用前景.

針對(duì)非線性流固耦合問題，基于黏性理論數(shù)值求解流體場(chǎng)時(shí)，流體網(wǎng)格節(jié)點(diǎn)具有較強(qiáng)非線性，使得流體網(wǎng)格節(jié)點(diǎn)很容易出現(xiàn)數(shù)值噪聲(包括異常值).對(duì)含有數(shù)值噪聲的流場(chǎng)網(wǎng)格節(jié)點(diǎn)進(jìn)行徑向基函數(shù)插值，可能導(dǎo)致插值結(jié)果的失真.近年來，關(guān)于邊界數(shù)據(jù)傳遞的徑向基函數(shù)插值研究主要集中在徑向基核函數(shù)和算法效率方面[3-4].

常用的數(shù)值噪聲處理方法有統(tǒng)計(jì)分析法、特征提取法和密度聚類法等[5-6].統(tǒng)計(jì)分析法，適用于噪聲的判斷，噪聲處理需配合其他算法；特征提取法，如小波去噪，則過分依賴于閾值的選??；密度聚類算法，如DBSCAN等，計(jì)算復(fù)雜度高，對(duì)輸入?yún)?shù)敏感.上述方法均不適用于流固耦合邊界數(shù)據(jù)的去噪計(jì)算，而圖像濾波方法無需先驗(yàn)知識(shí)，原理簡單，應(yīng)用廣泛、靈活，在三維數(shù)據(jù)去噪方面效果明顯[7-10].文中嘗試將圖像處理領(lǐng)域的濾波算法應(yīng)用于CFD流場(chǎng)數(shù)據(jù)處理，分析高斯濾波、雙邊濾波和導(dǎo)向?yàn)V波三種濾波方法的去噪效果，確定合適的濾波方法，對(duì)CFD計(jì)算數(shù)據(jù)進(jìn)行去噪處理，并將其運(yùn)用到流固耦合邊界數(shù)據(jù)的插值計(jì)算中，進(jìn)而形成基于徑向基函數(shù)插值的邊界數(shù)據(jù)傳遞方法.

1 邊界數(shù)據(jù)的徑向基函數(shù)插值

1.1 流固耦合數(shù)據(jù)傳遞矩陣

在流固耦合問題的流體動(dòng)力學(xué)計(jì)算中，由于結(jié)構(gòu)網(wǎng)格與流場(chǎng)網(wǎng)格的差異，需要在流固耦合邊界上進(jìn)行數(shù)據(jù)傳遞.流固兩相交界處一般滿足虛功相等原理，即：

(1)

式中：W為虛功；uf和us分別為流固邊界處的虛位移；ff和fs分別為流體表面壓力和固體表面壓力。以H表示傳遞矩陣，經(jīng)推導(dǎo)，可以得到：

(2)

由于流固耦合數(shù)據(jù)傳遞是雙向的過程，在此為了簡化敘述，僅以流體區(qū)域信息傳遞到固體結(jié)構(gòu)區(qū)域?yàn)槔?。假定流固耦合邊界處兩?cè)分別存在nf(歐拉點(diǎn)數(shù))和ns(拉格朗日點(diǎn)數(shù))個(gè)節(jié)點(diǎn)，以Zf和Zs分別表示流體和固體節(jié)點(diǎn)坐標(biāo)，數(shù)據(jù)傳遞中引入臨時(shí)傳遞變量Uf，則徑向基函數(shù)插值的基本形式為

(3)

(4)

式中：λ、α均為徑向基函數(shù)中的系數(shù)矩陣，式中徑向基函數(shù)φ的變量是歐拉點(diǎn)(即流體域上的點(diǎn))之間的距離.則待求系數(shù)：

(5)

對(duì)于流體傳遞給固體的情況，有臨時(shí)傳遞變量：

(5)

(6)

式中：徑向基函數(shù)φ的變量是拉格朗點(diǎn)(即固體結(jié)構(gòu)域上的點(diǎn))與歐拉點(diǎn)(即流體域上的點(diǎn))之間的距離。由以上公式可得：

(8)

由此，流體區(qū)域數(shù)據(jù)交換到固體結(jié)構(gòu)區(qū)域的傳遞矩陣為

(9)

同理，固體區(qū)域數(shù)據(jù)交換到流體區(qū)域的傳遞矩陣為

(10)

由式(9)～式(10)，結(jié)合本節(jié)中討論的數(shù)據(jù)交換能量守恒，可實(shí)現(xiàn)流固耦合邊界上速度、壓力等高效精確的傳遞.從推導(dǎo)過程可知，徑向基插值對(duì)邊界數(shù)據(jù)傳遞十分必要.

1.2 徑向基函數(shù)插值

在流固耦合邊界數(shù)據(jù)傳遞過程中，基于徑向基函數(shù)的插值不依賴網(wǎng)絡(luò)節(jié)點(diǎn)的拓?fù)浣Y(jié)構(gòu)，適用于大量散亂數(shù)據(jù)問題的求解.理論上，徑向基函數(shù)表示為一類以歐式距離為變量的函數(shù)集合，其數(shù)學(xué)表達(dá)式為

(11)

式中：f(x),x∈R為給定高維空間的未知函數(shù)；λi為待求的加權(quán)系數(shù)；Φ(x)為徑向基函數(shù)；?！瑇-xi‖為插值點(diǎn)x到被插值點(diǎn)xi的歐式距離；n為插值點(diǎn)個(gè)數(shù).對(duì)條件正定的徑向基函數(shù)，增加一個(gè)多項(xiàng)式，以確保方程解唯一.多項(xiàng)式pi(x)需滿足如下約束條件.

(12)

由式(11)～式(12)，可得線性方程組：

(13)

有ΘW=Ω，因此為求系數(shù)矩陣W，需要求解系數(shù)矩陣Θ-1，若插值點(diǎn)x存在噪聲，系數(shù)矩陣Θ易出現(xiàn)變態(tài)，導(dǎo)致徑向基函數(shù)插值結(jié)果的失真，甚至無法求解.因此，有必要對(duì)插值點(diǎn)進(jìn)行去噪處理，以提高徑向基函數(shù)插值的穩(wěn)健性.

2 濾波算法

2.1 高斯濾波

一維高斯函數(shù)為

(14)

高斯濾波使用的核函數(shù)為x和y兩個(gè)方向上一維高斯的乘積，且兩個(gè)維度上的標(biāo)準(zhǔn)差σ相同，其中心點(diǎn)始終為濾波窗口中心坐標(biāo)，故取均值μ=0，因此，濾波所用到的核函數(shù)為

(15)

高斯濾波器寬度是由參數(shù)σ控制的，σ越大，高斯濾波器的頻帶就越寬，平滑程度就越大，因此可以通過調(diào)節(jié)參數(shù)σ的取值，在過平滑和欠平滑之間取得折衷.然而，由于高斯濾波的高斯核只考慮了數(shù)據(jù)間的空間距離，沒有考慮數(shù)據(jù)取值上的差異，因而容易丟失數(shù)據(jù)的細(xì)節(jié)信息.

2.2 雙邊濾波

雙邊濾波是一種非線性的濾波方法，綜合了高斯濾波器和α-截尾均值濾波器的特點(diǎn)，同時(shí)考慮濾波數(shù)據(jù)的空間信息及其數(shù)值上的相似性，以達(dá)到保邊去噪的目的，一定程度上彌補(bǔ)了高斯濾波的缺陷.雙邊濾波的核函數(shù)為

w(i,j,k,l)=

(16)

式中：(i,j)為濾波窗口中數(shù)據(jù)的坐標(biāo)；(k,l)為濾波窗口中心點(diǎn)坐標(biāo)；I為對(duì)應(yīng)數(shù)據(jù)的大小；σd、σr均為高斯函數(shù)的標(biāo)準(zhǔn)差，但在數(shù)值上一般不等.

則w(i,j,k,l)=d(i,j,k,l)×r(i,j,k,l)，即雙邊濾波的核函數(shù)為兩個(gè)子函數(shù)的乘積，d(i,j,k,l)代表數(shù)據(jù)的空間域信息，決定濾波器模板系數(shù)的空間鄰近因子，其大小隨著坐標(biāo)點(diǎn)與中心點(diǎn)之間距離的增加而減??；r(i,j,k,l)代表數(shù)據(jù)取值的范圍域信息，決定濾波器模板系數(shù)的數(shù)據(jù)數(shù)值相似因子，其大小隨著數(shù)據(jù)之差的增大而減小.

現(xiàn)考慮兩種極限，①濾波窗口中數(shù)據(jù)與中心數(shù)據(jù)基本相同，鄰域內(nèi)各數(shù)據(jù)的差趨近于0，有相似因子r(i,j,k,l)≈1，此時(shí)空間域權(quán)重起主要作用，相當(dāng)于對(duì)數(shù)據(jù)進(jìn)行高斯濾波；②窗口中其他數(shù)據(jù)與中心數(shù)據(jù)相差極大的情況下，相似因子r(i,j,k,l)→0，此時(shí)雙邊濾波的核函數(shù)w(i,j,k,l)≈0，不對(duì)數(shù)據(jù)進(jìn)行濾波操作，從而達(dá)到保留數(shù)據(jù)細(xì)節(jié)信息的目的.也正因此，當(dāng)出現(xiàn)高強(qiáng)度噪聲時(shí)，濾波算法的范圍域核函數(shù)權(quán)值的計(jì)算將受到影響，導(dǎo)致其去噪能力下降，甚至出現(xiàn)“梯度反轉(zhuǎn)”現(xiàn)象[12].

2.3 導(dǎo)向?yàn)V波

根據(jù)文獻(xiàn)[10]，導(dǎo)向?yàn)V波需滿足以下條件，

1)輸出數(shù)據(jù)Q與原始數(shù)據(jù)P盡可能的相似，即(Q|P)取得最小.

2)輸出數(shù)據(jù)Q與引導(dǎo)數(shù)據(jù)I的梯度相似，有Q=aI，即滿足Q=aI+b.

現(xiàn)假設(shè)在濾波窗口ωk內(nèi)輸出數(shù)據(jù)Q與引導(dǎo)數(shù)據(jù)I滿足線性關(guān)系，有qi=akIi+bk，?i∈ωk，其中：k表示窗口ωk的中心位置，ak、bk在窗口中為常數(shù).為滿足第一個(gè)條件，通過求解該線性函數(shù)的系數(shù)，使得輸出值q與輸入值p之間差值最小化，即存在

上式可通過最小二乘求解，可得：

(17)

(18)

(19)

(20)

由此可見，在ε≠0的情況下，若引導(dǎo)數(shù)據(jù)I和輸入數(shù)據(jù)P相同，導(dǎo)向?yàn)V波可在處理數(shù)據(jù)的同時(shí)保留數(shù)據(jù)的細(xì)節(jié)信息.此外，由于輸出數(shù)據(jù)與引導(dǎo)數(shù)據(jù)梯度相似，濾波結(jié)果不會(huì)出現(xiàn)雙邊濾波的“梯度反轉(zhuǎn)”問題.

2.4 仿真驗(yàn)算分析

為了比較上述三種濾波算法的濾波效果，采用二元函數(shù)z=xexp(-x2-y2)進(jìn)行驗(yàn)算.對(duì)函數(shù)進(jìn)行加噪，結(jié)果見圖1a)，利用三種濾波算法對(duì)其進(jìn)行去噪，比較濾波前后徑向基函數(shù)插值誤差，由此選取最適合的濾波方法對(duì)CFD數(shù)據(jù)進(jìn)行去噪.

仿真過程中，各濾波器參數(shù)的取值經(jīng)過多次驗(yàn)算獲得：高斯濾波使用5×5窗口、標(biāo)準(zhǔn)差σ=1.5的核函數(shù)對(duì)數(shù)據(jù)進(jìn)行濾波處理，結(jié)果見圖1b).

雙邊濾波同樣取5×5窗口，核函數(shù)中取σd=1.5，σr=0.9.此外，由于雙邊濾波要求輸入數(shù)據(jù)P∈[0,1]，需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，預(yù)處理過程如下：

1)記Pmin為輸入數(shù)據(jù)P最小值，將數(shù)據(jù)P中所有元素加上|Pmin|，得新數(shù)據(jù)P′=P+|Pmin|，記Pmax為其最大值。

2)對(duì)數(shù)據(jù)P′做歸一化處理，得最終輸入雙邊濾波器的輸入數(shù)據(jù)P″.

同理，對(duì)濾波后數(shù)據(jù)進(jìn)行恢復(fù)，最終得到雙邊濾波的結(jié)果見圖1c).

導(dǎo)向?yàn)V波將原始數(shù)據(jù)同時(shí)作為輸入數(shù)據(jù)和引導(dǎo)數(shù)據(jù)進(jìn)行濾波，同樣選取5×5的濾波窗口，參數(shù)ε取1.5，濾波結(jié)果見圖1d).

圖1 加噪數(shù)據(jù)與三種濾波方法所得數(shù)據(jù)對(duì)比

由圖1可知，高斯濾波對(duì)邊緣數(shù)據(jù)野值的處理效果最差，且與雙邊濾波一樣，在梯度變化大的區(qū)域(如濾波數(shù)據(jù)中心部分)濾波效果差，后者主要由“梯度反轉(zhuǎn)”現(xiàn)象所致，而導(dǎo)向?yàn)V波并沒有出現(xiàn)此類現(xiàn)象，在處理此類數(shù)據(jù)時(shí)濾波效果較為理想.從算法計(jì)算效率看，見表1。高斯濾波計(jì)算時(shí)間最短，這得益于其簡單的計(jì)算原理；同時(shí)，盡管導(dǎo)向?yàn)V波和雙邊濾波均考慮了數(shù)據(jù)細(xì)節(jié)信息的特性，導(dǎo)向?yàn)V波計(jì)算所用時(shí)間仍在可接受范圍內(nèi)，而雙邊濾波頗為耗時(shí)，不適用于處理大數(shù)據(jù)，如流場(chǎng)數(shù)值計(jì)算.在平均絕對(duì)誤差(MAE)上，導(dǎo)向?yàn)V波數(shù)值最大，其原因在于為比較三種算法的濾波效果，而選取了同樣大小的窗口，多次仿真計(jì)算可知，減小濾波窗口可進(jìn)一步提高導(dǎo)向?yàn)V波的濾波效果，如導(dǎo)向?yàn)V波在3×3窗口中MAE=1.894，而高斯濾波和雙邊濾波已在5×5窗口取得最優(yōu)濾波效果.

表1 三種濾波方法計(jì)算效率對(duì)比

將濾波前后的數(shù)據(jù)分別進(jìn)行徑向基函數(shù)插值計(jì)算，比較插值精度.插值核函數(shù)采用多元二次核函數(shù)，徑向基函數(shù)插值平均絕對(duì)誤差結(jié)果見表2。噪聲對(duì)插值精度有明顯影響；濾波后的插值精度都得到一定程度的提高，其中導(dǎo)向?yàn)V波后邊界數(shù)據(jù)的插值效果得到大幅度改善，為三種濾波算法中的最優(yōu)，因此本文將采用導(dǎo)向?yàn)V波對(duì)CFD數(shù)據(jù)進(jìn)行濾波處理.

表2 五類數(shù)據(jù)插值誤差對(duì)比

3 導(dǎo)向?yàn)V波流場(chǎng)數(shù)據(jù)處理應(yīng)用

3.1 流場(chǎng)流速數(shù)據(jù)濾波效果對(duì)比

文中以CFD數(shù)值計(jì)算所得某一截面上的速度場(chǎng)數(shù)據(jù)為例進(jìn)行濾波計(jì)算，見圖2a)，速度場(chǎng)以矢量圖的形式給出，圖2b)為由矢量圖轉(zhuǎn)換得到的可視化圖形，由此可見流場(chǎng)數(shù)據(jù)與圖像像素值之間存在著相似性.基于這一方面的考慮，本文將導(dǎo)向?yàn)V波方法用于流場(chǎng)數(shù)據(jù)的去噪處理上，以提高徑向基函數(shù)插值的精度.

文獻(xiàn)[13]采用移動(dòng)最小二乘算法對(duì)耦合前CFD數(shù)據(jù)進(jìn)行去噪，取得了較好效果，并已應(yīng)用于實(shí)際的自主CFD計(jì)算平臺(tái).下面將本文的研究對(duì)象——導(dǎo)向?yàn)V波方法，與移動(dòng)最小二乘方法在流場(chǎng)數(shù)據(jù)去噪方面進(jìn)行對(duì)比，結(jié)果見圖3.結(jié)合圖2a)速度場(chǎng)數(shù)據(jù)的變化，可見移動(dòng)最小二乘算法雖能有效剔除流場(chǎng)數(shù)據(jù)存在的噪點(diǎn)，但其并不考慮鄰域(濾波窗口)內(nèi)數(shù)據(jù)的變化趨勢(shì)，對(duì)連續(xù)變化物理場(chǎng)(如速度場(chǎng))的描述顯然不如導(dǎo)向?yàn)V波，而導(dǎo)向?yàn)V波所得結(jié)果更符合實(shí)際速度場(chǎng)的變化規(guī)律，這使得相較于傳統(tǒng)的移動(dòng)最小二乘算法，導(dǎo)向?yàn)V波方法在CFD數(shù)據(jù)去噪方面更具優(yōu)勢(shì).

圖2 實(shí)驗(yàn)數(shù)據(jù)及其模信息

圖3 移動(dòng)最小二乘與導(dǎo)向?yàn)V波去噪效果對(duì)比

3.2 流固耦合邊界數(shù)據(jù)插值驗(yàn)證

為對(duì)比導(dǎo)向?yàn)V波前后邊界數(shù)據(jù)的徑向基函數(shù)插值精度，以2.4中導(dǎo)向?yàn)V波前后數(shù)據(jù)為例，分別在不同流體網(wǎng)格密度條件下，對(duì)3種密度的邊界數(shù)據(jù)進(jìn)行插值，以此分析導(dǎo)向?yàn)V波方法對(duì)流固耦合邊界數(shù)據(jù)插值效果的影響。仿真中，根據(jù)流體網(wǎng)格密度分成2組：G1，G2，每組網(wǎng)格中又有著不同密度的待插值邊界數(shù)據(jù)點(diǎn)，見圖4，由此可得六種組合方式。為分析導(dǎo)向?yàn)V波方法在流固耦合邊界數(shù)據(jù)插值精度上的影響，表3～4為導(dǎo)向?yàn)V波前后對(duì)應(yīng)網(wǎng)格密度下邊界數(shù)據(jù)的插值誤差對(duì)比.

圖4 不同邊界點(diǎn)密度網(wǎng)格示意圖

表3 濾波前數(shù)據(jù)在不同邊界點(diǎn)密度下插值誤差

表4 濾波后數(shù)據(jù)在不同邊界點(diǎn)密度下插值誤差

由表3～4可知，無論是對(duì)濾波前還是濾波后的邊界數(shù)據(jù)進(jìn)行插值，G2分組所得插值誤差都明顯大于G1分組。根據(jù)數(shù)據(jù)插值或擬合理論，插值誤差和數(shù)據(jù)密度存在一定的關(guān)系，當(dāng)填充距離很小，即數(shù)據(jù)密度很大時(shí)，插值誤差也很小.

一般地，用填充距離h描述數(shù)據(jù){xj}的密度：

(21)

可以證明[14]：

(22)

同時(shí)，在每一分組內(nèi)部，不同的邊界數(shù)據(jù)密度對(duì)應(yīng)的插值誤差十分接近，這是因?yàn)椴逯稻戎饕c插值數(shù)據(jù)密度有關(guān)，而與被插值點(diǎn)密度無關(guān)。另外，對(duì)比表3～4誤差，濾波后數(shù)據(jù)的插值誤差明顯小于原始數(shù)據(jù)的插值誤差，可知導(dǎo)向?yàn)V波處理有效提高了徑向基函數(shù)插值的精度，對(duì)實(shí)現(xiàn)高效的流固耦合邊界數(shù)據(jù)傳遞具有積極作用.

4 結(jié) 論

1)導(dǎo)向?yàn)V波既能保留數(shù)據(jù)的細(xì)節(jié)信息，又具有較高的計(jì)算效率，濾波效果較高斯濾波和雙邊濾波方法更優(yōu)。

2)相較于移動(dòng)最小二乘算法，導(dǎo)向?yàn)V波在處理流場(chǎng)數(shù)據(jù)的過程中，綜合考慮了鄰域內(nèi)數(shù)據(jù)的變化，更適用于連續(xù)變化物理場(chǎng)數(shù)據(jù)的處理。

3)對(duì)數(shù)據(jù)進(jìn)行導(dǎo)向?yàn)V波處理后，邊界數(shù)據(jù)的徑向基函數(shù)插值精度明顯提高，且差值精度受邊界數(shù)據(jù)密度的影響較小.

綜上，利用導(dǎo)向?yàn)V波對(duì)數(shù)據(jù)進(jìn)行去噪后插值，可提高徑向基函數(shù)插值精度，有利于實(shí)現(xiàn)高效的流固耦合邊界數(shù)據(jù)傳遞.