?

基于壓電分流技術的PT對稱梁散射特性研究

2022-11-14 18:39王剛羅彩明張曉東
湖南大學學報·自然科學版 2022年8期

王剛 羅彩明 張曉東

摘要:為了解決宇稱-時間(Parity-Time,PT)對稱結構存在結構復雜、奇異點難以調諧的問題,基于壓電分流技術,設計一種針對彎曲波的PT對稱梁.推導PT對稱條件,基于等效介質法和有限元法,驗證所設計增益單元、損耗單元的等效參數滿足該PT對稱條件,并通過改變諧振頻率和分流電阻研究了奇異點的可調性.通過傳遞矩陣法和有限元法,對PT對稱梁的散射特性進行分析,討論單向無反射點與奇異點之間的關系.理論計算和仿真結果表明,PT對稱梁有包括511 Hz和520.5 Hz在內的多個奇異點.彎曲波頻率為511 Hz的彎曲波從右端入射時反射系數趨于零,而彎曲波頻率為520.5 Hz時,它必須施加在左端從而達到完全透射而不產生反射.

關鍵詞:壓電材料;彎曲波;宇稱-時間對稱;奇異點;單向無反射;散射矩陣

中圖分類號:TH113文獻標志碼:A

Study on Scattering Characteristic of PT Symmetric Beam Based on Piezoelectric Shunting Technology

WANG Gang,LUO Caiming,ZHANG Xiaodong

(State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,Hunan University,Changsha 410082,China)

Abstract:To solve the problems of complex structure and difficulty in tuning the exceptional points in the existing PT-symmetric structures,a PT-symmetric beam for flexural waves is designed based on piezoelectric shunting technology. Firstly,the PT-symmetric condition is derived. Then,based on the effective medium method and finite element simulation,it is verified that the effective parameters of the gain and loss unit meet the PT-symmetric condition,and the tunability of exception points is studied by changing the resonant frequency and the shunting resistance. Finally,the scattering property of the PT-symmetric beam is derived by the transfer matrix method and finite element simulation,and the relationship between exceptional points and unidirectional non-reflection is illustrated. The calculated and simulated results show that the PT-symmetric beam has several exceptional points including 511 Hz and 520.5 Hz. When the incident flexural wave of 511 Hz is applied at the right side of the PT-symmetric beam,the reflection coefficient is close to zero. However,when the frequency of the incident flexural waves changes to 520.5 Hz,it should be applied on the left side of the PT-symmetric beam to gain an entire transmission without reflection.

Key words:piezoelectric materials;flexural wave;Parity-Time(PT)symmetric;exceptional point;unidirectional non-reflection;scattering matrix

宇稱-時間(Parity-Time,PT)對稱由Bender和Boettcher于1998年提出[1],他們發現在滿足PT對稱的條件下,系統的哈密頓量即使是非厄米的也可以具有實數的本征值.系統具有PT對稱性的條件是勢函數共軛對稱,即勢函數的實部偶對稱、虛部奇對稱.PT對稱系統的一個重要性質為PT對稱性的自發性破缺,具體表現是:當勢函數虛部在某一閾值以下時,哈密頓量的本征值全部為實數,系統處于PT對稱相;當勢函數虛部越過該閾值時,哈密頓量的本征值開始有復數產生,此時對應于PT破缺相.這一閾值也被稱為相變點或奇異點(Exceptional Point).PT對稱將量子力學拓展到了非厄米范圍,此后它便成了量子力學中的一個重要研究方向.現階段非厄米PT對稱的研究工作已經拓展到光學和聲學領域中,表現出許多奇特的物理現象.在光學領域的研究中,調節折射率的分布,使折射率在傳播方向上滿足實部偶對稱、虛部奇對稱,即可構造PT對稱光學系統[2-6];而在聲學領域的研究中,當聲學系統滿足一定條件時,它也會有PT對稱性,例如Fleury等人[7]通過調節兩個麥克風連接的阻抗使得麥克風的等效質量密度共軛對稱,從而利用PT對稱性實現了一種聲學隱身傳感器.

