999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

一種平衡銜鐵受話器數值仿真分析方法

2022-12-05 07:49:18溫周斌
聲學技術 2022年5期
關鍵詞:磁場分析模型

蔣 鈺,陸 曉,溫周斌

(1.中國科學院聲學研究所東海研究站,上海 201815;2.中國科學院聲學研究所,北京 100049;3.中國科學院大學,北京 100190;4.浙江中科電聲研發中心,浙江 嘉善 314115)

0 引言

平衡銜鐵受話器(Balanced Armature Receiver,BAR),具有尺寸小、電聲轉換效率高和靈敏度高等特點[1-2],已被廣泛用于助聽設備、入耳式耳機和各種軍用頭戴式耳機等產品中[3]。近年來,隨著消費者對聲品質要求提高,平衡銜鐵受話器的應用更為廣泛。

數值仿真分析是一種通過計算機進行數值計算求解一組滿足一定邊界條件的微分方程,并得到其近似解的方法[4],已廣泛應用于動圈式揚聲器單體及系統的主要特性的仿真分析[5-8]。近10年來,國內外學者圍繞平衡銜鐵受話器的非線性、聲壓級和失真等主要特性進行了深入研究。Jensen等基于集總參數模型深入研究了平衡銜鐵受話器的非線性和失真等特性[9]。基于集總參數模型的仿真分析方法概念清晰、方法簡便,但所得到的仿真結果在中高頻段精度不夠高[10]。Jensen還采用數值仿真分析方法研究了不同平衡銜鐵磁導率對應的平衡銜鐵上的磁通密度分布[11]。Hu等創建了平衡銜鐵受話器的振動系統與聲場的數值仿真分析模型,驗證了其所提出的多模態模型的有效性[10]。Ziolkowski等將平衡銜鐵受話器磁路簡化為2D模型,通過耦合電磁場與結構,得到了平衡銜鐵的位移[12]。Xu等提出了平衡銜鐵受話器電磁場的瞬態分析和基于模態疊加法的平衡銜鐵受話器結構與聲場的頻域分析相結合的方法,利用前者得到的非線性參數和后者得到的力與聲壓級的傳遞函數,再通過數值迭代最終得到平衡銜鐵受話器輸出的聲壓級[13-14]。

本文首先介紹平衡銜鐵受話器及其集總參數模型,分析指出平衡銜鐵運動會對揚聲器特性帶來重大影響;然后提出一種基于COMSOL Multiphysics數值仿真分析軟件的平衡銜鐵受話器聲壓級和阻抗的仿真分析方法。該方法首先對平衡銜鐵受話器的磁場特性進行穩態分析,計算得到可隨平衡銜鐵運動變化的驅動力和反電動勢;然后對電磁場、結構和聲場耦合模型進行頻域分析,最終計算得到平衡銜鐵受話器的聲壓級和阻抗。

1 平衡銜鐵受話器及其集總參數模型

圖1是一款平衡銜鐵受話器的半剖幾何模型[1]。由圖1可知,平衡銜鐵受話器由線圈(Coil)、傳動桿(Drive pin)、振膜(Membrane)、永磁體(Magnets)、磁軛(Magnet house)、平衡銜鐵(Armature)、前蓋(Cover)、外殼(Case)和出聲孔(Sound outlet)等組成。平衡銜鐵一端被固定在外殼上,另一端懸空。平衡銜鐵通過驅動桿與振膜連接。振膜被粘于上蓋,其中三邊有折環,另外一邊是金屬。平衡銜鐵和振膜的振動都是類似于懸臂梁的振動。振膜和上蓋之間的空氣域定義為前腔,振膜另一面的空氣域為后腔。

圖1 平衡銜鐵受話器半剖幾何模型[1]Fig.1 Section model of BAR[1]

當平衡銜鐵受話器音圈中無電流通過時,平衡銜鐵在磁場中處于受力平衡狀態,因而穩定在兩個永磁體中間。當線圈中有電流通過時,在線圈周圍除了穩定的磁場之外,還有因電流變化所感應出的擾動磁場。平衡銜鐵會被擾動磁場磁化,它的磁化方向也會隨著流經線圈中電流方向的改變而改變。因此,平衡銜鐵所受到的力將不斷變化,這個力會使得平衡銜鐵上下往復運動。平衡銜鐵往復運動將通過傳動桿帶動振膜上下振動,而振膜的上下振動又會使得前腔的空氣壓縮或膨脹,從而產生聲波,聲波再從出聲孔傳遞出去。

