基于水平集的二維半空間聲屏障優化設計

張 放，姜富航，陳海波

（中國科學技術大學中國科學院材料力學行為與設計重點實驗室，合肥 230027）

近年來隨著工業的發展，噪聲控制已經逐漸成為產品設計、城市規劃等領域重點考慮的課題。在眾多降噪手段中，利用隔聲設備隔斷聲波的傳播路徑是一種簡單有效且被廣泛使用的方法，比如都市中常見的聲屏障。而隔聲體的隔聲效果與其形狀、分布及材料屬性有關，同時針對不同的聲場也不盡相同。為了達到預期的降噪效果，可通過聲學結構優化來設計出滿足要求的隔聲體。

以聲屏障為例，現有的研究包括參數優化[1]、形狀優化[2]及表面吸聲材料的拓撲優化[3]等。將聲屏障簡化成二維結構，在給定聲源的情況下，分別優化各尺寸參數[1]、邊界形狀[2]以及表面吸聲材料分布[3]。但以上方法都需要給定初始結構，且僅能對結構進行局部調整，而無法充分利用設計空間。若不限制障礙物的具體形狀，則可利用拓撲優化的方法對結構本身進行設計。拓撲優化的最終目的是得到材料的最佳分布，而不受結構外形的約束，因此可充分利用設計區間，具有更高的優化自由度。

水平集方法是一種用于追蹤界面演化的數值技術，可以方便地表征物體的拓撲結構改變[4]。與傳統拓撲優化相比，其優點在于避免了以單元密度作為設計變量而導致的如棋盤格、中間密度等數值問題。當前階段，水平集理論經過發展也產生了多個分支，如標準水平集、參數化水平集和反應擴散方程水平集等方法。盡管求解手段各異，水平集方法的核心思想都是用高一維度的水平基函數的0等高線表示結構邊界，也稱為零水平集，通過迭代更新水平集函數實現結構邊界的追蹤。

反應擴散方程水平集由Yamada 等[4]提出，根據虛擬界面能的思想，利用反應擴散方程實現水平集函數的更新，并被應用于最小柔度優化問題[5]中。這一方法避免了標準水平集方法中求解Hamilton-Jacobi 方程中遇到的諸多數值困難，但需要使用Fréchet 導數作為拓撲導數進行計算。在聲學問題中，Matsumoto等[6]基于伴隨變量法給出了拓撲導數的高效計算方法。Ⅰsakari 等[7]利用快速多極邊界元法進行了三維聲學結構拓撲優化。Sugihara 等[8]則對Robin邊界條件下二維結構進行優化。最近，Gao等[9]針對邊界積分形式的目標函數，給出了二維結構聲學拓撲導數的具體推導，并用于聲子晶體的優化研究中。以上研究均是在全空間自由聲場中進行，而實際生活中地面對聲波的散射往往不容忽略。在計算聲屏障的聲學性能時，通常可用半空間格林函數模擬地面對聲波的反射效果[2]。

本文將針對半空間的聲學問題，以邊界元法對聲屏障的降噪性能進行計算分析，并借助Yamada等[4]提出的反應擴散方程水平集方法，對聲屏障截面進行優化設計。為了避免優化中產生孔洞而破壞結構的拓撲關系，僅取零水平集的主體部分作為其截面的實際形狀，最終達成形狀優化的目的。同時為了區別于Gao等[9]的研究，將以域積分作為目標函數的定義形式，并給出對應的伴隨方程。該方法適用于任何域積分形式的目標函數，也可用于更大規模的半空間聲學優化問題中。

1 二維半空間聲學計算

理想流體介質中頻域聲場由Helmholtz 方程及邊界條件控制：

以二維聲學邊界元法對環境中的聲場進行求解，入射波作用下的邊界積分方程如式（2）所示：