PT對稱性的關鍵點在于引入均衡的增益和損耗,而目前,固體介質中針對彎曲波的PT對稱性的研究相對較少,其中一個重要的原因是自然介質難以實現能量損耗和能量增益之間的平衡.壓電材料由于能夠實現機械能與電能之間的轉化,并且具有易于調控的特點,有望解決這一問題.Christensen等人[4]利用壓電半導體中的聲電效應構造了聲子PT對稱系統,并分析了其中的單向無反射現象.Hou等人[8]基于壓電分流單元,提出了一種可調PT對稱系統,通過改變分流電路的阻抗便可以調節PT對稱系統奇異點出現的頻率.上述研究工作所提出的PT對稱系統均是針對縱波而言的,所獲得的研究成果為PT對稱系統的應用奠定了理論基礎.而與縱波相比,彎曲波的理論模型更為復雜,因此針對彎曲波的PT對稱系統同樣值得研究.

本文首先基于歐拉梁的假設,推導了針對彎曲波的PT對稱條件;然后基于該條件,利用壓電分流單元設計一種PT對稱梁,并采用等效介質法和有限元仿真,得到增益單元、損耗單元的等效質量密度和等效彎曲剛度;最后通過傳遞矩陣法和有限元仿真證明了PT對稱梁中的單向無反射現象,并通過計算散射矩陣的特征值和特征向量說明了單向無反射現象源于奇異點的存在.本文所提出的PT對稱梁對于彎曲波有許多潛在的應用,包括增強傳感、彎曲波放大和非對稱控制等.

1PT對稱梁的設計與分析

1.1彎曲波PT對稱條件

根據歐拉梁假設,彎曲波沿x軸傳播時滿足控制方程[9]:

式中:w(x)、Deff(x)、ρeff(x)和hb分別表示梁的橫向位移、等效彎曲剛度、等效質量密度和厚度.若梁滿足PT對稱性,則根據式(1)可以得到:

公式(3)和公式(4)表明,若滿足PT對稱性,則要求等效質量密度ρeff(x)和等效彎曲剛度Deff(x)關于原點共軛對稱.

1.2PT對稱梁設計

壓電分流單元已經被證明可以有效地改變基體材料的等效彎曲剛度而幾乎不改變其等效質量密度[10],因此,可以通過設計合理的壓電分流單元使其滿足公式(3)和(4)所提出的等效質量密度和等效彎曲剛度共軛對稱條件.圖1(a)為所設計的PT對稱梁示意圖,基體梁上表面貼有兩塊壓電片,兩塊壓電片分別連接由正、負電阻與電感并聯組成的分流電路,從而與基體梁部分構成損耗單元、增益單元.分流電路中,正、負電阻的作用分別是衰減和放大彎曲波. 正電阻Rsh和電感Lsh并聯后的分流阻抗ZL、負電阻-Rsh與電感Lsh并聯后的分流阻抗ZG分別為:

式中:ω為角頻率.

PT對稱梁中的負電阻可以通過圖1(b)所示的Non-Foster電路來實現[11-12],大電感可以用圖1(c)所示的Antoniou's模擬電感電路實現[13-14].此外,分流電路中的電感會與壓電片的固有電容產生諧振,從而使得增益單元、損耗單元的等效彎曲剛度Deff(x)與彎曲波頻率相關,該現象在第2節進一步說明.

1.3PT對稱梁的特性分析

本節基于等效介質理論計算增益單元、損耗單元的等效參數,并通過傳遞矩陣法,針對圖1(a)所示PT對稱梁,計算彎曲波分別從左端和右端兩種入射情況下的透射系數和反射系數.同時,使用有限元軟件COMSOL Multiphysics進行有限元仿真,以驗證理論計算結果的準確性.

如圖2所示,令基體梁的厚度及彈性模量分別為hb、Eb,壓電片的厚度和長度分別為hp、Lp.假設壓電片除了垂直于x軸的端面外,其余表面均為自由狀態,且壓電片僅在z軸存在極化,因此壓電片的本構方程可簡化為[15-16]:

另外,圖2所示單元內任意一點x方向的位移u(x,z)為縱波位移u0(x)與彎曲波位移wp(x)的組合:

由于壓電片的厚度與基體梁的厚度處于同一個

數量級,因此電場強度E3在z方向不能被認為是恒定的.本文采用文獻[17]提出的假設,認為電場強度E3在z方向呈線性分布,即電場強度E3和電勢V(x,z)可以用如下公式表述:

電位移D3滿足控制方程[18]:

聯立公式(6)~公式(10),便可以求出系數a(x)、b(x)、c(x)和壓電片上表面的電壓Vup[17]:

以基體梁的中性面為基準面,對應力求積分,可以得到單元x-y截面上的彎矩M、法向力N和剪切力T:

其中系數I、J、K、F和G通過下式計算:

基于公式(15)~公式(17),圖2所示單元的運動控制方程為:

和彎曲波wp(x)的表達式分別為:

對于基體梁中未貼壓電片的部分及貼有壓電片的部分,本文分別定義向量Yb和Yp

根據公式(22)~公式(28),公式(15)~公式(17)可以用以下矩陣形式表述:

Yp=Bp(x)Ap(31)

Yb=BbAb(32)

其中:

根據連續性條件,結合公式(31)和公式(32),可以得到增益單元或損耗單元的傳遞矩陣Mi,其中i=G表示增益單元,i=L表示損耗單元:

式中:Pp為對角矩陣.

1.3.1等效參數

本文僅對波長遠大于單元長度Lp的低頻彎曲波進行分析,采用等效介質理論計算結構的等效參數[20].因為增益單元和損耗單元的等效參數求解原理相同,對應的等效參數計算過程相同,所以本文僅以損耗單元為例介紹等效參數求解過程.

令x=0,根據公式(36),Yb(Lp)與Yb(0)之間滿足:

求解等效質量密度ρeff時,需要獲取單元的質量,這要求單元任意一點的橫向加速度均相同.如圖2 (a)所示,損耗單元整體相對于虛線方框所示的原始位置作幅值為A1的橫向簡諧振動,并且限制其左、右邊界的縱向位移和旋轉角度為0.依據上述邊界條件,基于公式(38)可求解向量中Yb(0)、Yb(Lp)的未知元素,從而得到損耗單元的等效質量密度ρeff

為了驗證等效介質法計算結果的正確性,本文針對損耗單元和增益單元,使用COMSOL分別建立了如圖2所示2種情況下的二維有限元模型.有限元模型中,認為單元滿足平面引力假設,施加的邊界條件與等效介質法中的邊界條件一致.

求解等效質量密度ρeff時,分別在左、右邊界上,對單元所受的z方向力求線積分即可得到對應邊界上的剪切力T(0)、T(Lp);而求解等效彎曲剛度Deff時,在右邊界上對單元所受的y方向力矩求線積分可得到力矩M(Lp).得到剪切力和力矩之后,便可以通過公式(39)和公式(40)分別計算得到單元的等效質量密度ρeff及等效彎曲剛度Deff.

等效參數法和有限元仿真中,損耗單元和增益單元的幾何參數及材料參數如表1所示.另外,電感與壓電片的固有電容會產生諧振,因此取諧振頻率f0=500Hz,分流電阻Rsh=1 000Ω.等效參數法和有限元仿真的計算結果在2.1節進行討論.

1.3.2透射系數及反射系數

為了表述方便,本文用下標l表示彎曲波從左端入射,用下標r表示彎曲波從右端入射,將彎曲波從左端入射時的透射系數和反射系數分別定義為左透射系數t1和左反射系數r1.同理,定義彎曲波從右端入射時的右透射系數tr和右反射系數rr.

在PT對稱梁的左端或右端施加一個單位位移,即

式中:L表示兩塊壓電片之間的距離.求解公式(41)便可以分別得到左透射系數t1、左反射系數r1、右反射系數r1、右透射系數tr,進而可以得到彎曲波的散射矩陣S[21]:

為驗證傳遞矩陣法計算結果的正確性,本文使用COMSOL軟件建立如圖3所示的有限元仿真模型,基體梁上含有增益單元、損耗單元.分別在左端激勵點和右端激勵點施加力載荷,產生彎曲波.為了抑制邊界反射,基體梁的兩端各連接一個完美匹配層,其他邊界保持自由狀態.完美匹配層的長度至少為研究波長的1/2[22],在本文中,波長的最大值約為0.384 4 m,則本文中完美匹配層的長度設為0.2 m. 由于激勵點關于原點。對稱,所以觀測點相對原點也要對稱分布.為了避免受到激勵點的影響,觀測點與激勵點之間要保持一定距離.觀測點A和觀測點B的坐標分別為(-0.35 m,0)和(0.35 m,0),左端激勵點和右端激勵點的坐標分別為(-0.5 m,8.0×10-4m)和(0.5m,8.0×10-4m).圖4為試驗裝置圖,在圖4中,梁的兩端附著有藍丁膠,用來減少反射.信號發生器產生的信號經過功率放大器后作用在壓電激勵器上,產生彎曲激勵.激光測振儀用來測量指定點的位移響應,它與信號調理器、數據采集卡、電腦共同組成信號采集處理系統.試驗中所需要的負電阻和電感可以分別通過圖1(b)和圖1(c)所示電路來實現.

