一維波在水下多層介質中傳播的時域解析解

金澤宇, 殷彩玉, 諶 勇, 黃修長, 華宏星

(上海交通大學 振動、沖擊、噪聲研究所,上海 200240;上海交通大學 機械系統與振動國家重點實驗室,上海 200240)

隨著現代武器的發展和對艦艇生命力要求的提高，對艦艇低輻射噪聲和高抗爆炸沖擊性能的研究越來越受到各國海軍的重視。目前通常在艦艇的艇體表面敷設一層多孔介質覆蓋層以達到抗沖隔聲的目的。因此研究在水下爆炸沖擊波輸入情況下，表面敷設覆蓋層后對艦艇艇體表面壓力波的影響顯得極為重要。

波在多層介質傳播的問題有過很多研究。書籍[1-2]介紹了許多計算簡諧問題的公式，得到了解析和半解析的解。然而這些解不能解決考慮到全部復雜反射和透射情況的波在多層介質中傳播的問題。在本文中邊值問題是主要考慮的對象，這個問題的解通常采用頻域方法或者時域方法。在頻域(傅里葉對時間做變換)方法中，波有限元法(WFE method)[3-7]是很好的選擇，雖然這些方法會產生數值困難，但可以通過使用一些代數技巧避免。

本文的重點是時域方法，它主要包含了兩種方法，一種是數值方法，一種是解析方法。前人做過了許多工作，例如有限時域差分法[8-11]，有限元法(FEM)[12-15]。Desceliers等[16]發展了一種新的數值雜交方法來仿真在給定瞬態載荷的情況下瞬態波在有多層的可以是液體或者固體的半無限介質中傳播的問題。對于有限元法，Brasek[17-18]得到了一維有覆蓋層的鋁結構在階躍載荷下的瞬態響應，并且研究了在指數波輸入情況下的金屬結構的響應特性。他的研究顯示不管覆蓋層的類型和邊界介質如何，結構的應力和節點速度隨著覆蓋層剛度的下降而上升。Neilson等[19]運用有限元和雙漸進逼近法(DAA)得到了在平面階躍壓力下沉浸雙層彈性球殼的瞬態響應，兩層殼之間充滿流體。Bergersen[20]得到了在水下爆炸沖擊波下的貼有橡膠材料的金屬圓柱的動態響應。他的研究發現在一定的幾何條件和覆蓋層材料特性條件下，貼覆蓋層在動態響應方面有負面的影響。雖然在理論上時域數值方法能解決復雜的問題，但是一些缺點仍然存在。一是數值方法本身不精確，更重要的是我們可能不會完全理解問題的機理，而且在對參數分析時，需要大量的工作。

對于解析解，雖然它們不能獲得所有問題特別是復雜問題的精確解，但是它們對于相對簡單的問題的結果能夠幫助人們透徹地理解機理。此外，還能夠為復雜問題提供指導。在過去的研究中，一些人做了解析方面的工作。Taylor[21]研究了弱沖擊與氣背板的作用，得到了一維解析解。Snay等[22]運用特征線法考慮了強沖擊和氣背板作用，改進了Taylor的結果。Schechter等[23]考慮的與弱沖擊作用的板可以是氣背或者水背的。Huang[24]運用Laplace和Hankel變換研究了在球面沖擊波作用下的大彈性平板的初始瞬態響應。童宗鵬[25]采用波動理論，拓展Taylor平板模型，研究了水下爆炸沖擊波在水、覆蓋層、空氣層、船體鋼中的傳播過程，建立了實心覆蓋層結構和含空腔覆蓋層結構的流固耦合模型，但其只考慮了一次反射波。諶勇等[26]運用多自由度系統模擬艦船覆蓋層的動態模型。然而這些解存在一些假設，并且沒有給出考慮多重反射透射的解析解。

本文的目的是基于波動理論提出一種解決一維波在多層介質中傳播的時域解析解。這種新方法通過考慮在多層介質中任意輸入下的多次反射透射獲得準確的瞬態解。通過無量綱參數分析獲得峰值壓力隨覆蓋層阻抗、厚度和指數波衰減常數的變化。

