劉 程， 趙文暢， 陳磊磊， 陳海波＊

（1．中國科學(xué)技術(shù)大學(xué) 近代力學(xué)系 中國科學(xué)院材料力學(xué)行為和設(shè)計(jì)重點(diǎn)實(shí)驗(yàn)室，合肥230027；2．信陽師范學(xué)院 土木工程學(xué)院，信陽464000）

1 引 言

隨著計(jì)算機(jī)技術(shù)的應(yīng)用與普及，CAD與CAE技術(shù)在工程領(lǐng)域得到了長足的發(fā)展。非均勻有理B樣條（NURBS）在過去數(shù)十年成為了幾何造型領(lǐng)域應(yīng)用最廣泛的造型工具，在曲線曲面建模中有突出優(yōu)點(diǎn)?？梢跃_表示二次規(guī)則曲線曲面，從而能用統(tǒng)一的數(shù)學(xué)形式表示規(guī)則曲面與自由曲面；具有可影響曲線曲面形狀的權(quán)因子，使形狀更易于控制和實(shí)現(xiàn)。在CAE領(lǐng)域，運(yùn)用最廣泛的模型描述語言多是網(wǎng)格語言，隨著工程問題呈現(xiàn)出大而復(fù)雜的特點(diǎn)，傳統(tǒng)的網(wǎng)格數(shù)據(jù)越來越龐大，使得產(chǎn)品設(shè)計(jì)與分析割裂開來。Hughes等［1］提出了等幾何分析方法，利用NURBS基函數(shù)替代傳統(tǒng)的拉格朗日插值函數(shù)，統(tǒng)一CAD和CAE的描述語言，建立起產(chǎn)品設(shè)計(jì)與分析的橋梁。

作為經(jīng)典的CAE數(shù)值方法之一，邊界元法在某些特定的領(lǐng)域如無限域問題中展現(xiàn)出了卓越的性能。傳統(tǒng)邊界元法在結(jié)構(gòu)邊界處劃分網(wǎng)格，同樣會(huì)遇到設(shè)計(jì)和分析模型描述語言不融合的問題。近些年，等幾何邊界元法得到了快速發(fā)展，已在彈性力學(xué)［2］、位勢(shì)問題［3，4］、斷裂力學(xué)［5，6］、聲學(xué)［7－9］和形狀優(yōu)化［10，11］等領(lǐng)域取得了重要的研究成果。在聲學(xué)領(lǐng)域，Burton－Miller方 法［12，13］可 以 避 免 邊 界元法在求解外聲場問題中遇到解的非唯一性現(xiàn)象，但該方法具有超奇異積分問題，文獻(xiàn)［7，8］用正則化算子消除了等幾何邊界積分方程中的超奇異性。本文不同于以往，選取奇異性相消技術(shù)［14］，并結(jié)合Cauchy主值積分和Hadamard有限部分積分，推導(dǎo)出直接計(jì)算超奇異積分項(xiàng)的半解析表達(dá)式。

聲學(xué)敏感度［15］表征聲學(xué)物理量對(duì)設(shè)計(jì)變量的變化率，可以用來決定設(shè)計(jì)變量的優(yōu)化方向，因此敏感度信息是基于梯度優(yōu)化方法的重要基礎(chǔ)。傳統(tǒng)網(wǎng)格劃分中，網(wǎng)格節(jié)點(diǎn)的敏感度信息往往難以實(shí)現(xiàn)結(jié)構(gòu)的整體形狀優(yōu)化，而等幾何分析方法將形狀優(yōu)化對(duì)象轉(zhuǎn)化為建模時(shí)的控制點(diǎn)坐標(biāo)，優(yōu)化過程中可以靈活控制幾何形狀，避免了網(wǎng)格重新劃分等問題。本文以NURBS控制點(diǎn)坐標(biāo)為設(shè)計(jì)變量，利用直接微分法［16－18］推導(dǎo)出等幾何敏感度邊界積分方程，直接求出域內(nèi)聲壓對(duì)控制點(diǎn)坐標(biāo)的敏感度，從而找到結(jié)構(gòu)形狀與域內(nèi)聲壓分布之間的對(duì)應(yīng)關(guān)系，為進(jìn)一步的聲學(xué)結(jié)構(gòu)形狀優(yōu)化做好準(zhǔn)備。針對(duì)敏感度分析的超奇異積分問題，將采用等幾何邊界積分方程中相同的處理方式，以保證計(jì)算效率與精度。