傳遞矩陣法和有限元仿真中,壓電片長度Lp=45 mm,兩塊壓電片之間的距離L = 300 mm,分流電阻Rsh=1 900Q,分流電路的諧振頻率f0=500 Hz,其他參數與表1中所列參數相同.有限元仿真模型的長度為2 005 mm.

2結果和討論

2.1等效參數

基于上述等效介質法和有限元仿真,計算得到增益單元和損耗單元的等效質量密度ρeff和等效彎曲剛度Deff.為了表達方便,本文采用歸一化等效質量密度ρeffb和歸一化等效彎曲剛度Deff/Db來表征增益單元、損耗單元的等效參數.等效介質法和有限元仿真計算得到的歸一化等效參數分別如圖5和圖6所示.由圖5(a)可知,通過2種方法得到的增益單元和損耗單元的歸一化等效質量密度實部Re(ρeffb)在計算頻率范圍內幾乎保持不變,并且與分流阻抗無關.由圖5(b)可知,歸一化等效質量密度虛部Im(ρeffb)數值很小,可以認為,等效質量密度ρeff在計算的頻率范圍內為實數.由圖6(a)和圖6(b)可以發現,增益單元和損耗單元的等效彎曲剛度Deff共軛.因此,本文提出的PT對稱梁滿足公式(3)和公式(4)所提出的等效質量密度ρeff和等效彎曲剛度Deff共軛對稱條件.由于電感和分流電阻Rsh的存在,增益單元和損耗單元的等效彎曲剛度Deff與彎曲波頻率相關.等效介質法和有限元仿真獲取的結果總體上相吻合,存在部分差異的原因主要在于等效介質法存在的局限性,如壓電片z方向上電場線性分布的假設等.

為了進一步表征歸一化等效彎曲剛度Deff/Db與分流電阻Rsh、彎曲波頻率的關系,本文采用等效介質法計算不同分流電阻Rsh及不同頻率的彎曲波激勵下損耗單元的歸一化等效彎曲剛度Deff/Db,計算結果分別如圖7和圖8所示.由圖7可知,當彎曲波的頻率遠離530 Hz時,分流電阻Rsh的取值不會影響歸一化等效彎曲剛度實部Re(Deff/Db)和虛部Im(Deff/Db),并且虛部Im(Deff/Db)接近0,表明此時損耗單元幾乎沒有消耗彎曲波能量;當彎曲波頻率為530 Hz附近時,歸一化等效彎曲剛度實部Re(Deff/Db)出現極值,并且增大分流電阻Rsh會使得歸一化等效彎曲剛度虛部Im(Deff/Db)數值減小,從而加強損耗作用.當分流電路的諧振頻率f0=800 Hz時,可以發現圖8(a)和圖8(b)所示結果分別與圖7(a)和圖7(b)所示結果相近,只是歸一化等效彎曲剛度Deff/Db的梯度變大,并且隨著分流電路諧振頻率f0的改變,極值出現的頻率變為850 Hz附近.

由圖5~圖8可知,本文提出的PT對稱梁滿足公式(3)和公式(4)所提出的等效質量密度ρeff和等效彎曲剛度Deff共軛對稱條件.通過改變分流電阻Rsh和諧振頻率f0,增益單元、損耗單元的彎曲剛度Deff會產生較大變化.因此,本文所提出的PT對稱梁具有易于調控的特點.

2.2透射系數、反射系數及散射矩陣