1 時域解析解

水中壓力波經覆蓋層-鋼板后，會產生一系列的反射和透射現象，隨著反射和透射次數的增多，情況越來越復雜。本文所提出的方法不考慮壓力波在覆蓋層-鋼板中反射透射的具體傳播方式，而是從某一特定時刻出發，分析在該時刻壓力波在覆蓋層和鋼板中可能的傳播方式，將這些傳播方式組合即可得到壓力波在多層介質中的響應。具體的思路如下：對于某一特定時刻tm,n，從覆蓋層透射到水中的壓力增量記為f(tm,n,t)，該透射波在覆蓋層中傳播m次，在鋼板中傳播n次。對這一特定時刻，m和n是定值，但從整體的時間上看，m可以取1→∞，n可以取0→∞，m和n的任意組合都有可能，對于這一特定時刻，就是所有組合中的一種?；谏鲜鱿敕?，若波在覆蓋層以及在鋼板中傳播一次的時間已知，則可計算出所有m和n的組合所對應的時刻，對這些時刻從小到大排序，就得到每次從覆蓋層中透射到水中的壓力波的情況，獲得波在多層介質中傳播的時域解析解。

1.1 基本假設

在壓力波波幅不大的情況下，將水、覆蓋層、鋼板簡化為一維問題，并做以下假設：

(1) 假設應力波在各介質中的傳播形式為縱波；

(2) 假設各個介質為均勻介質，有固定線性阻抗；

(3) 忽略波傳播過程中阻尼的影響；

(4) 忽略覆蓋層、鋼板的剛體運動；

(5) 忽略波傳播引起的介質變形對波在每層介質中傳播時間的影響。

1.2 系統簡化模型及模型參數

對于有水、覆蓋層和鋼板三層介質的系統，簡化后的模型為：介質1是水，且左邊是無限的；中間介質2是有限厚度的覆蓋層；右側介質3是鋼板，鋼板右側邊界Ⅲ可以是不同的邊界條件，如氣背、水背、剛性邊界等。系統簡化模型圖如圖1所示。

簡化后的系統的物理參數如下：密度：ρ1，ρ2，ρ3；波速：c1，c2，c3；阻抗率：R1，R2，R3；介質厚度：h2，h3。

波在不同介質中傳播參數如下：反射系數與透射系數：r12，s12，r21，s21，r23，s23，r32，s32，rⅢ(其中rab為波由介質a到介質b的反射系數，sab為波由介質a到介質b的透射系數，rⅢ為邊界Ⅲ的反射系數)。

圖1 系統簡化模型

1.3 時域解析解

根據聲學理論知，均勻理想彈性媒質中一維小振幅線性聲波方程為:

(1)

由之前的假設可知，本文考慮的波滿足此方程。在邊界Ⅱ處反射波返回到邊界Ⅰ之前，介質1中存在行波解：

(2)

介質2中的透射波：

(3)

粒子振動速度滿足歐拉方程：

(4)

可得：

(5)

由此可知：對于任意F、G，壓力p與粒子振動速度v的關系與簡諧入射波是等同的。所以輸入可取簡諧波、指數衰減波、三角波等任意波。

(1) 特定時刻濕表面壓力增量

設壓力波在介質2和介質3傳播一次的時間分別為：

(6)

對于某一特定時刻，有一壓力波從覆蓋層透射到水中，該壓力波在覆蓋層中傳播m次，在鋼板中傳播n次，則對應的時刻tm,n為：

tm,n=m×ts2+n×ts3

(7)

設在特定時刻tm,n由介質2透射到介質1中的壓力波增量為f(tm,n,t)，pin(t)為初始時刻入射波，f(t0,0,t)為初始時刻的反射波。其中初始時刻反射波可以表示為：

f(t0,0,t)=r12·pin(t)

(8)