2 NURBS基函數(shù)

NURBS插值基函數(shù)通常建立在B樣條基函數(shù)的基礎(chǔ)上，二維形式定義如下［3］，

式中 N 為一維B樣條基函數(shù)，可以看出，Ri，j（ξ，υ）是由兩個(gè)一維B樣條基函數(shù)通過張量乘積的形式組合得到，這兩個(gè)B樣條基函數(shù)分別定義在節(jié)點(diǎn)矢量Ξ＝［ξ1，ξ2，…，ξn＋p＋1］和Υ＝［υ1，υ2，…，υm＋l＋1］上，其中p和l分別為兩個(gè)方向上的基函數(shù)階次，n和m分別為兩個(gè)方向上的基函數(shù)數(shù)目，wi，j為權(quán)重。

NURBS基函數(shù)繼承了B樣條基函數(shù)非負(fù)性、緊支性以及歸一性等良好的數(shù)學(xué)特性，給定控制點(diǎn)Pi，j，則 NURBS曲面可表示為

NURBS可以精確構(gòu)造二次規(guī)則曲面，如典型的球面?？刂泣c(diǎn)Pi，j和權(quán)重wi，j一一對(duì)應(yīng)，可靈活控制結(jié)構(gòu)形狀。另外NURBS曲面的法向?qū)?shù)和雅可比系數(shù)均可以解析表達(dá)，因此幾何信息是完全精確的，有助于提高計(jì)算精度。

3 三維聲場等幾何邊界元法

3．1 等幾何邊界積分方程

三維聲場的控制方程為Helmholtz方程，通過格林第二等式變換，并將源點(diǎn)置于邊界上，得到邊界積分方程

式中 x為源點(diǎn)，y為場點(diǎn)，∮為聲壓，q＝?∮／?n為聲壓法向通量；若邊界光滑，c（x）＝1／2；核函數(shù)G（x，y）＝eikr／（4πr），k為波數(shù)，r＝｜x－y｜；核函數(shù)F（x，y）＝?G（x，y）／?n（y）。采用 Burton－Miller法處理解的非唯一性問題［12］，則法向?qū)?shù)邊界積分方程為

式中 G1（x，y）＝?G（x，y）／?n（x）

F1（x，y）＝?F（x，y）／?n（x）

采取配點(diǎn)法求解方程（4），配點(diǎn)坐標(biāo)為Greville坐標(biāo)［2］。邊界上的∮和q用NURBS基函數(shù)插值得到。

式中

其中 中的V（x，y）依次表示核函數(shù)G，F(xiàn)，G1和F×表示一個(gè)積分單元。用Burton－Miller法組合式（6，7），得到線性代數(shù)方程組

式中 H和G為積分項(xiàng)求得的系數(shù)矩陣，Φ～和q～分別為配點(diǎn)處聲壓系數(shù)與聲壓通量系數(shù)列向量，將所有未知值移到方程左邊，已知值移到方程右邊，求解該方程，即可得到邊界未知物理量。

3．2 奇異積分處理

式（6，7）的積分可由 Gauss－Legendre積分計(jì)算，但當(dāng)源點(diǎn)參數(shù)坐標(biāo)滿足×［υe，υe＋1］時(shí)，式中的積分出現(xiàn)奇異現(xiàn)象，這時(shí)需要特殊的處理方式才能正確計(jì)算這些奇異項(xiàng)。

式（6）中核函數(shù)G為弱奇異性，極坐標(biāo)變換即可消除奇異性，自然坐標(biāo) （ξ，υ）與極坐標(biāo) （ρ，θ）對(duì)應(yīng)關(guān)系為

坐標(biāo)轉(zhuǎn)換完成后，式中積分項(xiàng)簡寫為

式（11）右邊可直接用Gauss－Legendre積分計(jì)算。

核函數(shù)F和G1形式上具有強(qiáng)奇異性，但均含有?r／?n項(xiàng)，該項(xiàng)在r→0時(shí)也趨于0，因此也是弱奇異性的，相應(yīng)積分經(jīng)極坐標(biāo)變換后可直接由Gauss－Legendre積分計(jì)算。