建立平衡銜鐵受話器的集總參數模型,可分析得到諸多重要參數,進而還可分析得到這些重要參數對揚聲器特性的影響。圖2(a)給出了一個平衡銜鐵受話器的簡化的磁回路。將磁阻、磁場力和磁通量分別等效為電阻Rg、電壓源FM和電流Φ,可建立該磁回路的集總參數模型,圖2(b)是圖2(a)的等效電路圖[10]。利用該等效電路圖,可以得到平衡銜鐵中的磁通量和上下磁隙中的磁通量,進而計算得到線圈上所受到的反電動勢和平衡銜鐵所受的力。

圖2 平衡銜鐵受話器簡化的磁回路及其等效電路[9]Fig.2 Simplified magnetic circuit and equivalent circuit of BAR[9]

反電動勢Uback可表示為[11]

平衡銜鐵所受的力,該力F可表示為[11]

式中,kΦ(x)為非線性補償剛度系數,Tme(x)和Tme,d(x,i)分別為非線性傳遞系數和失真傳遞系數。由式(2)可知,當揚聲器工作時,平衡銜鐵所受的力不僅與流經線圈的電流有關,還與平衡銜鐵的位移相關。與位移相關的力可稱為動生磁場力。由于平衡銜鐵是導磁材料,其運動會導致揚聲器磁場特性變化,磁場特性的變化會導致平衡銜鐵受到動生磁場力的作用。

平衡銜鐵所受的力隨位移變化的特性會對揚聲器性能產生重要影響。平衡銜鐵的有效剛度是其機械剛度與非線性補償剛度系數的差值,非線性補償剛度系數增加,會減小平衡銜鐵的有效剛度,從而減小揚聲器的共振頻率;該系數的減小則會使振膜振幅增大,聲壓也會增大[11]。而該系數的增加會增加揚聲器的不穩定性,這是因為有效剛度減小會使得平衡銜鐵可能無法回到平衡位置。

2 數值仿真分析方法

由前文平衡銜鐵受話器集總參數模型分析可知,平衡銜鐵的運動會導致揚聲器的磁場特性發生變化,在線圈上產生一個額外的反電動勢,并使得平衡銜鐵受到一個額外的磁場力作用。采用數值仿真分析方法仿真計算揚聲器的聲場特性時一般選用頻域分析,直接的頻域分析則無法考慮平衡銜鐵運動對磁場特性的影響。

為此,本文提出一種基于COMSOL Multiphysics數值仿真分析軟件的平衡銜鐵受話器數值仿真分析方法,分為兩個步驟:

(1)對平衡銜鐵受話器的磁場特性進行穩態分析,并利用動網格技術,計算得到可隨平衡銜鐵運動變化的作用在線圈上的反電動勢(速度電動勢)和施加在平衡銜鐵上的動生磁場力。

(2)創建平衡銜鐵受話器電磁場、結構和聲場的耦合分析模型,并利用前文穩態求解得到的可隨平衡銜鐵速度變化的速度電動勢和可隨平衡銜鐵位移變化的動生磁場力,再進行頻率分析,最終計算得到揚聲器的聲壓級和阻抗。

無論是穩態分析,還是頻域分析,數值仿真分析主要包括建立幾何模型、定義材料參數、設置邊界條件、網格劃分、計算和結果后處理等步驟。

由于影響仿真分析結果準確性的主要因素包括幾何模型、材料參數和仿真分析方法等三個部分,因此,下文以一款平衡銜鐵受話器為例,首先簡要介紹它的幾何模型和材料參數,然后再詳細介紹它的聲壓級和阻抗的數值仿真分析方法。

2.1 幾何模型和材料參數

圖3給出了一款平衡銜鐵受話器的幾何模型。平衡銜鐵受話器具有幾何對稱性,因此只需建立其1/2幾何模型即可。為了提高計算效率,模型中將帶有花紋的金屬球頂簡化為平面,磁隙間距為0.1mm。

圖3 一款平衡銜鐵受話器的幾何模型(1/2模型)Fig.3 Half geometry model of BAR