在初始時刻以后的時間里(即m≥1)，當n=0時，波沒有透射到介質3中，此時從介質2透射到介質1中的波記為f(tm,0,t)：

f(tm,0,t)=s12s21(r21r23)m-1r23·pin(t-tm,0)

(9)

其過程可以簡述為入射波透射進入介質2，在介質2中來回傳播m-1次，最終透射回到介質1中。

當n≠0時，此時波從介質2透射到介質3的次數從1→k都有可能。

波從介質2到介質3透射1次

f1(tm,n,t)=

pin(t-tm,n)

(10)

波從介質2到介質3透射2次

f2(tm,n,t)=

pin(t-tm,n)

(11)

更一般地，波從介質2到介質3透射k次：

fk(tm,n,t)=

pin(t-tm,n)

(12)

綜上所述：當m≥1,n≠0時，tm,n時刻由介質2透射到介質1的波為上述k種可能性的疊加，即：

(13)

當m≥1,n=0時，tm,n時刻由介質2透射到介質1的波由方程(9)給出結果，當m=0，n=0，即初始時刻時的反射波由方程(8)給出結果。

(2) 每一時刻濕表面的壓力

(1)中介紹了對于特定時刻tm,n從介質2透射到介質1中的壓力波增量的計算方法，即m、n所有組合中的一種。現在算出m、n所有組合所對應的時刻tm,n和該時刻從介質2透射到介質1的壓力增量f(tm,n,t)。先對tm,n進行升序排列，根據tm,n的升序排列，將對應該時刻的f(tm,n,t)重新排列，就得到順序的每一時刻從介質2透射到介質1的壓力增量。而對于每一時刻濕表面壓力，則為該時刻及以前所有順序的時刻從介質2透射到介質1的壓力增量的和。

(14)

2 數值驗證和算例

為了驗證時域解析解的正確性，將圖1中各介質材料參數賦值如表1所示。取入射波為指數衰減波p0exp(-t/θ)，取衰減時間常數θ= 0.6 ms，壓力峰值p0為1 Pa。波在水、覆蓋層和鋼板中傳播，邊界Ⅲ可以是剛性邊界、氣背邊界和水背邊界。有限元模型由8節點的線性聲學的方形六面體單元構成。在ABAQUS中的幾何模型如圖2所示，覆蓋層和鋼板包括740個線性聲學單元，水和空氣各包括50個單元，整個系統包括3 376個節點。

表1 三種介質的材料參數

圖2 一維系統在ABAQUS中的模型

圖3、圖4和圖5分別為濕表面壓力時域解析解和有限元解在剛性邊界、氣背邊界和水背邊界下的對比。圖中有限元解開始時有一個下降，而解析解沒有，是因為有限元解采樣間隔小，在開始時采樣點比解析解多導致的。開始時刻有限元解的下降是由于指數衰減波的衰減作用造成的，而不是有透射波到水中造成的。

從圖中可以看出解析解與有限元解吻合得很好，解析解每次突變的時間和壓力與有限元解對比都相差不大。不大的壓力差是由于有限元存在數值阻尼和誤差造成的。解析解與有限元解的對比可以說明解析解的正確性。

圖6為指數衰減波對于不同邊界濕表面壓力的變化，從圖中可以看出：① 三種邊界下，濕表面總壓力最終都趨于0；② 衰減速度：氣背最快、水背次之、剛性邊界最慢；③ 氣背邊界時，總壓力先衰減到負壓再趨于0。如果壓力充分大，氣背產生的負壓大于靜水壓，板將和水脫離。水背和剛性邊界從正壓直接衰減為0。

圖7為指數衰減波對于三種邊界濕表面振動速度的變化，振動速度可以由壓力通過方程(4)求得，從圖中可以看出剛性邊界振動速度基本在0附近變化，而水背邊界、氣背邊界振動速度則有較大變化。這是由于邊界Ⅲ的約束能力影響的，如果邊界Ⅲ自由，約束能力弱，振動速度有較大變化。而水背、和剛性邊界對邊界Ⅲ約束能力強，使得振動速度不易有大變化。