由于F1為超奇異，單純的極坐標(biāo)變換無法消除奇異性，借助Guiggiani等［14］提出的結(jié)合Cauchy主值積分和Hadamard有限部分積分的奇異相消技術(shù)，可直接計(jì)算該類積分。核函數(shù)F1中超奇異項(xiàng)為nl（x）nl（y）／4πr3，極坐標(biāo)轉(zhuǎn)換后該項(xiàng)可展開為

根據(jù)Cauchy主值積分和Hadamard有限部分積分

法，包含F(xiàn)1的奇異積分項(xiàng)可表示為

式（13）各項(xiàng)均為無奇異項(xiàng)，可由 Gauss－Legendre積分直接計(jì)算，其中β（θ）和γ（θ）的表達(dá)式可參考文獻(xiàn)［14］，接下來確定f1（θ）和f2（θ）的具體表達(dá)式即可。

將NURBS插值基函數(shù)在源點(diǎn)處泰勒展開為

同理nl（x）nl（y）J（y）在源點(diǎn)處泰勒展開為

另外，r－3可展開為

式中 S3（θ）和S2（θ）的具體表達(dá)式與文獻(xiàn)［14］相同。將式（14～16）代入式（12）等號(hào)左邊得，

4 對(duì)控制點(diǎn)的敏感度分析

4．1 等幾何敏感度邊界積分方程

邊界物理量的敏感度由NURBS函數(shù)插值得

將式（3，4）對(duì)設(shè)計(jì)變量求導(dǎo)，并將式（2，5，19）代入，得到離散形式的等幾何敏感度邊界積分方程

用Burton－Miller法組合式（20，21），得到線性代數(shù)方程組，寫成矩陣形式為

4．2 奇異積分處理

對(duì)于式（20，21）的奇異積分，同樣采用3．2節(jié)的處理方式。包含核函數(shù)G，F(xiàn)，，，G1以及的奇異積分項(xiàng)均為弱奇異，采用極坐標(biāo)變換方式即可得到正確數(shù)值積分結(jié)果。包含F(xiàn)1的超奇異積分項(xiàng)處理方式同3．2節(jié)，包含的超奇異積分項(xiàng)的具體處理方式如下。

式中各項(xiàng)可由Gauss－Legendre積分直接計(jì)算，β（θ）和γ（θ）與3．2節(jié)一致，u1（θ）和u2（θ）同樣由泰勒展開的方式確定。

將式（26）左邊各項(xiàng)作泰勒展開，與文獻(xiàn)［14］主要不同之處是r·的展開表達(dá)式，

式中 Ai，Bi和A均為θ的函數(shù)，與文獻(xiàn)［14］一致，C＝AiBi。A·i和B·i分別為Ai和Bi對(duì)設(shè)計(jì)變量的敏感度，均可通過NURBS基函數(shù)插值顯式表達(dá)。S0和S1分別為式（27）中ρ對(duì)應(yīng)階次前的系數(shù)。另外，r－4可展開為［14］

式中

將式（14，15，27，28）代入式（25）左邊，最終得到u1（θ）和u2（θ）的表達(dá)式為

5 數(shù)值算例

5．1 脈動(dòng)球的聲輻射

以典型的脈動(dòng)球輻射算例來驗(yàn)證本文推導(dǎo)的公式。本文所有的計(jì)算過程均在臺(tái)式計(jì)算機(jī)上用Fortran90編程實(shí)現(xiàn)，計(jì)算機(jī)配置為Intel Core i7 CPU和16GB內(nèi)存，為了提高效率與保證精度，本文用于Gauss－Legendre積分的點(diǎn)數(shù)目在非奇異積分情況下為6×6，奇異積分時(shí)為10×10。

聲傳播介質(zhì)為空氣，密度為1．2kg／m3，聲波速度為340．0m／s。如圖1所示，該球面上各點(diǎn)以速度vn做同振幅同相位的振動(dòng)，球體振動(dòng)會(huì)在周圍介質(zhì)中引起聲輻射。設(shè)已知邊界條件為Neumann邊界條件，邊界聲壓法向通量值q＝1509．6i。整個(gè)球面通過NURBS曲面構(gòu)造，半徑為1．0m?？刂泣c(diǎn)信息如圖2所示，共26個(gè)控制點(diǎn)，參數(shù)ξ與υ方向的階次均為2，節(jié)點(diǎn)矢量分別為Ξ＝｛0，0，0，0．5，0．5，1，1，1｝和 Υ ＝｛0，0，0，0．25，0．25，0．5，0．5，0．75，0．75，1，1，1｝。

