• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    非均勻波導(dǎo)中的聲聚焦

    2020-04-30 08:33:20郭威楊德森1
    物理學(xué)報 2020年7期
    關(guān)鍵詞:聚焦點散射體入射波

    郭威 楊德森1)2)?

    1) (哈爾濱工程大學(xué), 水聲技術(shù)重點實驗室, 哈爾濱 150001)

    2) (海洋信息獲取與安全工業(yè)和信息化部重點實驗室(哈爾濱工程大學(xué)), 哈爾濱 150001)

    3) (哈爾濱工程大學(xué), 水聲工程學(xué)院, 哈爾濱 150001)

    理論研究了聲波在非均勻波導(dǎo)中的空間聚焦問題, 利用多模態(tài)導(dǎo)納法構(gòu)建波導(dǎo)內(nèi)任意位置處聲壓與入射聲壓在模態(tài)域的映射關(guān)系, 計算使聲波聚焦于空間某位置時的最佳入射波, 并畫出了相應(yīng)的聚焦聲場.研究了三種非均勻情況: 水平變截面波導(dǎo)、含散射體波導(dǎo)以及聲速垂直變化波導(dǎo).結(jié)果表明, 當(dāng)輸入最佳入射波時, 在非均勻波導(dǎo)中可以產(chǎn)生良好的單點或多點聲聚焦效果, 聲波的聚焦過程充分地利用了波導(dǎo)結(jié)構(gòu)及介質(zhì)非均勻性對聲波的散射作用.

    1 引 言

    聲波聚焦在聲通信、聲成像、無損檢測以及超聲醫(yī)學(xué)治療等領(lǐng)域具有重要的應(yīng)用價值.目前實現(xiàn)聲波聚焦的一種手段是設(shè)計特殊的材料和結(jié)構(gòu), 例如二維聲子晶體[1]、聲透鏡[2,3]、負(fù)折射聲超構(gòu)材料[4]等, 或借助熱聲學(xué)手段[5,6], 改變傳播介質(zhì)的折射率, 控制聲波的傳播路徑, 實現(xiàn)空間聚焦.該方法已在均勻自由場中產(chǎn)生了良好的聲聚焦效應(yīng).然而, 對于聲波在非均勻波導(dǎo)中的傳播, 由于介質(zhì)本身的不均勻性及邊界的散射作用, 傳播過程較自由場中復(fù)雜, 控制聲波的傳播路徑也相對困難, 較難應(yīng)用該手段實現(xiàn)非均勻波導(dǎo)中的聲聚焦.另一種研究波聚焦的方法為相位共軛法(phase conjugation),首先應(yīng)用于光學(xué)領(lǐng)域[7], 以實現(xiàn)光波穿過復(fù)雜散射區(qū)域后的空間聚焦.相位共軛法是基于絕熱物理系統(tǒng)的時間反轉(zhuǎn)不變性, 對點光源穿過散射區(qū)域后產(chǎn)生的透射波進行相位共軛(時間反轉(zhuǎn))處理, 再將處理后的光波作為入射波, 反向穿過散射區(qū)域, 最終在原本的點光源位置處形成聚焦[8].此外, 很多學(xué)者將相位共軛法與透射矩陣[9,10]相結(jié)合, 實驗測量透射矩陣, 進而給出入射波, 使其穿過散射區(qū)域后直接產(chǎn)生空間聚焦[11,12]或時間空間同時聚焦[13?16].Fink等[17]將相位共軛法引申到聲學(xué)領(lǐng)域中, 稱為時間反轉(zhuǎn)法 (time reversal), 可實現(xiàn)聲波穿過非均勻介質(zhì)后的空間聚焦[18], 其優(yōu)點是無需介質(zhì)的先驗信息即可完成自適應(yīng)聚焦[19].但相位共軛法多為研究波穿過散射區(qū)域后的聚焦[8?18], 且依賴于實驗手段, 相對缺乏系統(tǒng)的理論分析, 而且對在散射區(qū)域內(nèi)部實現(xiàn)聚焦的相關(guān)研究也較少.

    鑒于上述問題, 本文從均勻剛硬波導(dǎo)入手, 進而研究水平變截面、含散射體以及聲速垂直變化三種非均勻波導(dǎo)內(nèi)的聲聚焦問題.利用一種求解非均勻波導(dǎo)內(nèi)聲場的半解析方法—多模態(tài)導(dǎo)納法[20],構(gòu)建波導(dǎo)中任意位置處的總聲場與入射聲壓的映射關(guān)系, 給出聲波在該位置處聚焦, 即在該位置聲壓幅值達到最大值時的入射條件, 并畫出對應(yīng)的聲場分布.

    2 均勻波導(dǎo)內(nèi)的聲聚焦

    首先考慮二維笛卡爾坐標(biāo)系下均勻剛硬波導(dǎo)中的聲聚焦問題.簡諧聲源(時間因子 e?iωt)位于x=0處.經(jīng)由無量綱化處理其中 ρ ,c,h,ω,k 分別為介質(zhì)密度、聲速、波導(dǎo)高度、角頻率及無量綱波數(shù), x,y分別表示水平方向及垂直方向, 且(x,y)∈ [0,∞)×[0,1], 聲壓滿足無量綱亥姆霍茲方程?2p/?x2+?2p/?y2+k2p=0及邊界條件 ?yp|y=0,1=0.聲壓場的通解可寫為

    本文研究的目的是計算單位輻射聲功率(入射聲能流)條件下, 在空間任意位置 ( x0,y0) 處獲得最大聲壓幅值的最佳入射波, 即計算 pi或 pi,n以獲得在均勻波導(dǎo)中, 由 (1)式可知, 當(dāng)時, 聲壓在 ( x0,y0) 處獲得最大幅值, 其中上標(biāo)“*”表示共軛, L 是聲能流的歸一化系數(shù),即為最佳入射波.無量綱的入射聲能流表達式為

    (3)式中 Np表示傳播模態(tài)數(shù)目, 為使不等式Np? 1+k/π成立的最大整數(shù).聚焦點處的聲壓為

    該聲壓即為給定頻率時, 在 ( x0,y0) 處可達到的聲壓幅值上界.由(4)式可知, 聚焦點處的聲壓主要來自傳播模態(tài)的貢獻, 衰逝模態(tài)的作用為保證該聲壓值在有限頻率下收斂.

    圖1為均勻波導(dǎo)中, 入射波頻率為k=19.1π時, 在 ( x0,y0)=(3.2,0.4) 處產(chǎn)生聚焦的聲壓幅值分布.通過觀察發(fā)現(xiàn), 當(dāng)聚焦點距離邊界較遠(yuǎn)時, 聚焦點附近的聲場與分離常量為0時的奇對稱Weber波束非常相似[21].聲波在均勻波導(dǎo)中的聚焦過程充分利用了邊界的反射作用.

    圖1 均勻波導(dǎo)中的聲聚焦Fig.1.Sound focusing in homogeneous waveguides.

    另外, 從 (4)式可知, 當(dāng)固定頻率以及 x0時,聚焦點處的聲壓僅由決定.根據(jù)余弦函數(shù)的周期性和偶對稱性可得

    其中m為任意整數(shù).( x0,±y0+2m) 可看作是聚焦點 ( x0,y0) 關(guān)于波導(dǎo)上下壁面對稱產(chǎn)生的虛聚焦點(類比虛源), 選擇虛聚焦點中任意一點研究聲聚焦時, 聲波仍會在 ( x0,y0) 處聚焦.

    3 非均勻波導(dǎo)中的聲聚焦

    在均勻波導(dǎo)的聲聚焦問題中, 首先建立了任意位置處聲壓與入射聲壓的展開系數(shù)的映射關(guān)系((1)式), 進而求解使 p (x0,y0) 達到極值的最佳入射波.在非均勻波導(dǎo)中, 聲場與入射波之間不再是簡單的單向傳播的映射關(guān)系, 本文采用多模態(tài)導(dǎo)納法 (multimodal admittance method)[20]建立任意位置聲壓與入射波的展開系數(shù)的映射算子, 進而研究產(chǎn)生聲聚焦的最佳入射波條件并求解聲場.本節(jié)討論水平變截面、含散射體以及聲速垂直變化三種非均勻波導(dǎo)中的聲聚焦現(xiàn)象.

    3.1 水平變截面波導(dǎo)中的聲聚焦

    圖2 水平變截面剛硬波導(dǎo)示意圖Fig.2.Configuration of rigid waveguides with varying cross-sections.

    聲壓滿足的亥姆霍茲方程及邊界條件為

    這里 ψn(x,y) 為局部本征函數(shù), 對應(yīng)均勻波導(dǎo)條件下上邊界等于局部高度 h (x) 時的本征函數(shù), 其表達式為

    其中 εn=2? δn0.然后, 將 (8)式代入 (6)式并投影到局部本征函數(shù)利用本征函數(shù)的正交性可得耦合偏微分方程組:

    其中, p (x) 和 s (x) 分別為模態(tài)展開系數(shù) pn和 sn組成的向量, 上標(biāo)“T”表示轉(zhuǎn)置; I 為單位陣, K(x)為對角陣, 各元素為矩陣A 的表達式為

    由(10)式, 聲場的求解問題轉(zhuǎn)化為耦合偏微分方程組的初值問題, 兩個初值條件為入射條件和輻射條件.然而根據(jù)輻射條件從右向左積分求解(10)式時, 由于衰逝模態(tài)從右向左是指數(shù)發(fā)散的,積分結(jié)果不收斂.因此引入導(dǎo)納矩陣Y, 以規(guī)避衰逝模態(tài)的指數(shù)發(fā)散問題[22,23], 穩(wěn)定求解(10)式.導(dǎo)納矩陣的定義為 s =Yp , 代入(10)式, 求得導(dǎo)納矩陣滿足的Riccati方程:

    Y的初條件為Y(x?xmax)=iK(xmax), 可從xmax至0積分求得任意x處的導(dǎo)納矩陣.為建立x處聲壓與入射波聲壓之間的映射關(guān)系, 引入傳播算子M(x), 定義為 p (x)=M(x)p(0).代入 (10)式可得傳播算子滿足的方程為 ?xM=(Y?A)M , 初條件為 M (0)=I , 從 0 至 xmax積分求解.導(dǎo)納矩陣Y和傳播算子M可利用四階Magnus積分格式進行數(shù)值計算, 詳細(xì)過程見文獻[24]和[25].

    定義模態(tài)域的反射矩陣R滿足pr(0)=Rpi(0), 其中 pi(0),pr(0) 是 x =0 處的入射波和反射波聲壓的模態(tài)展開系數(shù).x =0 處總聲壓的模態(tài)展開系數(shù)可寫為 p (0)=pi(0)+pr(0) , 又由s(0)=iK(0)(I?R)pi(0)=Y(0)p(0), 可得

    定義透射矩陣T滿足pt(xmax)=p(xmax)=Tpi(0),其中 pt(xmax) 為 xmax處透射波聲壓的模態(tài)展開系數(shù), 根據(jù)傳播算子 M (x) 的定義式, 可得

    至此, 反射矩陣R, 透射矩陣T以及傳播算子M(x)均可計算.接著可以給出散射區(qū)域總聲壓的模態(tài)展開系數(shù) pn(x)=M(x)(I+R)pi(0) , 以及透射區(qū)域聲壓的模態(tài)展開系數(shù)

    再結(jié)合(8)式, 最終得到變截面波導(dǎo)內(nèi)任意點(x0,y0)處聲壓與入射聲壓的模態(tài)展開系數(shù) pi(0) 的映射關(guān)系:

    其中, QT(x0,y0) 為映射向量; Ψ(x0,y0)=[ψ0(x0,y0),ψ1(x0,y0),···,ψN?1(x0,y0)]T為 ( x0,y0) 處局部本征函數(shù)構(gòu)成的向量; N為模態(tài)的截斷數(shù), 要求N>max[Np(x)], 即截斷數(shù)的選取需大于波導(dǎo)中的傳播模態(tài)數(shù)量.(15)式可寫作希爾伯特空間中的向量內(nèi)積形式, 即 p (x0,y0)= ?Q?(x0,y0), pi(0)? ,其中上標(biāo)“H”代表共軛轉(zhuǎn)置.使內(nèi)積結(jié)果最大的 pi(0) 就是在 ( x0,y0) 處產(chǎn)生聲聚焦的最佳入射波的模態(tài)展開系數(shù), 即

    其中 L為聲能流的歸一化系數(shù), 計算方法見(2)式.因此, 最佳入射波可寫為

    綜上, 可采用多模態(tài)導(dǎo)納法將復(fù)雜波導(dǎo)中的聲場求解問題轉(zhuǎn)化為計算模態(tài)域的展開系數(shù)問題, 再引入導(dǎo)納矩陣 Y (x) 和傳播算子 M (x) , 可以直接構(gòu)建總聲場與入射波的展開系數(shù)的映射關(guān)系, 這種直接明確的映射關(guān)系將空間任意位置處的聲壓表示為映射向量 Q?(x0,y0) 與入射波模態(tài)展開系數(shù)pi(0)的內(nèi)積的形式, 求解聲聚焦的問題最終轉(zhuǎn)化為求解使該內(nèi)積((15)式)達到極值的最佳入射波問題.需要指出的是, 對于變截面波導(dǎo), 本文選取的局部本征函數(shù) ψn(x,y) 并不滿足實際的上邊界條件, 也就 是 ψn(x,y) 在 上邊 界滿 足的 邊 界 條 件 是?ψn/?y=0, 而不是真正的邊界條件 ? ψn/?n=0 ,這使(8)式的收斂速度不高( 1 /N2).文獻[26]提出了改進多模態(tài)導(dǎo)納法, 用以提高聲場解的收斂性(提高為 1 /N4).方法為構(gòu)造一階新的本征函數(shù)ψ?1(x,y), 使之既與原局部本征函數(shù) ψn(x,y)(n=0,1,···,N? 1)正交, 又在上邊界滿足非齊次Neumann邊界條件.但是, 由于構(gòu)建的本征函數(shù)ψ?1(x,y)在傳播方向上表現(xiàn)為衰逝模態(tài)[26], 且由(4)式, 發(fā)現(xiàn)衰逝模態(tài)對聲聚焦的作用很小, 因此可忽略 ψ?1(x,y) 對聲聚焦的貢獻, 故本文選取的局部本征函數(shù) ψn(x,y) 足夠滿足聲聚焦分析的要求.

    總而言之, 對于非均勻波導(dǎo), 只要求得散射區(qū)域內(nèi)的導(dǎo)納矩陣 Y (x) 和傳播算子 M (x) , 進而得到反射矩陣R及透射矩陣T, 即可給出聲壓與入射聲壓的模態(tài)展開系數(shù)的映射關(guān)系, 最終計算在任意點產(chǎn)生聚焦的最佳入射波.

    圖3給出了利用上述方法計算變截面波導(dǎo)中, 在不同位置處產(chǎn)生聚焦的聲場.入射波頻率為 k =29.1π , 波導(dǎo)上邊界表達式為h(x)=0.8+0.2cos(2πx/3), xmax=3.圖3(a)和 3(b)中的聚焦點分別位于透射區(qū)域和散射區(qū)域, 坐標(biāo)為(x0,y0)=(3.2,0.9)及 ( 1.6,0.2).圖3(c)及圖3(d)插圖中的藍(lán)色實線分別對應(yīng)圖3(a)及圖3(b)中的入射波形,均為由(16)式計算得到的最佳入射聲波, 黑色點線為平面入射波; 主圖中的藍(lán)色實線為固定 x0時,聲壓幅值隨高度的變化曲線, 即 | p (x0,y)| ; 黑色點線為入射波是平面波時, 對應(yīng)的 | p (x0,y)|.如圖3所示, 當(dāng)入射聲波是最佳入射波時, 不論聚焦點在散射區(qū)域還是透射區(qū)域, 聲波利用邊界的散射作用, 均在對應(yīng)點處發(fā)生了聚焦, 并且聚焦點處的聲壓幅值明顯大于入射波是平面波時的聲壓幅值, 聚焦效果良好.

    圖3 (a)和 (b)為變截面波導(dǎo)分別在 ( x0,y0)=(3.2,0.9) (透射區(qū)域)及 ( 1.6,0.2) (散射區(qū)域)處產(chǎn)生聚焦的聲場; (c)和 (d)主圖中的藍(lán)色實線分別為(a)和(b)中 x 0 處的聲壓幅值隨高度方向的分布, 黑色點線為 pi= Λψ0(y) (平面波)時 x0 處的聲壓幅值分布; 插圖中的藍(lán)色曲線和黑色點線分別為最佳入射波及平面波的幅值曲線Fig.3.Acoustic focusing field in the waveguide as calculated by the present method.The foci are located at (a) (x0,y0)=(3.2, 0.9) in transmission region and (b) ( 1.6,0.2) in scattering region, respectively.The blue solid lines in (c) and (d) are |p(x0,y)|corresponding to (a) and (b), respectively, and the black dotted lines are | p (x0,y)| generated by pi= Λψ0(y) (plane wave).The insets plot the modulus of the corresponding incident waves.

    根據(jù)線性疊加原理, 可以實現(xiàn)非均勻波導(dǎo)中的多點聲聚焦.選取多個聚焦點位置, 利用(16)式分別獲得對應(yīng)的最佳入射波.然后將這些入射波求和構(gòu)建新的入射波并將其輸入至波導(dǎo)中, 對應(yīng)聲場則產(chǎn)生多點聲聚焦效應(yīng).圖4(a)給出多點聲聚焦的聲場, 其中頻率的選取及波導(dǎo)結(jié)構(gòu)與圖3一致, 聚焦點位于 ( 3.2,0.9) 及 ( 3.2,0.1).圖4(b)畫出了疊加后的總?cè)肷渎晧旱姆捣植?圖4(c)中的黑色虛線、紅色點劃線和藍(lán)色實線分別為只在 ( 3.2,0.1) 處產(chǎn)生單點聚焦、只在 ( 3.2,0.9) 處產(chǎn)生單點聚焦和同時在 ( 3.2,0.1) 和 ( 3.2,0.9) 處產(chǎn)生雙點聚焦時x=3.2處的聲壓幅值隨y的分布.可以看出雙點聚焦時各個聚焦點處的聲壓幅值均低于單點聚焦的情況, 這是符合能量守恒定律的.雙點聚焦時各個聚焦點處的聲壓幅值均明顯大于其他位置處的聲壓幅值, 說明輸入計算得到的總?cè)肷洳? 可以實現(xiàn)良好的多點聲聚焦效果.

    3.2 含散射體波導(dǎo)中的聲聚焦

    簡單起見, 考慮波導(dǎo)中僅存在一個散射體時的聲聚焦問題, 波導(dǎo)模型如圖5所示, 散射體邊界為[y=a(x),y=b(x)], 波導(dǎo)主介質(zhì)的密度及聲速為ρ1,c1, 散射體內(nèi)部介質(zhì)的密度及聲速為 ρ2,c2, 且ρ1,ρ2,c1,c2均為常量.波導(dǎo)介質(zhì)區(qū)域與散射體占據(jù)區(qū)域分別用 ?1和 ?2表示.聲壓滿足的無量綱亥姆霍茲方程為

    其中, 波數(shù) k =ωh/c1, 1 +η=c1/c2為折射率.邊界條件和連續(xù)性條件分別為

    圖4 (a) 變截面波導(dǎo)中的雙點聚焦聲場, 聚焦點為 ( 3.2,0.9) 及 ( 3.2,0.1) ; (b) 最佳入射聲壓幅值分布; (c) 藍(lán)色實線為 (a) 中聲場在 x =3.2 處的聲壓幅值分布; 紅色點劃線表示聲波在 ( 3.2,0.9) 處單點聚焦時的聲壓幅值分布, 與圖3(c)中藍(lán)色曲線一致; 黑色虛線為聲波在 ( 3.2,0.1) 處單點聚焦時的聲壓幅值分布.頻率和波導(dǎo)幾何參數(shù)與圖3一致Fig.4.(a) Sound two-point focusing field in the waveguide with varying cross-section, the foci are located at ( 3.2,0.9) and(3.2,0.1); (b) modulus of the optimal incident pressure; (c) the blue solid line represents | p (3.2,y)| in (a); the red dot-dashed line shows | p (3.2,y)| when the wave focus only at ( 3.2,0.9) , which is same as the blue solid line in Fig.3(c); and the black dashed line shows | p (3.2,y)| when the wave focus only at ( 3.2,0.1).The frequency and geometries of the waveguide are same as Fig.3.

    圖5 含散射體剛硬波導(dǎo)示意圖Fig.5.Configuration of rigid waveguides involving a scatterer.

    然后將方程(18)變換為模態(tài)展開系數(shù) pn(x) 滿足的模態(tài)域方程[20]:

    (19)式中K與(10)式中一致, E和F的表達式為

    其中

    令 s =E?p/?x=Yp , (19)式可寫成耦合偏微分方程組形式:

    此時導(dǎo)納矩陣Y與傳播算子M滿足的方程和初條件分別為

    二者依然可利用四階Magnus積分格式進行數(shù)值求解[24,25].同樣地, 依照 (13)—(16)式, 可計算得到含散射體波導(dǎo)內(nèi)任意位置處產(chǎn)生聲聚焦的最佳入射聲壓.

    圖6 (a) 含散射體波導(dǎo)中的聚焦聲場, 聚焦點位于 ( 3,0.35) ; (b) pi= Λψ0(y) (平面波) 入射時的聲場; (c) 產(chǎn)生 (a)中聲聚焦的最佳入射波(藍(lán)色實線)與入射平面波(黑色點線)的幅值分布; (d) 藍(lán)色實線與黑色點線分別為(a)與(b)中 x =3 處的聲壓幅值分布Fig.6.(a) Sound focusing field in the waveguide with a scatterer.The focus is located at ( 3,0.35) ; (b) sound field generated by a plane wave pi= Λψ0(y) ; (c) modulus of the pressure of optimal incident wave in (a) and that of the plane incident wave in (b);(d) | p (3,y)| in (a) (blue solid line) and (b) (black dotted line).

    圖6給出了計算含散射體波導(dǎo)中聲聚焦的算例.圖6(a)為在 ( 3,0.35) 處產(chǎn)生聚焦的聲場(聲壓幅值分布), 輸入的入射波是由(16)式計算得到的最佳入射波(圖6(c)中藍(lán)色實線), 圖6(b)是平面波(圖6(c)中黑色點線)入射時的聲場.圖6(d)中的藍(lán)色實線為圖6(a)中聲場在 x =3 處的聲壓幅值分布, 黑色點線對應(yīng)圖6(b)在 x =3 處的聲壓幅值分布.入射波頻率均為 k =29.1π , 散射體為圓形, 圓心位于 ( 1.6,0.45) , 半徑為 0.3.散射體的介質(zhì)參 數(shù) 設(shè) 為 ρ1/ρ2=0.01 , η =?0.99 ( c1/c2=0.01 ),近似剛硬, 因此散射體近似看作對聲波僅有反射作用.可以看出, 剛硬散射體的尺寸較大, 對聲波有強烈的反向散射作用, 但是當(dāng)輸入計算得到的最佳入射波時, 聲波利用散射體的散射和邊界的反射作用, 仍然可在目標(biāo)點處發(fā)生聚焦.

    3.3 聲速垂直變化波導(dǎo)中的聲聚焦

    考慮波導(dǎo)中介質(zhì)參數(shù)隨空間位置連續(xù)變化的非均勻情況, 對應(yīng)的典型波導(dǎo)為淺海波導(dǎo), 表現(xiàn)為聲速沿垂直方向的連續(xù)變化, 且海面為軟邊界.此時聲壓滿足的亥姆霍茲方程及邊界條件為

    對聲場解進行模態(tài)展開, 令

    其 中 ψn(y) 變 為將展開式代入(24)式并投影到局部本征函數(shù)利用正交性, 可得入射聲壓在任意點(x0,y0)處產(chǎn)生的聲壓為

    其中 pi(0) 為入射聲壓 pi的模態(tài)展開系數(shù), K的表達式為

    其計算方法參照文獻[27].由(25)式可知, 在(x0,y0)產(chǎn)生聲聚焦的最佳入射波的模態(tài)展開系數(shù)為

    由于聲速垂直變化波導(dǎo)在傳播方向上介質(zhì)參數(shù)不變, 聲聚焦現(xiàn)象與均勻波導(dǎo)中的聚焦現(xiàn)象類似, 這里不再贅述.

    接下來考慮在聲速垂直變化及散射體共同作用下的復(fù)雜波導(dǎo)聲聚焦問題, 模型類似圖5, 其中聲速 c1(y) 與 c2(y) 變?yōu)閥的函數(shù).此時聲壓滿足的亥姆霍茲方程為

    對(28)式的求解依然基于多模態(tài)導(dǎo)納法, 詳細(xì)的方程推導(dǎo), 映射算子 M (x) 及透射矩陣T的構(gòu)建和計算見文獻[27].利用映射算子及透射矩陣, 參照(13)—(15)式的分析過程, 構(gòu)建任意點 ( x0,y0) 處聲壓與入射聲壓的模態(tài)展開系數(shù) pi(0) 的映射關(guān)系,其中矩陣K的表達式由(26)式給出.最后根據(jù)(16)式計算產(chǎn)生聲聚焦的最佳入射條件.

    圖7(a)給出利用上述方法計算的淺海波導(dǎo)中同時存在聲速垂直變化和散射體時的聚焦聲場, 聚焦點位于 ( 35h,0.1h) , 對應(yīng)最佳入射波幅值分布見圖7(c)插圖.入射波頻率為 f = ω/2π=200Hz , 波導(dǎo)深度 h =100m , [散射體成山狀, 邊]界表達式為a(x)=h?0.5hexp?(x?20h)2/18h2, b (x)=h ,范圍為 [ 8 h,32h].海水介質(zhì)密度 ρ1=1000kg/m3, 聲速為典型淺海負(fù)聲速梯度 c1=1530?0.1y(m/s)[28],如圖7(b)所示, 散射體介質(zhì)的密度和聲速為ρ2=2ρ1及 c2=2c1(或 η =?0.5 ).圖7(c)給出了圖7(a)中聲場在 x =35h 處的聲壓幅值隨深度的變化曲線(藍(lán)色實線), 并與輸入第一階展開模態(tài)pi= Λψ0(y)(圖7(c)插 圖 內(nèi) 黑 色 虛 線)時 在x=35h處獲得的聲壓(黑色虛線)作對比, 此時聲能流的歸一化系數(shù)L由(29)式計算:

    在負(fù)聲速梯度的作用下, 聲波在傳播時會向海底彎曲, 并且由于海底山狀散射體的反向散射作用, 部分聲能會被限制在 x ∈[0,20h]區(qū)域, 而選擇的聚焦點 ( 35h,0.1h) 既位于透射區(qū)域, 又處在海面附近,直覺上聲波似乎難以傳播至聚焦點所在區(qū)域, 但輸入計算得到的對應(yīng)目標(biāo)點的最佳入射波時, 實現(xiàn)了在該點良好的聲聚焦現(xiàn)象.由于散射體不再近似剛硬, 在圖7(a)中可以發(fā)現(xiàn)散射體內(nèi)部存在聲壓分布, 且密度和聲速的不同改變了聲波的波長.最佳入射聲波在實現(xiàn)聚焦的傳播過程中充分地利用了邊界的反射、散射體的散射以及聲速梯度變化引起的介質(zhì)折射率變化等因素, 最終在目標(biāo)點處獲得了最大聲壓.

    圖7 (a) 負(fù)聲速梯度含散射體淺海波導(dǎo)中的聲聚焦, 波導(dǎo)深度 h =100m , 聚焦點 ( 35h,0.1h) 位于透射區(qū)域; (b) 聲速垂直變化情況; (c) x =35h 處聲壓幅值隨深度的變化曲線(藍(lán)色實線)以及 pi= Λψ0(y) (非平面波)入射時產(chǎn)生的聲壓(黑色虛線)對比,插圖中為最佳入射波及 pi= Λψ0(y) 的幅值分布Fig.7.(a) Sound focusing field in the waveguide with negative sound-speed gradient and a scatterer.The focus is located at(35h,0.1h), where h =100m is the depth; (b)the sound speed profile; (c) | p (35h,y)| (blue solid line) compared with that generated by pi= Λψ0(y).The insetplotsthe modulus of the optimal incident pressure (blue solid line) and Λ ψ0(y) (black dashed line).

    4 討 論

    基于多模態(tài)導(dǎo)納法, 本文提出在非均勻波導(dǎo)中的任意位置處實現(xiàn)聲聚焦的理論分析方法.多模態(tài)導(dǎo)納法可以簡單直觀地將亥姆霍茲方程變換為模態(tài)域的耦合偏微分方程組, 并能夠有效地規(guī)避數(shù)值積分時衰逝模態(tài)可能帶來的指數(shù)發(fā)散問題, 是一種準(zhǔn)確高效求解非均勻波導(dǎo)聲場的手段.

    由于實際中難以把計算得到的精確最佳入射波作為輸入, 故討論當(dāng)稀疏輸入最佳入射波時的聲聚焦效果.以圖3(a)中變截面波導(dǎo)的聲聚焦為例,分別以半波長和單倍波長為間隔對理論得到的最佳入射波(圖3(c))進行空間采樣, 獲得稀疏入射波(圖8(c)).通過數(shù)值積分及適當(dāng)補零得到對應(yīng)的模態(tài)展開系數(shù), 進而計算聲場, 結(jié)果分別對應(yīng)圖8(a)及圖8(b).兩個聲場在聚焦點 x0=3.2 處的幅值分布在圖8(d)中給出, 并與理論值進行對比.由圖8可知, 半波長離散得到的稀疏輸入產(chǎn)生的聚焦效果與理論情況基本一致, 原因是半波長離散后的輸入可包含絕大部分原輸入聲波的信息(奈奎斯特采樣定律), 從而能夠產(chǎn)生良好的聚焦效果.而以單倍波長離散得到的稀疏輸入作為聲源時, 該聲源只能包含理論最佳入射波的部分信息, 故產(chǎn)生的聚焦效果對比理論情況有所下降, 但依然可觀察到明顯的聚焦現(xiàn)象.因此, 在y方向稀疏輸入最佳入射波時,通過合理的離散, 仍可實現(xiàn)良好的聲聚焦效果.

    此外, 考慮當(dāng)最佳入射波的幅值或相位存在誤差時對聲聚焦效果產(chǎn)生的影響.依然以圖3(a)為例, 首先固定最佳入射波的相位, 構(gòu)造一個關(guān)于y的取值在 [ 0.5,1.5]均勻分布的隨機函數(shù), 將其與最佳入射波的幅度相乘(相當(dāng)于對最佳入射波的幅度疊加了一個上限為 ± 50% 的隨機偏差)作為輸入(圖9(c)紅色虛線), 進而求解聲場 (圖9(a)).接著,固定最佳入射波的幅值, 類似地, 構(gòu)造一個在[?π/2,π/2]區(qū)間內(nèi)均勻分布的隨機函數(shù)疊加到最佳入射波的相位上作為輸入(圖9(c)黑色點劃線)并求解聲場(圖9(b)).兩個聲場在聚焦點 x0=3.2 處的幅值分布在圖9(d)中給出, 并與無擾動時的理論值對比.從圖9(a)和圖9(d)可以看出, 當(dāng)最佳入射波的幅值存在較大范圍隨機擾動時, 其產(chǎn)生的聚焦效果與理論情況基本一致.而根據(jù)圖9(b)和圖9(d), 當(dāng)最佳入射波的相位存在隨機擾動時, 聚焦效果有所下降.聲源存在隨機擾動時依然可產(chǎn)生聚焦的原因是波導(dǎo)中的衰逝模態(tài)能夠抑制聲波高階振蕩(與波長對比)成分的傳播, 所以即便入射聲波的幅值或相位存在隨機誤差, 經(jīng)過衰逝模態(tài)的修正, 擾動后的入射波依然會在目標(biāo)位置處實現(xiàn)聲聚焦效應(yīng).圖9表明聲聚焦對幅度的擾動具有強魯棒性, 對相位的擾動兼具一定的穩(wěn)健性.

    圖8 y方向稀疏輸入對聚焦結(jié)果的影響 (a) 對最佳入射波進行半波長采樣后的聚焦聲場;(b) 對最佳入射波進行單倍波長采樣后的聚焦聲場;(c)采樣后的入射波幅值分布; (d) 聚焦點 x0 處的聲壓幅值分布.藍(lán)色實線為理論值, 與圖3(c)中的藍(lán)色曲線一致Fig.8.Sound focusing fields when the optimal incident wave is discretized: (a) Half-wavelength spacing; (b) single-wavelength spacing; (c) the moduli of the two spaced incident waves; (d) the red dashed line and the black dot-dashed line are the corresponding|p(3.2,y)| generated by the incident waves in (c).The blue solid line is the theoretical result which is same as that in Fig.3(c).

    圖9 (a) 最佳入射波的幅值存在隨機擾動時的聚焦聲場; (b) 最佳入射波的相位存在隨機擾動時的聚焦聲場; (c) 紅色虛線與黑色點劃線分別為幅值擾動與相位擾動后的入射波幅值分布; (d) 紅色虛線與黑色點劃線分別為(c)中的入射波在聚焦點 x0 處產(chǎn)生的聲壓幅值分布, 藍(lán)色實線為理論值, 與圖3(c)中的藍(lán)色曲線一致Fig.9.Sound focusing fields when (a) the moduli and (b) the arguments of the optimal incident wave are perturbed; (c) the red dashed line is the incident wave with perturbed moduli, and the black dot-dashed line is that with perturbed arguments; (d) the red dashed line and the black dot-dashed line are the corresponding | p (3.2,y)| generated by the incident waves in (c).The blue solid line is the result without perturbation which is same as that in Fig.3(c).

    另外, 本文實現(xiàn)的聲波空間聚焦是輸入最佳入射條件時, 相應(yīng)位置處聲壓的幅值達到最大.最佳入射條件通過推導(dǎo)聲壓場與入射波的映射關(guān)系, 即映射向量 QT(x0,y0) 給出, 該映射關(guān)系在透射區(qū)域用透射矩陣描述, 在散射區(qū)域由傳播算子表征.在光學(xué)領(lǐng)域, 有很多學(xué)者研究優(yōu)化透射理論(optimal wave transmission), 對透射矩陣進行處理, 實現(xiàn)波的最大能量透射[29,30].本文的方法結(jié)合優(yōu)化透射理論, 可以為實現(xiàn)非均勻波導(dǎo)中全透射聚焦提供可行性, 使波既能全透射(無反射)地穿過散射區(qū)域, 又能產(chǎn)生空間聚焦.

    5 結(jié) 論

    本文分析了水平變截面、含散射體以及聲速垂直變化三種非均勻波導(dǎo)中的聲聚焦現(xiàn)象, 叢理論上推導(dǎo)出使聲波在波導(dǎo)中任意位置處聚焦的入射條件.采用多模態(tài)導(dǎo)納法, 將聲壓及其水平方向偏導(dǎo)數(shù)用(局部)本征模態(tài)作為基底展開, 代入亥姆霍茲方程中進行投影, 利用本征模態(tài)的正交性, 最終使聲壓滿足的亥姆霍茲方程轉(zhuǎn)化為模態(tài)展開系數(shù)的耦合偏微分方程組.引入導(dǎo)納矩陣和傳播算子,計算反射矩陣和透射矩陣, 給出任意一點處聲壓與入射聲壓的映射關(guān)系, 將聚焦問題轉(zhuǎn)變?yōu)榧墧?shù)求和或希爾伯特空間內(nèi)積的極值問題, 發(fā)生聚焦時的入射波模態(tài)展開系數(shù)為映射向量的共軛轉(zhuǎn)置.結(jié)果表明, 當(dāng)入射聲波是計算得到的精確的最佳入射聲波時, 可在均勻、非均勻波導(dǎo)中實現(xiàn)明顯的單點或多點聲波聚焦現(xiàn)象.并且當(dāng)稀疏輸入最佳入射波, 或最佳入射波的幅度或相位存在誤差時, 仍可實現(xiàn)較好的聲聚焦效果.聲聚焦對幅度的擾動具有強魯棒性, 對相位的擾動有一定的穩(wěn)健性.聲聚焦的過程充分利用了波導(dǎo)的結(jié)構(gòu)及非均勻性對聲波的散射作用, 例如邊界的反射、散射體或邊界起伏帶來的前向后向散射以及聲速梯度變化引起的介質(zhì)折射率連續(xù)變化等.該方法在聲保密通信、超聲醫(yī)療及無損檢測等領(lǐng)域具有潛在的應(yīng)用價值.

    猜你喜歡
    聚焦點散射體入射波
    SHPB入射波相似律與整形技術(shù)的試驗與數(shù)值研究
    振動與沖擊(2022年6期)2022-03-27 12:18:26
    一種基于單次散射體定位的TOA/AOA混合定位算法*
    二維結(jié)構(gòu)中亞波長缺陷的超聲特征
    無損檢測(2019年11期)2019-11-20 07:07:50
    關(guān)愛聚焦點
    關(guān)愛聚焦點
    關(guān)愛聚焦點
    關(guān)愛聚焦點
    瞬態(tài)激勵狀態(tài)下樁身速度以及樁身內(nèi)力計算
    高斯波包散射體成像方法
    城市建筑物永久散射體識別策略研究
    城市勘測(2016年2期)2016-08-16 05:58:24
    亚洲第一区二区三区不卡| 国产一区有黄有色的免费视频 | 爱豆传媒免费全集在线观看| 国语对白做爰xxxⅹ性视频网站| 国产一级毛片在线| 久久这里有精品视频免费| 免费大片黄手机在线观看| 亚洲av日韩在线播放| 久久久久久久久大av| 欧美激情在线99| 日本黄色片子视频| 天天一区二区日本电影三级| 男女边摸边吃奶| 亚洲综合色惰| 97热精品久久久久久| 亚洲av电影不卡..在线观看| 亚洲人成网站在线观看播放| 精品一区二区三区人妻视频| 亚洲自偷自拍三级| 一本久久精品| 丝袜喷水一区| 国产精品久久久久久久电影| 最新中文字幕久久久久| 91久久精品国产一区二区成人| 午夜精品在线福利| 男女那种视频在线观看| 身体一侧抽搐| 麻豆国产97在线/欧美| 久久国内精品自在自线图片| 成人综合一区亚洲| 国产熟女欧美一区二区| 精华霜和精华液先用哪个| 18禁在线播放成人免费| 国产一级毛片在线| 一个人看的www免费观看视频| 日本午夜av视频| 国产精品精品国产色婷婷| 水蜜桃什么品种好| 欧美日本视频| 久久99精品国语久久久| 日本爱情动作片www.在线观看| 麻豆av噜噜一区二区三区| 国产午夜精品一二区理论片| 国产极品天堂在线| 舔av片在线| 最近2019中文字幕mv第一页| 欧美日本视频| 美女主播在线视频| 成人无遮挡网站| 成年av动漫网址| www.av在线官网国产| 91av网一区二区| 综合色丁香网| 蜜桃久久精品国产亚洲av| 色综合色国产| av免费在线看不卡| 国产精品嫩草影院av在线观看| 久久久a久久爽久久v久久| 成人亚洲欧美一区二区av| 成人国产麻豆网| 午夜精品在线福利| 亚洲在线观看片| 亚洲电影在线观看av| av免费在线看不卡| 免费观看的影片在线观看| videos熟女内射| 老司机影院毛片| 九九久久精品国产亚洲av麻豆| 欧美激情在线99| 国产v大片淫在线免费观看| 国产伦一二天堂av在线观看| 又大又黄又爽视频免费| 插阴视频在线观看视频| 国产精品麻豆人妻色哟哟久久 | av女优亚洲男人天堂| 成人综合一区亚洲| 国产黄频视频在线观看| 一级毛片久久久久久久久女| 国产久久久一区二区三区| 人人妻人人澡人人爽人人夜夜 | 男人狂女人下面高潮的视频| 久久久久久久久久人人人人人人| 亚洲国产欧美在线一区| 久久人人爽人人片av| 国产成人午夜福利电影在线观看| 美女大奶头视频| 啦啦啦中文免费视频观看日本| 特大巨黑吊av在线直播| 插逼视频在线观看| 日韩三级伦理在线观看| 免费黄色在线免费观看| 日本熟妇午夜| 久久久久久伊人网av| 午夜免费激情av| 黄色欧美视频在线观看| 亚洲精品日韩在线中文字幕| 一区二区三区四区激情视频| 直男gayav资源| 精品国产一区二区三区久久久樱花 | 最近最新中文字幕大全电影3| 亚洲色图av天堂| 久久99热这里只有精品18| 欧美xxxx黑人xx丫x性爽| 国产综合懂色| 亚洲欧美一区二区三区黑人 | 国产久久久一区二区三区| 欧美日韩综合久久久久久| 国产 亚洲一区二区三区 | 久久韩国三级中文字幕| 国产在视频线在精品| 欧美日韩一区二区视频在线观看视频在线 | 国产精品一区二区在线观看99 | 人妻系列 视频| 欧美变态另类bdsm刘玥| 五月伊人婷婷丁香| 国产老妇伦熟女老妇高清| 午夜免费观看性视频| 久久热精品热| 69人妻影院| 精品99又大又爽又粗少妇毛片| 黄色日韩在线| 亚洲av免费在线观看| 能在线免费观看的黄片| 免费av观看视频| 国产真实伦视频高清在线观看| 中国国产av一级| 晚上一个人看的免费电影| 欧美激情久久久久久爽电影| 免费电影在线观看免费观看| 国产人妻一区二区三区在| 97超碰精品成人国产| 国产黄色免费在线视频| 国产精品久久视频播放| 人体艺术视频欧美日本| 久久精品综合一区二区三区| 中文字幕免费在线视频6| 欧美激情久久久久久爽电影| 亚洲18禁久久av| 国产精品美女特级片免费视频播放器| 久久久久免费精品人妻一区二区| 精品久久久久久久久亚洲| 亚洲婷婷狠狠爱综合网| 大又大粗又爽又黄少妇毛片口| 国产成人精品福利久久| 精品久久久久久久久av| 精品久久久久久成人av| 国产精品综合久久久久久久免费| 日日干狠狠操夜夜爽| 亚洲国产精品国产精品| 日韩精品青青久久久久久| 国产亚洲最大av| 男人狂女人下面高潮的视频| 国内精品美女久久久久久| 国产精品三级大全| 亚洲美女搞黄在线观看| 男女啪啪激烈高潮av片| 卡戴珊不雅视频在线播放| 尾随美女入室| 插阴视频在线观看视频| 嫩草影院新地址| 精品久久久精品久久久| 亚洲在线观看片| 亚洲精品日韩在线中文字幕| 一区二区三区免费毛片| 免费av观看视频| 欧美日韩精品成人综合77777| 色播亚洲综合网| 又爽又黄a免费视频| 亚洲第一区二区三区不卡| 麻豆国产97在线/欧美| 在现免费观看毛片| 99re6热这里在线精品视频| 婷婷色综合大香蕉| 男人爽女人下面视频在线观看| 97超碰精品成人国产| 男女下面进入的视频免费午夜| 国产老妇伦熟女老妇高清| 日韩三级伦理在线观看| 99热6这里只有精品| 国产探花极品一区二区| 精品午夜福利在线看| 简卡轻食公司| 伦理电影大哥的女人| 99久久中文字幕三级久久日本| 简卡轻食公司| 亚洲成色77777| 中文精品一卡2卡3卡4更新| 成人性生交大片免费视频hd| 国产一区二区三区综合在线观看 | 亚洲av国产av综合av卡| 视频中文字幕在线观看| 国产成人精品一,二区| 国产黄色免费在线视频| 最近手机中文字幕大全| 国产高潮美女av| 美女xxoo啪啪120秒动态图| 成年av动漫网址| 精品一区二区三区人妻视频| 卡戴珊不雅视频在线播放| 日韩视频在线欧美| 久久精品国产自在天天线| 国产精品熟女久久久久浪| 国产午夜精品论理片| 免费av毛片视频| 精品人妻熟女av久视频| 国产精品一区二区在线观看99 | av黄色大香蕉| 亚洲乱码一区二区免费版| 大片免费播放器 马上看| 大又大粗又爽又黄少妇毛片口| 亚洲成人精品中文字幕电影| 国内精品宾馆在线| 少妇高潮的动态图| 国产成人福利小说| 亚洲天堂国产精品一区在线| 亚洲乱码一区二区免费版| 性色avwww在线观看| 久久久精品免费免费高清| 美女主播在线视频| 国产淫片久久久久久久久| 99re6热这里在线精品视频| 亚洲av男天堂| 成人二区视频| 大话2 男鬼变身卡| 亚洲精品久久午夜乱码| 美女被艹到高潮喷水动态| 免费观看a级毛片全部| 女人十人毛片免费观看3o分钟| 亚洲美女搞黄在线观看| 你懂的网址亚洲精品在线观看| 欧美精品国产亚洲| 成人特级av手机在线观看| 99久久精品热视频| 91aial.com中文字幕在线观看| 少妇裸体淫交视频免费看高清| 嫩草影院新地址| 亚洲av电影不卡..在线观看| 色5月婷婷丁香| 纵有疾风起免费观看全集完整版 | 久久久精品欧美日韩精品| 午夜激情福利司机影院| 国产成人午夜福利电影在线观看| 在线观看av片永久免费下载| 我的女老师完整版在线观看| 婷婷色综合www| 精品久久久噜噜| 日韩一本色道免费dvd| 精品一区在线观看国产| 三级国产精品片| 欧美+日韩+精品| 黄片wwwwww| 美女内射精品一级片tv| 天堂影院成人在线观看| 日本-黄色视频高清免费观看| 天堂√8在线中文| 51国产日韩欧美| 日韩成人伦理影院| 日本与韩国留学比较| 日日摸夜夜添夜夜爱| 欧美成人午夜免费资源| 国产高清不卡午夜福利| xxx大片免费视频| 日本欧美国产在线视频| 欧美激情久久久久久爽电影| 我的女老师完整版在线观看| 别揉我奶头 嗯啊视频| 人人妻人人澡欧美一区二区| 水蜜桃什么品种好| 精品熟女少妇av免费看| 精品久久国产蜜桃| 免费黄网站久久成人精品| 麻豆久久精品国产亚洲av| 久久久精品欧美日韩精品| 亚洲国产精品国产精品| 成人亚洲欧美一区二区av| 国产精品1区2区在线观看.| 亚洲18禁久久av| 18禁裸乳无遮挡免费网站照片| av在线老鸭窝| 精品久久久久久久久亚洲| 午夜老司机福利剧场| 天堂av国产一区二区熟女人妻| 又大又黄又爽视频免费| 亚洲精品成人av观看孕妇| 欧美日本视频| 亚洲av免费在线观看| 成人二区视频| 三级男女做爰猛烈吃奶摸视频| 一级毛片电影观看| 老师上课跳d突然被开到最大视频| 人妻一区二区av| 亚洲精品成人久久久久久| 最近手机中文字幕大全| 日韩视频在线欧美| 3wmmmm亚洲av在线观看| 午夜免费激情av| 国产色婷婷99| 高清午夜精品一区二区三区| 一级毛片 在线播放| videossex国产| 联通29元200g的流量卡| 免费观看a级毛片全部| 美女内射精品一级片tv| 午夜久久久久精精品| 夜夜看夜夜爽夜夜摸| 精品久久国产蜜桃| 91精品国产九色| 夫妻性生交免费视频一级片| 亚洲av成人av| 久久人人爽人人片av| 日韩,欧美,国产一区二区三区| 日韩制服骚丝袜av| 男人和女人高潮做爰伦理| 18+在线观看网站| 免费播放大片免费观看视频在线观看| 毛片女人毛片| 日韩在线高清观看一区二区三区| 国内揄拍国产精品人妻在线| 中文乱码字字幕精品一区二区三区 | 亚洲国产欧美在线一区| 亚洲人与动物交配视频| 午夜精品一区二区三区免费看| 街头女战士在线观看网站| 久99久视频精品免费| 日本免费在线观看一区| 欧美xxxx黑人xx丫x性爽| 亚洲,欧美,日韩| 成人特级av手机在线观看| 免费观看精品视频网站| 男女啪啪激烈高潮av片| 国产精品一区www在线观看| 亚洲电影在线观看av| 亚洲国产最新在线播放| 最近中文字幕2019免费版| 91久久精品国产一区二区成人| 国产探花极品一区二区| 久久97久久精品| 国产伦在线观看视频一区| 激情五月婷婷亚洲| 久久这里有精品视频免费| 亚洲自偷自拍三级| 亚洲av成人精品一区久久| 国产伦精品一区二区三区视频9| 久热久热在线精品观看| 亚洲真实伦在线观看| 亚洲精品成人av观看孕妇| 日韩一本色道免费dvd| 久久久久精品性色| 精品国产露脸久久av麻豆 | 日韩欧美国产在线观看| 搡女人真爽免费视频火全软件| 亚洲av不卡在线观看| 91aial.com中文字幕在线观看| 青春草国产在线视频| 亚洲精品日韩在线中文字幕| 国产精品美女特级片免费视频播放器| 欧美 日韩 精品 国产| 两个人的视频大全免费| 一本一本综合久久| 国产精品一区二区性色av| 69人妻影院| 一级毛片aaaaaa免费看小| 日本黄大片高清| av一本久久久久| 成年免费大片在线观看| 亚洲经典国产精华液单| 久久久久久久午夜电影| 女的被弄到高潮叫床怎么办| 一级毛片aaaaaa免费看小| 欧美日韩精品成人综合77777| 国产淫片久久久久久久久| 中文精品一卡2卡3卡4更新| 91精品国产九色| av网站免费在线观看视频 | 女人十人毛片免费观看3o分钟| 亚洲丝袜综合中文字幕| 久久久久久九九精品二区国产| 国产av国产精品国产| a级一级毛片免费在线观看| 热99在线观看视频| 亚洲精品第二区| 永久免费av网站大全| 五月天丁香电影| 麻豆久久精品国产亚洲av| 五月天丁香电影| 纵有疾风起免费观看全集完整版 | 亚洲欧美成人综合另类久久久| 欧美丝袜亚洲另类| 亚洲电影在线观看av| 高清欧美精品videossex| 天堂中文最新版在线下载 | 国产精品麻豆人妻色哟哟久久 | 波多野结衣巨乳人妻| a级毛片免费高清观看在线播放| 欧美97在线视频| 一边亲一边摸免费视频| 在线观看人妻少妇| 亚洲内射少妇av| 秋霞伦理黄片| 狂野欧美激情性xxxx在线观看| 欧美日本视频| 欧美激情在线99| 岛国毛片在线播放| 丝袜美腿在线中文| 亚洲av二区三区四区| 天天躁日日操中文字幕| 熟女电影av网| 毛片一级片免费看久久久久| 美女脱内裤让男人舔精品视频| 国产亚洲av嫩草精品影院| 午夜福利在线观看免费完整高清在| 男插女下体视频免费在线播放| 身体一侧抽搐| 亚洲精品国产av蜜桃| 婷婷色综合大香蕉| 尤物成人国产欧美一区二区三区| 精品午夜福利在线看| 成年女人在线观看亚洲视频 | 国产精品久久久久久精品电影小说 | 国产亚洲最大av| 国产伦理片在线播放av一区| av在线天堂中文字幕| 国产色婷婷99| 男女国产视频网站| 小蜜桃在线观看免费完整版高清| 亚洲av福利一区| 夜夜爽夜夜爽视频| av一本久久久久| 国产伦一二天堂av在线观看| 一夜夜www| 少妇高潮的动态图| 亚洲国产精品成人综合色| 精品不卡国产一区二区三区| 在现免费观看毛片| 国产亚洲5aaaaa淫片| 免费看日本二区| 国产精品av视频在线免费观看| 中文字幕人妻熟人妻熟丝袜美| 99热网站在线观看| 蜜桃久久精品国产亚洲av| 人人妻人人澡欧美一区二区| 欧美高清性xxxxhd video| 久久久久久久久大av| 欧美97在线视频| 久久久久久国产a免费观看| 91狼人影院| 搡老妇女老女人老熟妇| 国产高清不卡午夜福利| 我要看日韩黄色一级片| 熟妇人妻久久中文字幕3abv| 国模一区二区三区四区视频| 国产高潮美女av| 秋霞在线观看毛片| 男女下面进入的视频免费午夜| 观看免费一级毛片| 男女边吃奶边做爰视频| 欧美激情国产日韩精品一区| 午夜激情欧美在线| 伦精品一区二区三区| 男女视频在线观看网站免费| 人妻夜夜爽99麻豆av| 麻豆乱淫一区二区| 亚洲美女视频黄频| 久久精品熟女亚洲av麻豆精品 | 日韩视频在线欧美| 欧美一区二区亚洲| 女人久久www免费人成看片| 最近手机中文字幕大全| 国产黄色免费在线视频| 99热这里只有精品一区| 国产精品一区二区三区四区久久| 18+在线观看网站| 日日摸夜夜添夜夜爱| 欧美三级亚洲精品| 久久精品国产亚洲网站| 国产伦精品一区二区三区视频9| 在线观看人妻少妇| 男女国产视频网站| 国产白丝娇喘喷水9色精品| 高清欧美精品videossex| 亚洲av中文av极速乱| 中文字幕亚洲精品专区| 在线免费观看的www视频| 日韩一区二区视频免费看| 一个人免费在线观看电影| 色网站视频免费| 国产精品久久久久久精品电影小说 | 波多野结衣巨乳人妻| 日韩成人av中文字幕在线观看| 亚洲成人av在线免费| 亚洲国产精品成人久久小说| 免费观看a级毛片全部| 亚洲精品aⅴ在线观看| 国产精品一二三区在线看| 亚洲美女视频黄频| 国产伦精品一区二区三区四那| 搡老妇女老女人老熟妇| 亚洲第一区二区三区不卡| 国产乱来视频区| 又大又黄又爽视频免费| 成年女人在线观看亚洲视频 | www.av在线官网国产| 国产亚洲精品久久久com| 亚洲精品自拍成人| 内地一区二区视频在线| 亚洲内射少妇av| 天天一区二区日本电影三级| 婷婷色综合www| 亚洲精品视频女| 精品99又大又爽又粗少妇毛片| 欧美精品国产亚洲| 国产一区二区三区综合在线观看 | 久久99热这里只频精品6学生| 全区人妻精品视频| 国产精品国产三级国产专区5o| 婷婷六月久久综合丁香| 80岁老熟妇乱子伦牲交| 69av精品久久久久久| 1000部很黄的大片| 女的被弄到高潮叫床怎么办| 国产亚洲91精品色在线| 欧美激情国产日韩精品一区| 亚洲最大成人手机在线| 亚洲乱码一区二区免费版| 亚洲天堂国产精品一区在线| 舔av片在线| 精品不卡国产一区二区三区| 国产乱来视频区| 亚洲综合色惰| 免费av观看视频| 欧美xxxx黑人xx丫x性爽| 亚洲精品亚洲一区二区| 久久精品国产鲁丝片午夜精品| a级毛片免费高清观看在线播放| 三级经典国产精品| 最近的中文字幕免费完整| 国产精品av视频在线免费观看| 国产伦精品一区二区三区四那| 一二三四中文在线观看免费高清| 我的老师免费观看完整版| 欧美丝袜亚洲另类| 亚洲一级一片aⅴ在线观看| 日韩成人伦理影院| 美女cb高潮喷水在线观看| 国产av在哪里看| 老女人水多毛片| av一本久久久久| 国产伦一二天堂av在线观看| 欧美激情在线99| 免费观看在线日韩| 亚洲精品中文字幕在线视频 | 一个人观看的视频www高清免费观看| 神马国产精品三级电影在线观看| 国产中年淑女户外野战色| 国产精品av视频在线免费观看| 亚洲人成网站在线观看播放| 六月丁香七月| 国产黄色免费在线视频| 在线观看人妻少妇| 日韩成人av中文字幕在线观看| 人妻夜夜爽99麻豆av| 亚洲精品国产av蜜桃| 久久精品综合一区二区三区| 亚洲成色77777| 久久精品国产亚洲网站| 日韩欧美三级三区| 亚洲,欧美,日韩| 国产亚洲av嫩草精品影院| 日本午夜av视频| 啦啦啦韩国在线观看视频| 国产免费视频播放在线视频 | 亚洲久久久久久中文字幕| 久久精品综合一区二区三区| 亚洲成人av在线免费| 美女高潮的动态| 国产精品人妻久久久影院| 十八禁网站网址无遮挡 | 少妇人妻一区二区三区视频| 一级爰片在线观看| 成人午夜高清在线视频| 男女边摸边吃奶| 丝瓜视频免费看黄片| 亚洲最大成人av| 插逼视频在线观看| 女人被狂操c到高潮| 国产男人的电影天堂91| 九草在线视频观看| 少妇人妻精品综合一区二区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 一级毛片aaaaaa免费看小| 男女视频在线观看网站免费| 汤姆久久久久久久影院中文字幕 | 久99久视频精品免费| 日韩制服骚丝袜av| 伦精品一区二区三区| 久久久久精品性色| 免费av毛片视频| 一个人看视频在线观看www免费| 精品人妻熟女av久视频| 国产成人免费观看mmmm| 一个人看视频在线观看www免费| 精品人妻熟女av久视频| 日日撸夜夜添| 日韩,欧美,国产一区二区三区| 能在线免费观看的黄片| 久久久久久久久大av| 日本熟妇午夜| 少妇人妻精品综合一区二区| 亚洲精品色激情综合| 精品酒店卫生间| 三级毛片av免费| 国产欧美另类精品又又久久亚洲欧美| 国产毛片a区久久久久| 国产午夜福利久久久久久| 欧美zozozo另类|