圖3 剛性邊界濕表面壓力解析解與有限元解對比曲線

圖6 指數衰減波輸入對于三種邊界濕表面壓力曲線

圖8為指數衰減波濕表面位移的變化，位移曲線是由振動速度曲線積分得到。從圖中可以看出：剛性邊界時，濕表面位移為0；水背和氣背邊界時，濕表面有剛體位移。因此，如果壓力波過大，則計算結果將不再準確。

3 指數衰減波輸入的近似結果比較

上面提出的解析方法可考慮所有反射及透射的影響，得到的濕表面壓力的精確結果。下面針對邊界條件為剛性邊界、氣背邊界以及水背邊界時，將精確結果與僅考慮前幾次反射透射的近似結果進行比較。

圖9、圖10、圖11分別為剛性邊界、氣背邊界、水背邊界的近似結果與解析結果的對比圖。從這3個圖中可以看出：

(1) 僅考慮濕表面第一次反射波而不考慮在介質2與3中波的情況，所得近似結果與實際結果相差很大；

(2) 僅考慮波在邊界Ⅱ上反射一次，不考慮波在介質3中的傳播情況，對于剛性邊界近似結果與實際結果相差不大，但是近似結果無法反映壓力由于多次反射、透射所造成的震蕩效應；對于氣背邊界和水背邊界，結果則相差很大，考慮多次反射、透射的結果衰減的速度遠大于僅考慮波在邊界Ⅱ上反射一次的結果。

所以僅考慮幾次反射、透射而認為后續反射波、透射波很小可以忽略的近似方法存在缺陷。想得到準確的結果，需要考慮足夠多的反射、透射波的作用。

圖9 剛性邊界下近似結果與解析結果比較

4 指數衰減波的無量綱化參數分析

對指數衰減波進行無量綱化，設td2=ts2/θ，td3=ts3/θ，n1=R1/R2，n2=R2/R3，四個變量加上r3五個變量，構成無量綱化后的5個參數。以下研究各個參數變化對結果的影響。

4.1 覆蓋層阻抗的影響

圖12為濕表面壓力峰值隨n1/n2比值變化曲線。從圖中可以看出：隨著覆蓋層阻抗減小(即n1/n2增大)，各個邊界壓力波峰值先減小后增大再減小。當n1/n2的比值近似為0.06時，剛性邊界的峰值出現最小值；比值近似為0.5時，氣背邊界、水背邊界的峰值出現最小值；當n1/n2的比值近似為3時，剛性邊界峰值出現最大值，超過了覆蓋層為剛性層時的峰值；比值近似為5時，氣背邊界、水背邊界峰值出現最大值，但小于覆蓋層為剛性層時的峰值。當n1/n2的比值繼續增大，三種邊界壓力波的峰值出現明顯減小。

峰值曲線先減小后增大再減小的原因可分析如下：① 對于剛性邊界，壓力的峰值主要是由第一次反射波、第二次透射波、第三次透射波(第三次透射是入射波進入介質2再進入介質3在邊界Ⅲ上反射后透射到介質2再透射回介質1的波)的疊加結果決定。在阻抗較大時，起決定作用的是第一次反射波，隨著阻抗的減小，第一次反射波幅值減小，峰值曲線下降。阻抗繼續減小達到最低點，此時起決定作用的是第二次和第三次透射波，隨著阻抗減小，第二次和第三次透射波幅值增強，峰值曲線上升到最大值。隨著阻抗繼續減小，覆蓋層的特性趨于空氣，第一次透射波趨于0，后續作用越來越小，而第一次入射波和反射波之和也趨于0，所以最終結果使峰值越來越?。虎?對于氣背邊界和水背邊界的情況，壓力峰值主要是由第一次反射波、第二次透射波的疊加結果決定，此處與剛性邊界不同，由于第三次透射波對于邊界Ⅲ剛性，壓力波變為原來2倍，而若是氣背、水背邊界反射波相位相反，疊加將減小壓力波峰值，其余討論與剛性邊界相同。