選取xoy平面內(nèi)半徑為2．0m的圓上等分的180點(diǎn)為參考點(diǎn)，將本文算法計(jì)算的數(shù)值解與脈動(dòng)球聲輻射解析解取相對(duì)誤差，以考察本文算法的精度。圖3為采用等幾何邊界元法時(shí)該算例計(jì)算誤差隨頻率的變化關(guān)系，BM代表Burton－Miller法，自由度數(shù)為121，考察頻率為50Hz～350Hz?？梢钥闯?，沒有使用Burton－Miller法時(shí)，在k＝π和k＝2π處的計(jì)算結(jié)果嚴(yán)重失真；而Burton－Miller法的計(jì)算精度盡管有所降低，但在整個(gè)頻段上均可以計(jì)算準(zhǔn)確，證明該方法可以克服解的非唯一性。

圖1 NURBS建模的脈動(dòng)球Fig．1 NURBS model of pulsating sphere

圖2 NURBS球的控制點(diǎn)信息Fig．2 Control points of the NURBS sphere

圖3 等幾何邊界元計(jì)算誤差隨頻率變化的關(guān)系Fig．3 Relative error of IGA BEM simulations with respect to frequencies

設(shè)計(jì)變量取圖2中控制點(diǎn)P22的y坐標(biāo)，在球體中其初始值為1．0m，考慮結(jié)構(gòu)不出現(xiàn)交叉重疊的現(xiàn)象，設(shè)該設(shè)計(jì)變量取值恒為正值。由于設(shè)計(jì)變量的變化，結(jié)構(gòu)將會(huì)成為不規(guī)則構(gòu)型，因此聲場聲壓敏感度值不存在解析解，本文借助有限差分法（finite difference method）來驗(yàn)證本文敏感度數(shù)值解的正確性。分析聲頻率為100Hz，計(jì)算自由度取121，有限差分步長為0．001m。圖4給出了點(diǎn)（0，2，0）處的聲壓敏感度實(shí)部與虛部值隨控制點(diǎn)P22的y坐標(biāo)的變化關(guān)系，控制點(diǎn)P22的y坐標(biāo)范圍為0．5m～1．5m?？梢钥闯觯疚挠弥苯游⒎址ㄓ?jì)算的結(jié)果與有限差分法的結(jié)果吻合良好，由于積分方程中的超奇異積分非常容易影響精度，而本文采取的奇異相消技術(shù)給出了半解析的直接計(jì)算公式，從直接微分法的結(jié)果來看，很好地消除了超奇異積分的影響。

5．2 不同形狀脈動(dòng)結(jié)構(gòu)的聲輻射

圖4 點(diǎn)（0，2，0）處的聲壓敏感度值隨控制點(diǎn)P22的y坐標(biāo)變化關(guān)系Fig．4 Pressure sensitivities at point（0，2，0）with respect to the ycoordinate of P22

圖5 不同幾何形狀脈動(dòng)結(jié)構(gòu)Fig．5 Different shapes of pulsating structures

NURBS構(gòu)造的構(gòu)型有眾多控制點(diǎn)，每個(gè)控制點(diǎn)的取值不同，結(jié)構(gòu)的形狀就會(huì)不同。因此，本文在圖2脈動(dòng)球控制點(diǎn)的基礎(chǔ)上，改變一些控制點(diǎn)的坐標(biāo)取值，構(gòu)造出不同幾何形狀的脈動(dòng)結(jié)構(gòu)，如圖5所示。保持幾種構(gòu)型的幾何節(jié)點(diǎn)矢量和權(quán)重一致，控制點(diǎn)坐標(biāo)取值如下。（a）保持球體模型不變。（b）P22坐標(biāo)?。?，0．5，0），其他控制點(diǎn)坐標(biāo)不變。（c）P22坐標(biāo)?。?，1．5，0），其他控制點(diǎn)坐標(biāo)不變。（d）P00坐標(biāo)?。?，0，1．5），P22坐標(biāo)?。?，1．5，0），P26坐標(biāo)取（0，－1．5，0），P48坐標(biāo)取（0，0，1．5），其他控制點(diǎn)坐標(biāo)不變。值得注意的是，若控制點(diǎn)重合，則重合點(diǎn)坐標(biāo)也隨之改變。