基于傳遞矩陣法和有限元仿真得到的透射系數和反射系數分別如圖9和圖10所示.圖9表示彎曲波分別從左端和右端入射的情況下透射系數的幅值,由圖9可知,2種方法得到的透射系數幅值一致,這是因為PT對稱梁具有互易性,即左透射系數與右透射系數相等.在圖9中的陰影部分內,透射系數的幅值大于1,這是由于PT對稱梁是非保守系統,與外界存在能量交換.圖10(a)和圖10(b)分別表示左反射系數幅值和右反射系數幅值.由圖10可以發現,PT對稱梁對彎曲波的非對稱散射特性,即當彎曲波分別從左、右兩端入射時,2種情況下的反射系數幅值不相等,2種情況下的反射系數都存在零點,該零點對應的彎曲波頻率即為單向無反射點,如圖10(a)中的520.5 Hz及圖10(b)中的511 Hz.由圖9可以發現,在單向無反射點511 Hz和520.5 Hz,透射系數的幅值接近1.

根據PT對稱性,左反射系數rl、右反射系數r1、透射系數t之間滿足:

根據公式(44),可以得到廣義守恒定律[23]:

為了更好地闡述單向無反射特性,本文從有限元仿真的結果中,提取PT對稱梁在2個單向無反射點511 Hz和520.5 Hz下的歸一化彎曲波位移場幅值,結果如圖11所示.由圖11可知,區域I和區域II能直觀地體現PT對稱梁的散射特性,區域I代表左端激勵點到損耗單元左邊界的范圍;區域II代表增益單元右邊界與右端激勵點之間的范圍.當511 Hz 的彎曲波從左端入射時,區域I中有明顯的干涉條紋,表明區域I中入射彎曲波與反射彎曲波形成了駐波,導致區域I中各點的彎曲位移幅值明顯不一致.當511 Hz的彎曲波從右端入射時,區域II中各點的彎曲位移幅值近似相等,表明區域II中只有向右傳播的彎曲波而幾乎不存在反射彎曲波,并且區域I和區域II的彎曲位移幅值相近,表明此時PT對稱梁對右端入射的彎曲波完全透射.當520.5 Hz的彎曲波從左端入射時,區域I中幾乎無反射彎曲波存在,并且彎曲波幾乎完全透射至區域II;而從右端入射時,區域II中反射彎曲波與入射彎曲波形成了駐波,導致區域II中各點之間的彎曲位移幅值有較大差別.

3結論

1)利用壓電分流單元設計一種針對彎曲波的PT對稱梁,其中壓電分流單元連接了正/負電阻并聯電感組成的分流電路,構成損耗單元、增益單元;采用等效介質法計算損耗單元、增益單元的等效彎曲剛度和等效質量密度,并通過有限元仿真對等效介質法進行驗證.結果表明,本文所提出的設計滿足等效質量密度和等效彎曲剛度共軛對稱.

2)采用傳遞矩陣法計算PT對稱梁的透射系數和反射系數,并與有限元仿真的結果進行比較,驗證傳遞矩陣法的正確性.傳遞矩陣法和有限元仿真的結果表明,所提出的PT梁針對彎曲波具有非對稱散射特性,并且存在單向無反射點.

3)通過計算散射矩陣的特征值與特征向量,分析PT對稱梁的奇異點與單向無反射點之間的關系,結果表明,PT對稱梁的奇異點即為單向無反射點.

參考文獻

[1] BENDER C M,BOETTCHER S. Real spectra in non-hermitianhamiltonians Having PT Symmetry [J]. Physical Review Letters,1998,80(24):5243-5246.

[2] GUO A,SALAMO G J,DUCHESNE D,et al. Observation of PT - symmetry breaking in complex optical potentials [J]. Physical Review Letters,2009,103(9):093902.

[3] RUTER C E,MAKRIS K G,EL-GANAINY R,et al. Observation of parity - time symmetry in optics [J]. Nature Physics,2010,6 (3):192-195.

[4] CHRISTENSEN J,WILLATZEN M,VELASCO V R,et al. Paritytime synthetic phononic media [J]. Physical Review Letters,2016,116(20):207601.

[5] LIN Z,RAMEZANI H,EICHELKRAUT T,et al. Unidirectional invisibility induced by PT-symmetric periodic structures [J]. Physical Review Letters,2011,106(21):213901.

[6] KLAIMAN S,GUNTHER U,MOISEYEV N. Visualization of branch points in PT-symmetric waveguides[J]. Physical Review Letters,2008,101(8):080402.