式中：c（x）根據點x在域內、光滑邊界上和域外分別取1、1/2和0；p（x）為聲壓，q（x）為聲通量，且q（y）=為格林函數，也稱基本解；q*（x,y）為基本解的法向導數，且q*（x,y）=為點y處的外法向矢量；pi（x）為入射聲壓。

區別于全空間格林函數，對于二維半空間問題，格林函數包含鏡面項：

式中：k為波數；r=|x-y|代表場點x與源點y之間的距離；r′=|x′-y|代表場點x關于地面的鏡像點x′與源點y之間的距離；H0（1）為第一類0 階Hankel函數。

數值計算中，首先將結構邊界Γ劃分成一系列的邊界單元，再對式（2）進行離散，并計算出每個邊界單元上的積分系數后組裝成：

式中：p和q分別是邊界上節點對應的聲壓和聲通量組成的列向量；H和G是對應的系數矩陣；pi是入射聲壓在邊界節點上值的列向量。

2 水平集方法

2.1 水平集函數

水平集函數φ（x,t）是設計域內點x坐標（x，y）的函數。在二維問題中，φ（x,t）本身對應第三個維度的坐標，也即高度z。以其零水平集作為結構的邊界，以區分材料域與聲場區域：

式中：Ω為聲場域為材料域；Γ為兩域交界面。

2.2 優化模型

定義目標函數如式（6）所示：

式中：f是聲壓的函數，用于建立域積分形式的統一格式。優化的目的是得到最優的材料區域配置，使得目標函數達到全局或局部的最小值。水平集方法通過改變水平集函數φ（x,t），更新結構邊界，進而影響到觀測區域的聲壓分布。因此若希望使目標函數逐步減小，則需要使水平集函數沿著其負梯度方向演化。Yamada等[4]提出的反應擴散方程如下所示：

式中：T（x,t）是目標函數的拓撲導數，表明當前位置的拓撲改變對目標函數的影響。K＞0 是正則化參數，負號則滿足了優化的負梯度要求。第二項里，Δφ是曲率項，用于確保水平集函數的光滑性，系數τ用于人為指定迭代過程中水平集函數變化的幾何復雜度。

將時間項φ（x,t）/?t以有限差分格式表示，φ（x,t）以有限元法進行離散，最終可建立起水平集函數φ（x,t）的迭代關系式，詳細推導可參考Otomori等[5]對結構柔度優化問題給出的樣例程序，但須注意他們使用的數值模型里每個單元長度固定，即最終模型會隨單元數增大而變大，并不適用于指定區域內的優化計算。

2.3 拓撲導數

拓撲導數T（x,t）定義為，聲場在位置x處產生一個小孔后，目標函數的改變量與小孔面積（二維情況下）比值的極限，如式（8）所示，用于表示目標函數受當前位置聲場拓撲變化的影響。

式中：表示增設小孔后的目標函數值；JΩ則表示小孔產生前的目標函數值；ε為小孔半徑；最終結果T（x,t）為小孔圓心x處的拓撲導數值。當T（x,t）＞0 時，在x處鋪設材料將會使目標函數增大，反之則減小，因此可以作為φ（x,t）隨時間的迭代依據。

依據Matsumoto等[6]給出的拓撲導數推導方法，計算域積分形式的目標函數，拓撲導數最終表達式與Gao 等[9]全空間問題給出的表達式一致，如式（9）所示。

式中：p（x）為聲壓，由微分方程（1）控制；λ（x）為伴隨變量，由微分方程（10）控制；λ,j（x）與p,j（x）分別表示λ（x）和p（x）的梯度分量；Re 表示取復數的實部。此時，由于目標函數選取為特定區域的聲壓積分，伴隨場的控制方程中存在域內的非齊次項，如式（10）所示。

聲場和伴隨場的微分方程均通過邊界元方法求解。在積分計算系數矩陣時，盡管式（1）和式（10）并不完全相同，但由它們推導得到的邊界積分方程系數矩陣H和G仍是一致的，只是λ（x）計算中需要將式（2）中的入射波自由項替換成與f有關的域積分而后根據邊界條件組裝系數矩陣及自由項，就可求得p（x）和λ（x）的邊界值。