設(shè)這些結(jié)構(gòu)表面各點(diǎn)同樣以速度vn做同振幅同相位的振動(dòng)，邊界條件不變，以考察域內(nèi)聲場聲壓以及敏感度值的分布情況。選取半徑為2．0m的球面作為考察面，分析該球面上的聲壓分布，以及該球面上的聲壓模值對(duì)P22的y坐標(biāo)的敏感度值分布。分析頻率為100Hz，結(jié)構(gòu)自由度取121。圖6給出了脈動(dòng)結(jié)構(gòu)各個(gè)形狀下考察面上的聲壓模值分布，不同構(gòu)型的聲壓值分布顯然完全不同，除了數(shù)值差異較大外，各個(gè)構(gòu)型上出現(xiàn)聲壓最大與最小值的位置也不同，如圖6（b，c）所示，因此若要對(duì)某目標(biāo)區(qū)域的聲壓值進(jìn)行調(diào)控，在一定范圍內(nèi)可以通過改變結(jié)構(gòu)形狀來完成。形狀I(lǐng)為規(guī)則球體，因此圖6（a）考察面上聲壓模值完全相同。形狀I(lǐng)V稍復(fù)雜，因此圖6（d）中聲壓模值分布與結(jié)構(gòu)上凸出的部分相關(guān)。圖7給出了各個(gè)形狀下，考察面上聲壓模值對(duì)P22的y坐標(biāo)的敏感度值分布?？梢钥闯觯疾烀嫔下晧耗Ｖ档拿舾卸确植贾饕獮閿?shù)值上的差異，但均呈現(xiàn)對(duì)稱性帶狀分布。這些敏感度值表征了考察面上聲壓對(duì)設(shè)計(jì)變量的變化率。

5．3 水下大尺度殼體的聲散射

圖6 不同形狀脈動(dòng)結(jié)構(gòu)考察面上聲壓模值分布Fig．6 Distribution of sound pressure modulus on the reference surface with respect to different structure shapes

大尺度殼體模型的聲散射數(shù)值仿真具有重要的意義，如水下潛艇聲散射特性的數(shù)值仿真可以有效模擬潛艇的聲學(xué)性能，為潛艇的聲學(xué)設(shè)計(jì)提供必要的參考。2001年德國FWG提出標(biāo)準(zhǔn)潛艇BeTSSi－Sub（Benchmark Target Strength Simulation Submarine）的概念，Nell等［19］在這基礎(chǔ)上作出了相關(guān)改進(jìn)，艇身主要尺寸如圖8所示。艇艏為半橢球體，艇體為圓柱體，艇艉為圓錐體，笛卡爾坐標(biāo)系原點(diǎn)位于圓柱艇體的中心。將圖8的模型用NURBS建模，控制點(diǎn)信息如圖9所示，參數(shù)ξ與υ方向的階次均為2，節(jié)點(diǎn)矢量分別為Ξ＝｛0，0，0，0．25，0．25，0．5，0．5，0．75，0．75，1，1，1｝和 Υ ＝｛0，0，0，0．25，0．25，0．5，0．5，0．75，0．75，1，1，1｝。

圖7 不同形狀脈動(dòng)結(jié)構(gòu)考察面上聲壓模值的敏感度值分布Fig．7 Distribution of sound pressure sensitivities on the reference surface with respect to different structure shapes

圖8 NURBS構(gòu)建的簡化BeTSSi模型Fig．8 Simple BeTSSi model built by NURBS

圖9 簡化BeTSSi的NURBS控制點(diǎn)信息Fig．9 NURBS control points of the simple BeTSSi