4.2 覆蓋層厚度與衰減常數的影響

(1)改變td2、td3使之同比增大或減小

td2、td3同比增大或減小相當于介質2、3的厚度同比增加，或者波速同比減小，亦或θ減小。

圖13為同比例改變td2、td3，濕表面壓力峰值曲線，從圖中可以看出：隨著td2、td3同比例增大，壓力波峰值減??；衰減的速度越快。即同比增加介質2、3的厚度或者減小輸入波的衰減常數，可使濕表面上的總壓力峰值減小，衰減更快。

圖12 改變n1/n2濕表面無因次壓力峰值曲線

(2) 固定td3，改變td2

固定td3，改變td2中的參數，使td2，td3具有不同的比例。這樣做相當于改變覆蓋層厚度。

圖14為固定td3，改變td2中的參數濕表面壓力峰值的變化的曲線。從圖中可以看出對于剛性邊界隨著td2增大，壓力波峰值先減小再突變增大再減小；對于氣背和水背邊界隨著td2增大時，壓力波峰值減小。對于剛性邊界，壓力波峰值是由衰減作用和第三次透射作用共同決定的，曲線的前半段后后半段的下降趨勢是由于反射波與入射波疊加時的時間延長，入射波的衰減作用導致的。也可以理解為覆蓋層厚度變大，儲存能量的能力大了，等到釋放能量的時候，入射波已經衰減很多了，導致疊加后的波峰減小，而峰值曲線的突變增大則是由第三次透射作用導致的，厚度小時壓力波在介質2中傳播多次才有一次從介質3中傳播回來的波，但第一次從介質3傳播回來的波壓力峰值最大，而且由于邊界Ⅲ剛性，傳播回來的波是原壓力的二倍，隨著介質2的厚度增大，當這個波成為第三個傳播回介質1的波時，壓力峰值就會出現突變。對于氣背和水背邊界第一次從介質3傳播回來的波相位相反，疊加作用使波的峰值減小，不會出現壓力波峰值突變的作用，峰值隨厚度減小的原理同剛性邊界原理。

5 結 論

本文基于波動理論對一維波在覆蓋層-鋼板中的傳播情況進行了理論推導，得到了波在3層介質中傳播的一維時域解析解，計算了輸入為指數衰減波時濕表面壓力的時域解析結果。同時,針對輸入為指數衰減波的情況，該方法的結果通過與有限元的結果的對比證明了該方法的正確性。將該方法的結果與只計一、兩次反射透射影響的近似結果進行了對比，結果發現，對于剛性邊界條件，精確結果與近似結果吻合較好；對于氣背邊界和水背邊界，近似結果誤差較大。對指數衰減波進行無量綱化，分析了改變覆蓋層的厚度、特性阻抗以及指數衰減波的衰減時間常數對結果的影響。分析結果表明，當覆蓋層的特性阻抗減小時，壓力波峰值先減小后增大再減??；覆蓋層的厚度增大，會使壓力波峰值減小(對于剛性邊界壓力波峰值先減小后突變增大再減小)，衰減速度快，會使壓力波峰值減小。

[1] Graff K F. Wave motion in elastic solids [M]. Dover Publications Inc., New York, 1975.

[2] Cremer L, Heckel M, Petersson B A T. Structure-borne sound [M]. Third ed., Springer, Berlin, 2005.

[3] Mace B R, Duhamel D, Brennan M J, et al. Finite element prediction of wave motion instructural waveguides [J]. Journal of the Acoustical Society of America,2005, 117(5): 2835-2843.

[4] Duhamel D, Mace B R, Brennan M J. Finite element analysis of the vibrations of waveguides and periodic structures[J]. Journal of Sound and Vibration,2006, 294 (1-2):205-220.

[5] Mace B J, Manconi E. Modelling wave propagation in two-dimensional structures using finite element analysis [J]. Journal of Sound and Vibration,2008, 318(4-5):884-902.