域內點的計算由積分方程給出，如式（11）和式（12）所示。p,j（x）與λ,j（x）則可通過將其對坐標方向求導得到。

最后根據式（9），可直接得到點x處的拓撲導數值，并結合反應擴散方程對φ（x,t）進行更新。

在電氣工程中，首先要做的是確定電氣工程中電壓等級的選擇，以獲得節能效果。例如，當前通用電氣工程中的中壓配電電壓為10kV或者20kV，低壓配電電壓為220V或380V[6]。對于電壓的選擇，需要進行全面考量。電壓選擇與許多因素有關：一是國家電網在當地的供電設備的特點；二是電氣工程中的電源電路數量；三是電力用戶的計算電壓容量；四是電氣工程所在地區的發展水平。高壓電網和低壓電網在選擇時也有不同：選擇低壓配電網的情況是用戶的電力用于低功耗；選用高壓配電網的情況是用戶的計算能力，電氣設備，單電源儀表電源的距離和計算負荷超過一定的范圍內。在電壓選擇過程中，需要參考上述要求選擇正確的電壓以避免浪費。

2.4 優化流程

整個優化的流程如圖1所示。首先給出初始水平集函數，其零水平集對應結構的初始構型。在半空間邊界元計算時，需要對結構進行單元離散，用于積分系數的計算及系統方程的組裝，并結合邊界條件計算出所有邊界值。根據域內點積分表達式計算觀測點處的聲壓值之后，即可由式（9）得到設計區域內有限元節點處的拓撲導數值。進而依據式（7）所示的反應擴散方程更新水平集函數，對比目標函數的改變量，以判斷是否收斂。

圖1 水平集優化流程圖

若并未收斂，則需要根據更新后的水平集函數重新提取其零水平集，并劃分網格。值得注意的是，在結構的聲學拓撲優化中，很容易產生多個獨立懸浮的微小散射體，盡管確有一定的隔聲效果，但難以實際應用，因此可以考慮僅保留面積最大的主體部分作為實際聲屏障的截面。這一行為限制了結構拓撲關系的改變，最終表現為截面的形狀優化。

3 數值算例

3.1 Γ形聲屏障優化模型

以Γ 形聲屏障的優化為例，y=0 m 作為地面，用位于（0,1）處的單位點聲源模擬公路上直行的汽車噪聲。則環境中的入射波場可考慮為：

式中：k為波數，r代表當前位置x與點聲源之間的距離。

聲屏障寬為0.2 m，高為3 m，左下角置于（5,0）處。觀測點設置為x∈（6.1,16），y∈（0.1,3）之間均勻分布的100 個點，具體模型如圖2 所示，計算頻率為400 Hz，聲速為340 m/s。

圖2 聲屏障截面優化模型

半空間問題中，聲屏障底部與地面貼合，而地面作為鏡面，視作剛性邊界條件。在進行邊界元分析時，通過選取半空間格林函數作為基本解，模擬聲波在鏡面上的反射過程，可以抵消聲屏障在地面上的邊界積分，從而只需要對地面之上的部分進行建模計算[2]。

設計域高為4 m，寬為2 m。在求解反應擴散方程時，給定設計域邊界處φ（x,t）=-1 作為有限元法的邊界條件。

目標函數J可表示為：

則在微分方程（10）中有：

3.2 拓撲導數驗證

首先為展示式（9）所示的拓撲導數在半空間情況下的正確性，分別選擇設計域內y=3.2 m 和4.0 m直線上的有限元節點，利用式（8）以有限差分方法進行驗證。小孔半徑設置為ε=1×10-4m，對比結果繪圖如下：

由圖3 可見，式（9）計算結果與差分方法得到的拓撲導數值完全一致，證明了基于伴隨變量法計算得到的拓撲導數的正確性，其在半空間問題中也能夠很好地反映出當前位置的拓撲變化對目標函數的影響。