表1和表2分別給出了該揚聲器各個部件的磁性材料參數和力學材料參數。由于在模型中簡化了振膜的花紋,所以需要適當提高振膜的楊氏模量,使得簡化后的振膜剛度可以等效包含花紋的振膜剛度。圖4給出了平衡銜鐵受話器磁軛(軟磁材料)的BH曲線。其中,B為磁通密度,H為磁場強度,線圈匝數為350匝。

表1 平衡銜鐵受話器部件的磁性材料參數Table 1 Magnet material properties of BAR

表2 平衡銜鐵受話器部件的力學材料參數Table 2 Mechanical material properties of BAR

圖4 磁軛(軟磁材料)的BH曲線Fig.4 BH curve of soft iron

假設所有材料是各向同性的。對于磁性材料,需要給出磁通密度B和磁場強度H之間的本構關系,用剩磁Br和回復磁導率μrec定義永磁體本構關系[13-14]:

式中:μ0為真空磁導率。磁軛的磁通密度和磁場強度之間的關系是非線性的,需用完整的BH曲線來定義磁軛的本構關系。平衡銜鐵的磁通密度和磁場強度之間的關系也是非線性的,但本文將其簡化,僅用相對磁導率μr(常數值)來定義平衡銜鐵的本構關系[13-14]:

對于平衡銜鐵、傳動桿和振膜等振動系統部件,均采用線彈性材料模型定義它們的應力應變關系。

2.2 平衡銜鐵運動對揚聲器特性的影響

為了仿真分析平衡銜鐵運動對磁場特性的影響,平衡銜鐵受話器的幾何模型只需包含磁路和平衡銜鐵,忽略振膜和傳動桿等其他部件;假設金屬外殼有一定的磁屏蔽作用,因而忽略外殼以外的空氣域,僅建立外殼內的空氣域。材料參數也只需定義磁路系統各個部件的磁性材料參數和平衡銜鐵的力學材料參數。在平衡銜鐵與外殼相接面上施加固定約束條件,在對稱面上,通過定義磁矢勢的面內分量為0和位移的法向分量為0,來分別設定磁場和結構的對稱條件。在劃分網格時,需對平衡銜鐵適當加密網格。

在位于磁隙區域的平衡銜鐵表面上添加力載荷(如圖5所示),平衡銜鐵受力后會發生形變,可根據力的平衡方程計算得到平衡銜鐵的位移大小。平衡銜鐵作為導磁材料,其形變會改變磁場分布。為了模擬平衡銜鐵運動(形變)對磁場分布的影響,模型中需要使用動網格技術,即網格將隨著平衡銜鐵的運動(形變)而變化。基于麥克斯韋方程組計算得到磁通密度分布。在平衡銜鐵表面上加載不同大小的力,便可得到平衡銜鐵運動到不同位置時的磁通密度分布。

圖5 平衡銜鐵的加載面圖Fig.5 The loading surface diagram of the balanced armature

圖6(a)和圖6(b)分別給出了當平衡銜鐵處于初始位置時及其從初始位置向上移動至磁隙間距的1/10位置時線圈所包圍區域的磁通密度分布圖。由圖6可以看出,平衡銜鐵位置(位移)的變化,會導致該區域磁通密度發生明顯改變。磁通密度改變,會影響作用在線圈上感應電動勢和平衡銜鐵上所受到的力。

圖6 平衡銜鐵處于不同位置時的磁通密度分布圖Fig.6 Diagrams of the magnetic flux density distribution when the balanced armature is at different positions

2.2.1 作用在線圈上的速度電動勢

由式(1)可知,作用在線圈上的速度電動勢Uv為

根據電磁感應定律,速度電動勢Uv是磁通量Φ(x)對時間t的導數:

在平衡銜鐵的位移比較小的情況下,假設線圈中的磁通量Φ(x)與平衡銜鐵位移x呈線性關系;根據式(5)和式(6),可以將Φ(x)與x的關系表示為

式中,系數Tme和a均為常數。當位移x較小時,近似認為非線性傳遞系數Tme(x)是一個不隨x改變的常數Tme。

從前述穩態求解得到的仿真結果中選取幾個平衡銜鐵的位置,并根據平衡銜鐵處于這些位置時的磁場分布計算線圈所包圍面積的磁通量,再通過式(7)進行函數擬合,可以得到系數Tme和a。

將式(7)代入式(6)可得:

2.2.2 平衡銜鐵上所受到的動生磁場力

平衡銜鐵表面的電磁力密度矢量分量f為[14]:

式中,n是表面單位法向矢量,H2=H2x+H2y+H2z。

圖7給出了平衡銜鐵在平衡位置時其所受的表面電磁力密度矢量分量的分布圖,由圖7可知,該作用力主要集中在磁隙區域對應的表面上。

圖7 平衡銜鐵上表面應力分布Fig.7 Surface stress distribution of balanced armature

當平衡銜鐵運動時,它與上下永磁體之間的距離會發生變化,距離越小,磁通密度越大,靠近永磁體的平衡銜鐵表面應力會增大,遠離永磁體的平衡銜鐵表面應力會減小。

由式(9)可計算得到平衡銜鐵表面應力(表面上每個節點的應力)隨著位移變化的變化量,再對該表面應力的變化量求積分,即可得到(整個)平衡銜鐵上所受到的力的變化量,也即動生磁場力FΦ:

式中:Δf是應力的變化量,S表示平衡銜鐵表面。

由式(2)可知,平衡銜鐵上受到的動生磁場力與位移相關,當平衡銜鐵的位移非常小時,假設動生磁場力FΦ與位移x呈線性關系,即:

式中:系數kΦ為常數。當位移x較小時,近似認為非線性補償剛度系數kΦ(x)是一個不隨x改變的常數kΦ。

同理,從前述穩態求解得到的仿真結果選取幾個平衡銜鐵的位置,并根據式(10)計算對應的動生磁場力,再通過式(11)進行函數擬合,計算得到系數kΦ。

2.3 磁路、振動系統和聲場的耦合分析

創建包含磁路和振動系統所有部件的平衡銜鐵受話器的幾何模型,因為它具有幾何對稱性,所以只建立其1/2幾何模型即可。由于本案例的平衡銜鐵受話器是作為受話器使用,因此,在仿真模型中,其聲輸出方式為711耦合器(壓力場環境)。為提高計算效率,可采用特定的“阻抗邊界”條件替代“711耦合器聲腔(模型)”的簡化方法,因而在仿真模型中并不需要建立耦合器聲腔模型。圖8(a)是包含耦合器聲腔的全模型,圖8(b)則是簡化后的幾何模型。對于圖8(b),只需要包含平衡銜鐵受話器和小部分聲腔,這樣可以大幅降低計算量。采用表1和表2中的材料參數定義所創建模型中揚聲器各個部件的材料參數。

圖8 平衡銜鐵受話器幾何模型(含711耦合器)Fig.8 Geometry model of BAR(including 711 coupler)

在平衡銜鐵與外殼相接面上施加固定約束。在對稱面上,通過定義磁矢勢的面內分量為0、位移的法向分量為0和加速度的法向分量為0,分別設定磁場、結構和聲場的對稱條件。空氣域中除了與振動系統耦合面以外的邊界面,全部假設為全反射邊界。在小部分聲腔與711耦合器的連接面上設定聲阻抗,可等效711耦合器空氣域對邊界面上聲壓產生的影響,這個聲阻抗Zin可以預先利用圖8(a)全模型仿真分析得到,計算公式為

式中:Pin為輸入端面處的聲壓有效值,Zin是與輸入端面處空氣體積流速Qin相關的聲阻抗。合理劃分網格,以保證計算頻率范圍內一個波長至少包含4~5個網格。

平衡銜鐵受話器的磁路、振動系統和聲場是相互影響的,在模型中需要正確定義它們的耦合作用。

對于磁場和結構這兩個物理場,實際加載在線圈上的電壓U是由外部加載電壓U0和反電動勢Uback兩部分組成,而反電動勢包含由于電流變化帶來的反電動勢Ui和由于平衡銜鐵運動帶來的反電動勢(速度電動勢)Uv,即:

平衡銜鐵所受到的力F是由線圈電流產生的感生磁場力Fi和因其自身運動而產生的動生磁場力FΦ兩部分組成的,即:

對于結構和聲場這兩個物理場,振膜表面的加速度即為聲源。同時,作用在振膜表面上的空氣壓強會使得振動系統受到額外的表面應力。