[6] Renno J M, Mace B R. On the forced response of wave guides using the wave and finite element method [J]. Journal of Sound and Vibration,2010, 329(26):5474-5488.

[7] Renno J M, Mace B R. Calculating the forced response of two-dimensional homogeneous media using the wave and finite element method [J]. Journal of Sound and Vibration,2011, 330:5913-5927.

[8] Akyurtlu A, Werner D H. BI-FDTD: a novel finite-difference time-domain formulation for modeling wave propagation in bi-isotropic media [J]. IEEE Trans. Antenn. Propag,2004, 52(2):416-425.

[9] Gandhi O P, Gao B Q, Chen J Y. A frequency-dependent finite difference time-domain formulation for general dispersive media [J]. IEEE Trans. Microw. Theory Tech,1993, 41:658-664.

[10] Luebbers R, Hunsberger F P, Kunz K, et al. A frequency-dependent finite-difference time-domain formulation for dispersive materials [J]. IEEE Trans. Electromagn. Comp,1990, 32:222-227.

[11] Taflove A. Review of the formulation and applications of the finite-difference time-domain method for numerical modeling of electromagnetic wave interactions with arbitrary structures [J]. Wave Motion,1988,1(6):547-582.

[12] Bathe K J, Wilson E L. Numerical methods in finite element analysis[M]. Prentice-Hall, New York, 1976.

[13] Givoli D. Recent advances in the DtN FE method[J]. Arch. Comput. Methods Eng.,1999,6(2):71-116.

[14] Givoli D. Nonreflecting boundary-conditions [J]. J. Comput. Phys,1991,94(1):1-29.

[15] Zienkiewicz O C, Taylor R L. The finite element method [J]. McGraw-Hill, NY, 1991.

[16] Desceliers C, Soize C, Grimal Q, et al. A time-domain method to solve transient elastic wave propagation in a multilayer medium with a hybrid spectral-finite element space approximation[J]. Wave Motion,2008,45(4):383-399.

[17] Brasek T P. Effect of surface coating on one-dimensional system subjected to unit step pressure wave [D]. California: Naval postgraduate School, NPS-ME-94-006, 1994.

[18] Brasek T P. Response of dual-layered structures subjected to shock pressure wave [D]. California: Naval Postgraduate School, 1994.

[19] Neilson H C, Everstine G C, Wang Y E. Transient response of a submerged fluid-coupled double-walled shell structure to a pressure pulse [J]. Journal of the Acoustic Society of America, 1981,70:1776-1782.

[20] Bergersen J K L. Effect of surface coating on cylinders subjected to underwater shock [D]. Naval postgraduate School, Monterey, California, 1992.

[21] Tayor G I. The pressure and impulse of submarine explosion waves on plates [J]. In: Underwater Explosion Research, Vol. I, Office of Naval Research, 1950:1155-1173.

[22] Snay H G, Christian E A. The response of air-backed plates to high-amplitude underwater shockwaves[R]. NAVORD, Report 2462,1952.

[23] Schechter R S, Bort R L. The response of two fluid-coupled plates to an incident pressure pulse [R]. Naval Research Laboratory, Memorandum report 4647, 1981.

[24] Huang H. A qualitative appraisal of the doubly asymptotic approximation for transient analysis of submerged structures excited by weak shock waves [R]. AD-A015941, 1975.

[25] 童宗鵬.艦船抗沖瓦結構流固耦合建模及水下抗爆機理研究[D].上海：上海交通大學，2006.

[26] 諶勇,張志誼,華宏星,等.彈性泡沫夾芯結構的水下爆炸響應分析[J].振動與沖擊，2009, 28(11)：25-29.

CHEN Yong,ZHANG Zhi-yi,HUA Hong-xing, et al. Analytic model for the response to water blast of one-dimensional naval structure with elastic foam coatings [J]. Journal of Vibration and Shock, 2009, 28(11)：25-29.