圖3 半空間拓撲導數差分驗證

3.3 半空間水平集優化

將式（9）得到的拓撲導數作為水平集函數更新的梯度依據，可以使目標函數逐步減小。每次迭代中基于面積對零水平集進行的過濾處理，則能夠保證結構拓撲關系不發生改變，從而只調整其外形。本小節將采用了過濾處理的形狀優化結果與不采用過濾處理的拓撲優化結果進行對比，并利用商業軟件COMSOL計算出優化前后的聲場分布。

水平集函數的初始配置由離散后的有限元節點指定，在此將其指定為二值函數，位于聲屏障內部賦值為1，否則為-1。同時，為了固定聲屏障的位置，將0.5 m 作為聲屏障的基座尺寸，即y＜0.5 m 的有限元節點在迭代中固定不變，取為初始值。

依據式（9）得到的拓撲導數，在反應擴散方程式（7）中對水平集函數進行迭代更新。參數K依據Otomori 等[5]給出的樣例程序取為面積積分與靈敏度積分之比，τ取為1×10-5。優化步長dt初值取為0.06，并在目標函數增大較為明顯時適當減小。形狀優化和拓撲優化的目標函數的迭代收斂曲線如圖4 所示，縱軸以對數坐標表示。

圖4 目標函數迭代歷史

觀察迭代歷史，兩種處理下目標函數均能夠保持穩定下降直至收斂。最終收斂值相近，都相比初值降低了兩個多數量級，對應的結構對比如圖5 所示。同時與Liu 等[2]對聲屏障頂端的等幾何優化將近80%的降噪幅度相比，優化效果更為明顯。

圖5 優化結果對比

觀察圖5，拓撲優化結果中出現了多個微小且獨立的孔洞，而聲屏障主體結構僅發生了幾處局部的變形。形狀優化結果中，結構拓撲關系并未改變，但整體形狀與初始結構有了很大的區別。

在實際背景下，每一塊孤立的二維散射體都代表一個以此為截面的無限長桿狀結構。拓撲優化中出現的多個不規則孔洞對應著各自的桿狀體，其與主體結構分離，懸浮于空中，將給生產加工帶來極大的不便。而形狀優化后的結構外形更為復雜，為得到相當的隔聲性能使原本平整的結構外形變得起伏不定。

理論上拓撲優化由于設計自由度更高而應該能得到更好的優化結果，但若是希望只調整結構的外形，那么進行面積過濾的形狀優化則是一種值得嘗試的方法。這兩種處理并沒有優劣之分，只是按照實際需求進行取舍。最后利用商業軟件COMSOL對優化前后結構周圍的半空間聲場進行計算，對應的聲壓級云圖如圖6所示。

圖6 聲壓級云圖對比

圖6（a）是初始結構聲場分布，圖6（b）和圖6（c）分別代表拓撲與形狀優化后的聲場分布。與初始結構的聲場對比，觀測區域云圖由以40 dB 為主的金黃色區域變成以0～20 dB為主的藍綠色區域，同時也出現了負聲壓級的部分深藍色區域，整體聲壓級有了顯著的降低。

4 結 語

對于二維半空間問題而言，通過伴隨變量法計算的拓撲導數表達式與全空間一致，但需要利用半空間格林函數求解聲壓、伴隨變量及其梯度。在提取零水平集作為結構邊界時，可以取面積最大的部分作為目標結構，以限制其拓撲關系的改變。最終隨反應擴散方程更新結構外形，同樣能實現相當的優化效果。

基于反應擴散方程的水平集方法根據拓撲導數對水平集函數進行調整，不依賴于迭代歷史，僅針對當前結構進行優化，穩定可靠。對聲屏障截面進行優化設計時，降噪效果極為顯著，能夠為實際工程中的優化設計提供清晰具體的理論指導。