考慮入射波為平面波，沿y軸正向入射，幅值為1．0Pa，水中聲速1524．0m／s。為了驗(yàn)證本文算法對(duì)大尺度殼體聲散射仿真的正確性與有效性，首先用常量三角形單元離散該模型，形成103124個(gè)單元，將該方法得到的計(jì)算數(shù)值作為參考解。圖10給出了本文算法計(jì)算點(diǎn)聲壓相對(duì)于參考解的誤差，計(jì)算點(diǎn)等間距分布于xoy平面內(nèi)距離艇軸心點(diǎn)6．0m的圓上。圖10橫坐標(biāo)為計(jì)算點(diǎn)的極坐標(biāo)角度，分析頻率為100Hz?？梢钥闯觯S著自由度的增加，常量三角形單元邊界元相對(duì)誤差逐漸減小并趨于穩(wěn)定，因此取103124個(gè)單元時(shí)的解作為參考解是可行的。本文算法在自由度相對(duì)較少的情況下，取得了與常量三角形單元邊界元一致的計(jì)算精度，如本文算法881自由度與常量三角形單元5418自由度時(shí)的相對(duì)誤差較為接近。針對(duì)大尺度殼體的聲散射問題，本文算法利用NURBS建模，幾何信息精確，因此離散規(guī)模遠(yuǎn)小于常量三角形單元的規(guī)模。

圖10 等幾何邊界元法相比常量三角形邊界元法的相對(duì)計(jì)算誤差Fig．10 Relative error of IGA BEM simulations relative to constant triangle BEM

圖11 點(diǎn)（6，0，0）處的聲壓敏感度值隨控制點(diǎn)P04的x坐標(biāo)變化關(guān)系Fig．11 Pressure sensitivities at point（6，0，0）with respect to the xcoordinate of P04

設(shè)計(jì)變量取圖9中P04控制點(diǎn)的x坐標(biāo)，分析頻率為100Hz，計(jì)算自由度為315。圖11給出了點(diǎn)（6，0，0）處的聲壓敏感度隨設(shè)計(jì)變量的變化關(guān)系，有限差分法步長為0．001m，設(shè)計(jì)變量取值為2．5m～5．0m?？梢钥闯?，本文用直接微分法計(jì)算的結(jié)果與有限差分法的結(jié)果吻合良好，說明本文敏感度算法同樣適用于大尺度模型。隨著結(jié)構(gòu)形狀的變化，計(jì)算點(diǎn)聲壓敏感度呈現(xiàn)出不同的取值，這樣的梯度信息可以為結(jié)構(gòu)形狀優(yōu)化提供設(shè)計(jì)方向。圖12和圖13分別展示了xoy平面內(nèi)目標(biāo)區(qū)布情況，目標(biāo)區(qū)域?yàn)?0×20m的正方形。圖12目標(biāo)區(qū)域的聲壓模值分布表明艇身入射面聲壓較強(qiáng)，達(dá)到1．6Pa，而背面聲壓可降至0．4Pa。圖13表明控制點(diǎn)的影響區(qū)域是有限的，離控制點(diǎn)越近，聲壓對(duì)控制點(diǎn)變化越敏感，因此基于本文算法的后續(xù)形狀優(yōu)化設(shè)計(jì)必須考慮控制點(diǎn)的布局。

圖12 目標(biāo)區(qū)域聲壓模值分布（Pa）Fig．12 Distribution of sound pressure modulus on the reference surface

圖13 目標(biāo)區(qū)域聲壓敏感度值分布Fig．13 Distribution of sound pressure sensitivities on the reference surface

6 結(jié) 論

本文推導(dǎo)了基于NURBS基函數(shù)插值的三維聲場等幾何邊界積分方程及其敏感度邊界積分方程，針對(duì)邊界積分方程中的奇異性，結(jié)合Cauchy主值積分和Hadamard有限部分積分，成功運(yùn)用奇異相消技術(shù)給出超奇異積分的直接計(jì)算公式，并據(jù)此編寫了計(jì)算程序。采用脈動(dòng)球輻射算例和水下大尺度殼體的聲散射問題，驗(yàn)證了本文算法和程序的正確性與有效性，并且基于脈動(dòng)球給出其他不同形狀的結(jié)構(gòu)，初步考察了結(jié)構(gòu)形狀對(duì)聲場聲壓分布的影響，為下一階段結(jié)構(gòu)形狀優(yōu)化工作做了理論準(zhǔn)備。