結(jié)構(gòu)聲相互作用的二維結(jié)構(gòu)相互作用問題的一種多圓形隨機(jī)耦合方法

結(jié)構(gòu)聲相互作用的二維結(jié)構(gòu)相互作用問題的一種多圓形隨機(jī)耦合方法

0聲相互作用問題的數(shù)值求解人們一直關(guān)注混合體介質(zhì)中彈性結(jié)構(gòu)的振動(dòng)和聲源的傳播。在這種結(jié)構(gòu)中，潛水員和魚雷等水下結(jié)構(gòu)是典型的結(jié)合系統(tǒng)應(yīng)用形式。這類結(jié)構(gòu)不僅需要足夠的強(qiáng)度、剛度和穩(wěn)定性，還需要保證低的輻射噪聲、分散噪聲和自噪聲水平。另一方面，由于彈性結(jié)構(gòu)周圍流的存在，結(jié)構(gòu)和水體之間存在強(qiáng)烈的流固耦合作用，因此很難對(duì)理論分析和研究方法進(jìn)行。因此，研究復(fù)雜的水下結(jié)構(gòu)的結(jié)構(gòu)聲音具有很大的理論和現(xiàn)實(shí)意義。目前，結(jié)構(gòu)-聲音應(yīng)用的解決方案主要分為分析法和數(shù)值法。第一種方法通常僅適用于相對(duì)簡(jiǎn)單規(guī)則的彈性結(jié)構(gòu)，如球殼、無限長(zhǎng)柱、橢圓形殼等。bjarnason和aciebach使用frrier變換和行波展開法分析了兩個(gè)簡(jiǎn)單分支的聲音輸出，以便在非中、高頻率下通過frrier變換和光譜行動(dòng)法進(jìn)行分析。對(duì)于復(fù)雜結(jié)構(gòu)的聲音耦合，模型方法通常使用數(shù)值方法進(jìn)行。主要數(shù)值方法包括：殼的限制法、限制法、限制元法、限制元法和限制元法。這些方法具有共同特點(diǎn)，即新城模型結(jié)構(gòu)的聲音耦合，它可以用來解決復(fù)雜形狀結(jié)構(gòu)的聲音耦合。然而，由于數(shù)值方法本身的局限性和高頻率的較大干擾（包括質(zhì)量影響），數(shù)值方法僅適用于低頻范圍。對(duì)于具有相同性質(zhì)的復(fù)雜結(jié)構(gòu)的聲音矩陣，數(shù)值方法通常采用數(shù)值方法。主要數(shù)值方法包括：限制元法、限制元法、限制元法、限制元法、限制元法和限制元法。這些方法具有共同特點(diǎn)，即新城模型結(jié)構(gòu)與水體之間的相互作用稱為附加水質(zhì)量矩陣（或矩陣的組合）的形式。然而，由于數(shù)值方法本身的局限性和高頻光束的傳遞矩陣（包括質(zhì)量和衰減的影響），數(shù)值方法僅適用于低頻范圍。為了消除可變厚度波束運(yùn)動(dòng)方程的方程，不能將水動(dòng)壓以作為增加水位的形式引入數(shù)學(xué)中。因此，解算方法采用振動(dòng)函數(shù)迭代法，計(jì)算精度和效率低下。鑒于以上求解方法的不足,本文基于傳遞矩陣法和虛擬源強(qiáng)模擬技術(shù)提出了一種求解在諧激勵(lì)作用下二維結(jié)構(gòu)-聲相互作用問題的半數(shù)值半解析法.對(duì)于任意形狀彈性環(huán)建立了沿周向曲線坐標(biāo)的非齊次一階線性常微分方程組,并在足夠小的積分步長(zhǎng)內(nèi),將非齊次一階線性微分方程組進(jìn)行齊次擴(kuò)容化為齊次一階線性微分方程組,以便用精細(xì)積分法求解.對(duì)于外聲場(chǎng),本文基于虛擬源強(qiáng)模擬技術(shù)采用多圓形虛擬源強(qiáng)配置方案,并在每一條圓形配置曲線上將未知的源強(qiáng)密度函數(shù)用Fourier級(jí)數(shù)展開,同時(shí)結(jié)合快速Fourier逆變換法,提出了一種高精度、高效率求解任意形狀二維孔穴Helmholtz外問題的決速算法.在耦合方程的求解方面,本文根據(jù)疊加原理,將外激勵(lì)力和虛擬源強(qiáng)的Fourier級(jí)數(shù)展開項(xiàng)作為廣義力分別作用在彈性環(huán)上,借助齊次擴(kuò)容方法和精細(xì)積分法求得彈性環(huán)的狀態(tài)向量,再利用流固交接條件和最小二乘法直接建立了耦合系統(tǒng)的求解方程.最后,本文給出的幾個(gè)算例說明了本文方法的有效性.1任何形狀彈性環(huán)的狀態(tài)微分方程組及其序列擴(kuò)張的詳細(xì)方法1.1彈性環(huán)的齊次容錯(cuò)算法考慮一浸沒在流體介質(zhì)中的任意形狀彈性環(huán),橫截面積和慣性矩分別為A和I,其所界外域記為Ω(圖1),微元弧段的受力情況如圖2所示.當(dāng)彈性環(huán)在諧激勵(lì)外力作用下,即F(β,t)=f(β)exp(iωt)(其中i=√-1?ω表示外激勵(lì)圓頻率),由圖2根據(jù)平衡條件可直接導(dǎo)出運(yùn)動(dòng)微分方程{dΝdβ+1RcQ+ρAω2v=0-1RcΝ+dQdβ+ρAω2w+(f-pf)=0dΜdβ-Q=0(1)式中,β表示彈性環(huán)周向曲線坐標(biāo),pf(β)表示作用在環(huán)上的聲壓幅值.v(β)和w(β)表示切向和法向位移幅值,Rc(θ)和ρ分別為彈性環(huán)的曲率半徑和質(zhì)量密度.由廣義Hook定理,彈性環(huán)的內(nèi)力-位移關(guān)系為Ν=Κ(dvdβ+wRc),Φ=vRc-dwdβ,Μ=DdΦdβ(2)式中,Φ(θ)表示彈性環(huán)橫截面繞中性軸的轉(zhuǎn)角幅值,K=EA和D=EI分別為彈性環(huán)的抗拉剛度和抗彎剛度.為便于分析,引入下列無量綱物理量ˉv=vRa?ˉw=wRa?ˉΦ=Φ?ˉΝ=ΝΚ?ˉQ=QΚ?ˉΜ=RaΜD(3)式中,Ra=12π∫2π0√R2+(dRdθ)2dθ為與彈性環(huán)具有相同周長(zhǎng)圓的半徑.將(3)式代入(1)、(2)式中,可寫成下列一階常微分方程組的形式ddβ{ˉvˉwˉΦˉΜ-ˉQˉΝ}=1Ra[0A12000A16A210A23000000A34000000A4500A52000A56A61000A650]{ˉvˉwˉΦˉΜ-ˉQˉΝ}+1Κ{0000f(β)0}-1Κ{0000pf(β)0}(4)(4)式即為任意形狀彈性環(huán)的狀態(tài)微分方程組,其中:{Ζ}=[ˉvˉwˉΦˉΜ-ˉQˉΝ]Τ稱為彈性環(huán)橫截面的狀態(tài)向量,系數(shù)矩陣[Aij]=[A](i,j=1,2,…,6)中的非零元素為A12=-RaRc,A16=1,A21=RaRc,A23=-1?A34=1,A45=-ΚR2aD?A52=λ2,A56=-RaRc?A61=-λ2,A65=RaRc其中,頻率參數(shù):λ2=ρAω2R2aΚ.對(duì)于任意形狀的彈性環(huán),其曲率半徑Rc為曲線坐標(biāo)β(或極角θ)的函數(shù),因而(4)式是一非齊次變系數(shù)線性微分方程組,同時(shí)由于作用在彈性環(huán)上的表面聲壓pf(β)不能事先給定,所以對(duì)(4)式進(jìn)行直接求解是非常困難的.為了解決以上問題,我們提出了一種求解非齊次變系數(shù)線性微分方程組的齊次擴(kuò)容精細(xì)算法.1.2類型化并求解由(4)式的系數(shù)矩陣[A]可以看出,其中非零元素隨曲線坐標(biāo)β的變化關(guān)系僅取決于彈性環(huán)的曲率半徑Rc(β),因此,當(dāng)積分步長(zhǎng)取足夠小時(shí),曲率半徑Rc(β)可作為常數(shù)處理(即將矩陣[A]視為常數(shù)矩陣),這能保證要求的計(jì)算精度.另外,在這個(gè)足夠小的積分步長(zhǎng)內(nèi)將非齊次項(xiàng)用多項(xiàng)式展開同樣也能獲得很好的逼近精度.不失一般性令非齊次微分方程組為下列形式ddβ{Ζ}=1Ra[A(β)]{Ζ}+{F(β)}(5)當(dāng)積分步長(zhǎng)取足夠小時(shí),上述非齊次項(xiàng)可近似寫為(二項(xiàng)式逼近){F(β)}={F0}+{F1}β+{F2}β2(6)式中,{F0},{F1}和{F2}為常數(shù)向量,不同的積分段,其值不同.將(6)式代入(5)式中,并作簡(jiǎn)單整理,得非齊次微分方程組的齊次擴(kuò)容形式ddβ{{Ζ}1ββ2}=[1Ra[A]{F0}{F1}{F2}000001000020]{{Ζ}1ββ2}?ddβ{ˉΖ}=[B]{ˉΖ}(7)上式中,{Ζˉ}=[{Ζ}Τ1ββ2]Τ是{Z}的擴(kuò)容向量.經(jīng)擴(kuò)容處理后,(7)式已變成了齊次方程,它的解用精細(xì)積分法極易求得.如前一樣,此時(shí)矩陣[B]也可取為常數(shù)矩陣.根據(jù)矩陣?yán)碚?一階齊次微分方程組(7)的解可寫成{{Ζ}1ββ2}β2=e∫β1β2[B]dβ{{Ζ}1ββ2}β1=e[B](β2-β1){{Ζ}1ββ2}β1=[Τ]{{Ζ}1ββ2}β1矩陣[T]稱為傳遞矩陣,它表示了相鄰截面β2和β1處兩個(gè)狀態(tài)向量之間的關(guān)系,其值用精細(xì)積分法可精確求得.2任何形狀的二維孔和彈簧波2.1接觸問題的有限元離散解析對(duì)于任意形狀二維孔穴,在諧激勵(lì)作用下的外聲場(chǎng),輻射聲壓滿足下列Helmholtz方程?2pf+ω2c02pf=[?2?r2+1r??r+1r2?2?θ2+k2]pf=0(8)無窮遠(yuǎn)處滿足Sommerfeld輻射條件limr→∞r(nóng)1/2(?pf?r+ikpf)=0(9)式中,c0為流體介質(zhì)中的聲速,k=ω/c0為波數(shù).設(shè)在二維孔穴的內(nèi)域中設(shè)置一封閉曲線Γ′(虛擬邊界)如圖3所示,并在其上布置連續(xù)型源強(qiáng).根據(jù)疊加原理,滿足方程(8)和無窮遠(yuǎn)條件(9)的解可表示為以下封閉曲線積分的形式pf(Q)=∮Γ?′σ(Ρ)G(Q,Ρ)d?！?Ρ)(10)式中σ(P)為源強(qiáng)配置曲線Γ′上未知的源強(qiáng)密度函數(shù)(單位長(zhǎng)度的源強(qiáng)),G(Q,P)為二維自由空間的Green函數(shù),其表達(dá)式為G(Q?Ρ)=-i4Η0(2)(kr)(11)其中,r=ΡQˉ表示源點(diǎn)P和場(chǎng)點(diǎn)Q之間的距離,H0(2)(kr)為0階第二類Hankel函數(shù).因此,只要能求得虛擬配置曲線?！渖系脑磸?qiáng)密度函數(shù)σ(P),使得由(10)式描述的聲壓場(chǎng)在真實(shí)振動(dòng)邊界Γ上滿足給定的條件,那么根據(jù)解的唯一性定理可知,(10)式所描述的聲壓場(chǎng)即為原問題的解.然而,對(duì)于任意形狀的振動(dòng)孔穴,要在其振動(dòng)邊界上精確滿足給定條件是很困難的.在文中,采用在一條與真實(shí)振動(dòng)邊界相似的虛擬源強(qiáng)配置曲線Γ′上用單元離散插值(虛擬邊界元法),并結(jié)合最小二乘法逼近真實(shí)邊界條件,求解了二維孔穴的振動(dòng)聲輻射問題.但上述方法存在兩點(diǎn)不足:1)由于受離散單元低階形函數(shù)的限制,計(jì)算結(jié)果僅適用于中、低頻段;2)計(jì)算時(shí)間長(zhǎng),效率較低.為克服以上不足,本文提出采用多圓形虛擬源強(qiáng)配置方案,并在每一條圓形配置曲線上將源強(qiáng)密度函數(shù)σ(P)用Fourier級(jí)數(shù)展開,這時(shí)求未知函數(shù)σ(P)就化為求級(jí)數(shù)展開系數(shù),如果再結(jié)合快速Fourier逆變換法,一種高精度、高效率求解任意形狀二維孔穴Helmholtz外問題的快速算法就能形成.2.2快速furhen逆變換idft算法對(duì)于某些邊界Γ,虛擬源強(qiáng)配置區(qū)域有時(shí)可選多個(gè)圓形曲線,曲線間可以相互重疊和交叉(見圖4).在每一條圓形配置曲線上將源強(qiáng)密度函數(shù)σ(P)展開成Fourier級(jí)數(shù)σs(ηs)=∑n=-∞+∞cnseinηs?0≤ηs≤2π?s=1,2,?,Νc(12)式中,Nc為圓形配置曲線的個(gè)數(shù)(圖4中Nc=3),將上式代入(10)式,得聲壓表達(dá)式為pf(Q)=∑s=1Νc∑n=-∞+∞cns∫02πG(Q,ηs)einηsdηs(13)在積分∫2π0G(Q,ηs)einηsdηs中,由于含有周期振蕩因子einηs,因而當(dāng)級(jí)數(shù)展開項(xiàng)數(shù)n較大時(shí),采用一般的數(shù)值積分法(如Gauss積分法),會(huì)產(chǎn)生很大的計(jì)算誤差,使計(jì)算結(jié)果嚴(yán)重失真.如采用處理振蕩積分的Filon積分法,雖然能獲得較高的計(jì)算精度,但計(jì)算時(shí)間較長(zhǎng),不適合于結(jié)構(gòu)-聲耦合問題的計(jì)算.為此,采用快速Fourier逆變換法計(jì)算上述積分.數(shù)值Fourier逆變換(IDFT)的定義x(n)=1Ν(∑k=0Ν-1X(k)ei(2πΝ)nk)?n=0,1,?,Ν-1(14)其中,X(k)(k=0,1,2,N-1)為N個(gè)樣本值.將(13)式中的積分區(qū)間(0,2π)進(jìn)行N(N=2M)等分,即,Δηs=2π/N,ηs=kΔηs(k=0,1,…,N-1).根據(jù)梯形求積分公式有∫02πG(Q,ηs)einηsdηs=∑k=0Ν-1G(Q,kΔηs)einkΔηsΔηs=∑k=0Ν-1G(Q,kΔηs)ei(2πΝ)nk(2πΝ)=2π(1Ν∑k=0Ν-1G(Q,kΔηs)ei(2πΝ)nk)(15)從上式與數(shù)值Fourier逆變換的定義式(14)進(jìn)行比較可以看出:(13)式振蕩函數(shù)的數(shù)值積分恰恰能化為數(shù)值Fourier逆變換的格式加以計(jì)算,因而可以借助快速Fourier逆變換(IFFT)程序高效、精確地求得其積分值.根據(jù)抽樣定理:為避免頻率混疊(影響計(jì)算精度),(12)式中的Fourier級(jí)數(shù)展開項(xiàng)數(shù)n和函數(shù)G(Q,kΔηs)(k=0,1,…,N-1)在圓形配置曲線上樣本數(shù)N必需滿足:n≤N/2.在(13)式中,由于采用了多圓形虛擬源強(qiáng)配置曲線(當(dāng)真實(shí)邊界Γ成扁平形尤為必要),因而能有效地避免系統(tǒng)矩陣的病態(tài),有利于提高計(jì)算精度.當(dāng)孔穴邊界與圓形曲線相差不大時(shí),采用1～2個(gè)虛擬源強(qiáng)配置曲線(如圖3)即可獲得較滿意的計(jì)算精度.在同樣的計(jì)算條件下求解任意形狀空穴在已知法向振速激勵(lì)下的輻射外聲場(chǎng),本文方法的計(jì)算速度比一般虛擬邊界元法快近50倍.3彈性環(huán)法向位移計(jì)算在每條虛擬源強(qiáng)配置曲線上,取虛擬源強(qiáng)函數(shù)的Fourier級(jí)數(shù)求和項(xiàng)數(shù)為-Ms～+Ms.對(duì)于彈性環(huán)-聲耦合問題,若將表面聲壓(13)式代入彈性環(huán)運(yùn)動(dòng)方程(4)得ddβ{vˉwˉΦˉΜˉ-QˉΝˉ}=1Ra[A]{vˉwˉΦˉΜˉ-QˉΝˉ}+1Κ{0000f(β)0}-1Κ∑s=1Νc∑n=-Μs+Μscns{0000∫02πG(β,ηs)einηsdηs0}(16)利用前面第一節(jié)所述的齊次擴(kuò)容精細(xì)算法,分別計(jì)算彈性環(huán)在已知外激勵(lì)力f(β)和廣義力∫2π0G(β,ηs)einηsdηs作用下的法向位移,記為wˉf(β)、wˉnsp(β).根據(jù)疊加原理,在外激勵(lì)力和聲壓共同作用下彈性環(huán)的法向位移可寫為w(β)=Rawˉ(β)=RaΚ(wˉf(β)-∑s=1Νc∑n=-Μs+Μscnswˉnsp(β))(17)由彈性環(huán)和流體介質(zhì)的相容條件,有?pf(β)?nβ=ρ0ω2w(β)(18)式中,ρ0為流體介質(zhì)質(zhì)量密度,nβ為彈性環(huán)的外法線方向.將方程(13)和(17)代入方程(18),并整理成如下求解形式∑s=1Νc∑n=-Μs+Μscns(∫02π?G(β,ηs)?nβeinηsdηs+ρ0Raω2Κwˉnsp(β))=ρ0Raω2Κwˉf(β),β∈Γ(19)通過在彈性環(huán)上配點(diǎn)βi(i=1,2,…,M),其中,M≥Nc(2Ms+1),對(duì)方程(19)運(yùn)用最小二乘法,可得下列線性方程組[U]{c}={q}(20)一旦從上式中求得每條虛擬源強(qiáng)配置曲線上的Fourier級(jí)數(shù)展開系數(shù)cns({c})后,由(13)式就可計(jì)算任意點(diǎn)的聲壓,再?gòu)?17)式可計(jì)算彈性環(huán)上任意點(diǎn)的法向位移.另外,對(duì)(19)式的振蕩函數(shù)積分,只需將(15)式中的函數(shù)G(Q,ηs)換為函數(shù)?G(Q,ηs)?nβ按相同步驟進(jìn)行.4m的輻射聲壓特性算例1為驗(yàn)證本文方法的有效性,現(xiàn)分析一載流橢圓環(huán)在單位幅值f0=1N集中諧激勵(lì)力作用下的聲輻射問題(如圖5),記橢圓環(huán)長(zhǎng)半軸與短半軸之比為m=a/b,其幾何和物理參數(shù)如下:E=210GPa,ρ=7.84×103kg/m3,ρ0=1.0×103kg/m3,c0=1500m/s,Ra=1m,A=5cm2,I=10.42cm4,參考聲壓為:p0=2.0×10-5Pa.圖6、圖7和圖8中比較了圓環(huán)時(shí)(m=1)表面聲壓解析解與本文解的計(jì)算結(jié)果,在算例中采用單圓形虛擬源強(qiáng)配置,虛擬圓上樣本數(shù)取:N=2M=28.從圖中可看出本文方法的計(jì)算結(jié)果與解析解結(jié)果相當(dāng)吻合,說明本文的方法是有效的.圖9、圖10分別給出了橢圓環(huán)在觀察點(diǎn)(r,θ)=(5,π/2)和(5,π)處的輻射聲壓級(jí)(SPL)隨外激勵(lì)頻率的變化情況,從圖10中可看到:在外激勵(lì)頻率低于350Hz時(shí),輻射聲壓隨橢圓度m的增加而變大.而當(dāng)激勵(lì)頻率高于350Hz時(shí),輻射聲壓幾乎與m無關(guān),而且共振峰值對(duì)應(yīng)的激勵(lì)頻率也與圓環(huán)時(shí)相近.表1給出了相似虛擬邊界元法、單圓形配置法和三圓形配置法求解m=2的橢圓環(huán)在其長(zhǎng)軸作用有單位幅值集中諧激勵(lì)力的輻射聲壓比較.三種方法的相關(guān)參數(shù)如下:(1)在相似虛擬邊界元法中,虛擬橢圓邊界的短軸取b′=0.75b,用二次單元將其離散為40個(gè)單元.(2)在Fourier級(jí)數(shù)方法的單圓形配置方案中,虛擬源強(qiáng)配置圓的圓心取(x,y)=(0,0),半徑R′=0.75b,虛擬圓上Fourier級(jí)數(shù)展開項(xiàng)數(shù):-40～+40,樣本數(shù)取:N=2M=28.(3)在Fourier級(jí)數(shù)方法的三圓形配置方案中,三條虛擬源強(qiáng)配置圓的圓心位置和半徑分別取(x1,y1)=(-0.5,0)、R′1=0.65b;(x2,y2)=(0,0)、R′2=0.75b;(x3,y3)=(0.5,0)、R′3=0.65b,每條虛擬圓上Fourier級(jí)數(shù)展開項(xiàng)數(shù):-30～+30,樣本數(shù)取:N=2M=28.從表1中可以