在線圈兩端加載諧波擾動的電壓U0,電壓有效值為0.1 V。根據電磁感應定律,變化的電場會感應出磁場,磁路系統中除了原本的靜態磁場以外,會再疊加一個諧波擾動的磁場。根據麥克斯韋方程組求解得到電場分布和擾動的磁場分布。根據式(12)可計算得到平衡銜鐵所受到的感生磁場力Fi。平衡銜鐵受力產生振動,根據結構運動方程[15]求解得到振動幅值。根據振動系統與聲場的邊界面上振幅(或加速度),可計算得到聲場分布。一般空氣域的聲場分析建立在理想流體媒質的波動方程[15]之上。而對于平衡銜鐵受話器,在它的內部存在一些狹窄區域,不能忽略聲波在狹窄區域中的粘滯特性。

如前文所述,為提高計算效率,幾何模型中并沒有包含711耦合器空氣域,因此無法直接獲取麥克風面處的聲壓級。利用已建立的幾何模型聲腔邊界面(也即711耦合器輸入端面)上的空氣體積流速Qin和711耦合器空氣域的轉移阻抗Ztrans,可計算得到麥克風表面位置的聲壓pmic:

圖9給出了該平衡銜鐵受話器在頻率為1 kHz信號激勵下的聲壓級分布。由圖9可知,在振膜與上蓋之間的腔體內的聲壓最高,聲波通過出聲孔向外(711耦合器)傳遞。

圖9 平衡銜鐵受話器中的聲壓級分布圖(@1 kHz)Fig.9 Sound pressure level(SPL)of BAR(@1 kHz)

3 仿真分析結果

通過本文提出的數值仿真分析方法可以較快地仿真計算得到平衡銜鐵受話器的聲壓級和阻抗曲線。與此同時,作者還采用B&K Pulse系統測量了所研究的平衡銜鐵受話器的聲壓級和阻抗曲線,為保證測量與仿真分析條件的一致,測量是在IEC711耦合器(愛宏AWA6162耳模擬器)壓力場條件下進行的,符合IEC60318-4標準[16]。被測揚聲器通過夾具安裝在耳模擬器上,揚聲器和夾具之間使用少量橡皮泥以防止泄露,圖10(a)是被測平衡銜鐵受話器,外尺寸為6.3 mm×4.3 mm×3.0 mm。圖10(b)是被測揚聲器通過夾具安裝在耳模擬器上的照片。圖10(c)給出了測量示意圖,揚聲器出聲孔平面應與耳模擬器參考平面齊平。

圖10 平衡銜鐵受話器的測量Fig.10 Measurement of BAR

圖11(a)和圖11(b)分別給出了平衡銜鐵受話器的聲壓級和阻抗的仿真分析結果和測量結果。由圖11可知,在4 kHz以下,仿真分析結果和測量結果基本吻合;在高頻段,兩者趨勢一致,數值存在一定差異。

圖11 仿真結果和測量結果Fig.11 Comparison between measurement and simulation

由圖11(a)可知,在低頻段阻抗值基本與直流阻抗一致;而在高頻段阻抗值隨頻率增大而增大。這是因為由電流變化產生的感應電動勢增大的緣故。感應電動勢越大,電流越小,阻抗就越大。2 kHz頻率附近的峰以及4、7 kHz頻率附近的峰都與速度電動勢有關,由于振動系統在這些頻率附近發生共振(或諧振),平衡銜鐵的振幅相對較大,平衡銜鐵運動引起的位置改變所帶來的速度電動勢也增大。同樣,速度電動勢越大,電流越小,阻抗也就越大。

由圖11(b)可知,在低頻段聲壓級基本為一穩定值,當頻率大于1 kHz,聲壓級曲線出現多個諧振峰。低頻聲壓級與振動系統的有效剛度密切相關,而動生剛度是有效剛度的重要組成部分。動生剛度越大,有效剛度越小,振幅則越大,低頻聲壓級也就越高。2 kHz和7 kHz附近的諧振峰也與振動系統的共振(或諧振)相關,這些峰所對應的頻率高低與動生剛度有關。動生剛度越大,有效剛度則越小,這些峰對應的頻率也會越低。而14 kHz頻率附近的共振峰則與此無關,它是由于711耦合腔的共振所引起的。

本文所提方法較為準確地仿真了可隨平衡銜鐵速度變化的速度電動勢和可隨平衡銜鐵位移變化的動生磁場力。分析得到了阻抗共振峰和諧振峰,最終較為準確地仿真得到了低頻聲壓級和聲壓級頻響曲線上的各個峰。