[7] FLEURY R,SOUNAS D,ALU A. An invisible acoustic sensor based on parity-time symmetry [J]. Nature Communications,2015,6(1):5905.

[8] HOU Z,ASSOUAR B. Tunable elastic parity-time symmetric structure based on the shunted piezoelectric materials[J]. Journal of Applied Physics,2018,123(8):085101.

[9]張邦基,黃訓浩,張農,等.局域共振聲子晶體失諧梁的減振特性研究[J].湖南大學學報(自然科學版),2015,42(2):55-59.

ZHANG B J,HUANG X H,ZHANG N,et al. Research on properties of vibration attenuation in phononic crystal detuning beam with local resonant structure [J]. Journal of Hunan University (Natural Sciences),2015,42(2):55-59.(In Chinese)

[10]張浩.壓電分流型聲學超材料的隔聲特性研究[D].長沙:國防科學技術大學,2016:49-53.

ZHANG H. Research on the sound insulation properties of acoustic metamaterials with piezoelectric shunts [D]. Changsha:National University of Defense Technology,2016:49-53. (In Chinese)

[11] HRABAR S,KROIS I,BONIC I,et al. Negative capacitor paves the way to ultra-broadband metamaterials [J]. Applied Physics Letters,2011,99(25):254103.

[12] CHEN Y Y,HUANG G L,SUN C T. Band gap control in an active elastic metamaterial with negative capacitance piezoelectric shunting[J]. Journal of Vibration and Acoustics,2014,136(6):061008.

[13] ZHAN G H,XIAO Y,WEN J H,et al. Ultra-thin smart acoustic metasurface for low-frequency sound insulation [J]. Applied Physics Letters,2016,108(14):141902.

[14] ZHANG X D,CHEN F,CHEN Z S,et al. Membrane-type smart metamaterials for multi-modal sound insulation[J]. The Journal of the Acoustical Society of America,2018,144(6):3514-3524.

[15]羅松南,鄧慶田.帶壓電層圓柱形桿中波的傳播[J].湖南大學學報(自然科學版),2006,33(2):74-77.

LUO S N,DENG Q T. Wave propagation in the layered piezoelectric cylindrical bars[J]. Journal of Hunan University(Natural Sciences),2006,33(2):74-77(. In Chinese)

[16] WANG G,CHEN S B. Large low-frequency vibration attenuation induced by arrays of piezoelectric patches shunted with amplifierresonator feedback circuits [J]. Smart Materials and Structures,2016,25(1):015004.

[17] CHEN Y Y,HU G K,HUANG G L. A hybrid elastic metamaterial with negative mass density and tunable bending stiffness[J]. Journal of the Mechanics and Physics of Solids,2017,105:179-198.

[18] COLLET M,OUISSE M,ICHCHOU M N. Structural energy flow optimization through adaptive shunted piezoelectric metacompos- ites[J]. Journal of Intelligent Material Systems and Structures,2012,23(15):1661-1677.

[19]朱濱.彈性力學[M].合肥:中國科學技術大學出版社,2008:276-281.

ZHU B. Elastic mechanics[M]. Hefei:Press of University of Science and Technology of China,2008:276-281.(In Chinese)

[20] WANG G,CHENG J Q,CHEN J W,et al. Multi-resonant piezoelectric shunting induced by digital controllers for subwavelength elastic wave attenuation in smart metamaterial[J]. Smart Materials and Structures,2017,26(2):025031.

[21] LI D T,HUANG S B,CHENG Y,et al. Compact asymmetric sound absorber at the exceptional point [J]. Science China Physics,Mechanics & Astronomy,2021,64(4):1-7.

[22]楊小黎.結構聲學材料反向聲散射特性研究[D].哈爾濱:哈爾濱工程大學,2020:12-14.

YANG X L. Research on the characteristics of acoustic backscat-tering of structural materials[D]. Harbin:Harbin Engineering University,2020:12-14.(In Chinese)

[23] GE L,CHONG Y D,STONE A D. Conservation relations and anisotropic transmission resonances in one-dimensional PT- symmetric photonic heterostructures[J]. Physical Review A,2012,85(2):023802.

[24] CHONG Y D,GE L,STONE A D. PT-symmetry breaking and laser-absorber modes in optical scattering systems[J]. Physical Review Letters,2011,106(9):093902.

91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合