4 結論

本文基于平衡銜鐵受話器的集總參數模型分析,指出平衡銜鐵運動會產生額外的反電動勢和磁場力,提出了一種分步進行的平衡銜鐵受話器的數值仿真方法。首先對平衡銜鐵受話器的磁場特性進行穩態分析,計算得到可隨平衡銜鐵運動而變化的磁場力和反電動勢;然后,建立電磁場、結構和聲場耦合分析模型,并利用前面的穩態求解方法得到的可隨平衡銜鐵運動變化的速度電動勢和動生磁場力,再進行頻域分析,最終計算得到平衡銜鐵受話器的聲壓級和阻抗曲線。仿真結果和測量結果基本一致。

平衡銜鐵受話器的數值仿真分析方法還可以在以下幾個方面進一步優化和完善:(1)深入研究平衡銜鐵材料的磁特性,以建立更準確的非線性模型;(2)細致分析小體積聲腔中的空氣壓強對振動系統的影響;(3)非線性和失真特性機理和高效仿真分析方法研究。

猜你喜歡
磁場分析模型
一半模型
西安的“磁場”
當代陜西(2022年6期)2022-04-19 12:11:54
為什么地球有磁場呢
隱蔽失效適航要求符合性驗證分析
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
磁場的性質和描述檢測題
電力系統及其自動化發展趨勢分析
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产精品嫩草影院av| 国产精品妖精视频| 亚洲第一页在线观看| 五月天福利视频| 欧美va亚洲va香蕉在线| 久久久久青草线综合超碰| 亚洲AⅤ波多系列中文字幕| 女人18毛片一级毛片在线 | 国产一级毛片网站| 亚洲欧洲日韩久久狠狠爱 | аv天堂最新中文在线| 精品欧美一区二区三区久久久| 欧美综合区自拍亚洲综合绿色 | 欧美亚洲一区二区三区导航| 一级福利视频| 成人精品在线观看| 国产在线视频福利资源站| 国产人妖视频一区在线观看| 国产成人精品2021欧美日韩| 国产欧美视频一区二区三区| 国模私拍一区二区| 1024国产在线| 精品国产免费观看| 中文字幕人成乱码熟女免费| 波多野结衣一区二区三区四区| 欧美国产精品拍自| 国产亚洲成AⅤ人片在线观看| 三级国产在线观看| 欧美性爱精品一区二区三区| 免费国产一级 片内射老| 亚洲无线国产观看| 久久久精品无码一区二区三区| av在线无码浏览| 欧美三级视频网站| 人人看人人鲁狠狠高清| 亚洲av综合网| 亚洲精品综合一二三区在线| 亚洲精品爱草草视频在线| 人人妻人人澡人人爽欧美一区| 青青热久免费精品视频6| 青青操视频在线| 亚洲激情99| 亚洲无码一区在线观看| 欧洲日本亚洲中文字幕| 在线日韩日本国产亚洲| 99久久免费精品特色大片| 国产精品国产三级国产专业不| 亚洲AV无码一二区三区在线播放| 福利姬国产精品一区在线| 欧美另类图片视频无弹跳第一页| 成人毛片在线播放| 欧美日韩精品一区二区在线线| 亚洲成网777777国产精品| 麻豆AV网站免费进入| 久久久国产精品免费视频| 精品国产中文一级毛片在线看 | 精品一区国产精品| 午夜人性色福利无码视频在线观看| 久久久久青草大香线综合精品| 国产成人你懂的在线观看| 亚洲av无码久久无遮挡| 99久久精品国产精品亚洲 | 九九九九热精品视频| 一个色综合久久| 国产av剧情无码精品色午夜| 国产午夜人做人免费视频中文| 国产成人成人一区二区| 白浆免费视频国产精品视频| 黄色网址手机国内免费在线观看| 国产精品999在线| 国产97公开成人免费视频| 欧美午夜理伦三级在线观看 | 女人毛片a级大学毛片免费| 精品国产成人av免费| 久草视频中文| 精品久久高清| 国产精品真实对白精彩久久| 91免费国产高清观看| 99久久性生片| 亚洲AV成人一区二区三区AV| 一级一毛片a级毛片| 一级一级一片免费|