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

基于人工神經網絡的非結構網格尺度控制方法1)

2021-12-02 02:31:34王年華常興華張來平鄧小剛
力學學報 2021年10期
關鍵詞:背景效率方法

王年華 , 魯 鵬 * 常興華 張來平 鄧小剛

* (中國空氣動力研究與發展中心空氣動力學國家重點實驗室,四川綿陽 621000)

? (西南科技大學信息工程學院,四川綿陽 621010)

** (重慶文理學院智能制造工程學院,重慶 402160)

?? (國防科技創新研究院無人系統技術研究中心,北京 100071)

*** (軍事科學院,北京 100091)

引言

網格生成是計算流體力學(computational fluid dynamics,CFD)數值計算的第一步,也是未來CFD六大重要研究領域之一[1-2].在現代CFD 應用過程中,自動生成復雜構型的高質量網格(包括網格自適應)依然是一個重大挑戰性問題.自動化程度和網格質量是網格生成過程中最重要的兩個問題[3-4].據統計,網格生成通常占據整個計算周期大約60%的人力時間,高度自動化的網格生成方法無疑可以很大程度節約CFD 計算周期內的人工成本.

網格空間尺度分布控制在網格生成中至關重要,對于網格質量和流動求解精度影響較大.通常希望在幾何曲率大、流動梯度大等重點區域加密網格,而在非重點區域則希望網格盡可能稀疏并均勻過渡.傳統的非結構網格尺度控制方法主要有:函數指定法、插值類方法、背景網格法和根據流場特征進行自適應等.

函數指定法適用于簡單問題的全場網格尺度控制,比如可以指定線性函數控制翼型的網格尺度分布;或用于局部網格尺度控制,比如可以指定幾何級數或指數函數對局部網格進行加密[5].

采用背景網格法[6-9]進行網格尺度控制,需要先用規則的矩形結構網格或采用較稀疏的非結構網格覆蓋全計算域,在計算域內分布一定的“點源”“線源”或者“面源”等局部網格分布控制參數,將這些源視作離散的“熱源”,求解熱傳導方程(泊松方程),得到的穩態解即為全計算域的尺度控制參數分布[6].或者根據局部幾何特征(曲率、狹縫、窄邊等信息)確定局部網格尺度,再在背景網格上求解網格尺度滿足的梯度限制方程,將局部網格尺度光滑到全場,得到背景網格上的網格尺度分布[7-9].在生成網格時,根據控制空間中某點在矩形背景網格中所處的單元,通過背景網格的尺度分布插值得到該點處的網格尺度.

除此之外,還有一些學者提出一些其他形式的背景網格法.如Deister 等[10]提出由最大最小尺度及最大曲率角對幾何進行柵格化,計算得到局部網格尺度,并存儲在自適應的背景笛卡爾網格上.Quadros等[11-12]提出采用幾何體離散骨架的幾何臨近信息、特征尺寸、邊界曲率來測量幾何復雜度,并根據幾何復雜度生成點源,根據點源確定網格尺度分布,最終將網格尺度存儲在叉樹結構的笛卡爾背景網格上.Ruiz-Girones 等[13]提出通過在背景網格上求解一種新的非線性方程來控制四邊形網格尺度,等等.

徑向基函數(radial basis function,RBF)[14-15]可用于數據插值,在已知邊界上的網格尺度后,可用RBF 方法將邊界網格尺度插值到內場,采用貪婪算法還可以一定程度提高RBF 插值的效率,是一種簡單高效的插值類尺度控制方法.

根據流動特征物理量的梯度量等判據進行網格自適應[16-17]也是控制網格尺度分布的一種有效方法,能夠根據流場變化控制網格疏密,在梯度大的區域生成更密的網格,能夠更精細地捕捉流動特征,具有更高的計算精度和適應性,是來重要發展方向之一[18].

近年來,基于人工神經網絡的深度學習方法在工業社會甚至流體力學領域得到廣泛研究和應用[19-21].大數據驅動的人工智能方法成為除理論分析、數值計算和實驗技術以外一種新的研究范式,為各個領域帶來了新的研究思路和方法.

在網格生成領域,經過多年的工程實踐,已經積累了大量各種類型的網格數據,這些數據包含了網格生成規則及技術人員在網格生成方面的知識和經驗,是天然的機器學習訓練樣本.通過機器學習對網格生成規則進行學習,可以簡化傳統算法,提高網格生成效率[22].而網格尺度分布的控制也需要技術人員根據對流動問題的分析和經驗合理確定,采用機器學習方法對網格尺度分布進行控制,有望減少人工工作量和對人工經驗的依賴.

本文從網格質量、效率、靈活性和自動化程度4 個方面綜合分析各類網格尺度控制方法的優缺點.為了克服傳統背景網格插值法效率低、自動化程度不高等方面的不足,本文從效率和自動化程度角度提出兩種網格尺度控制方法,首先將RBF 插值方法應用于網格尺度控制,采用貪婪算法對RBF 插值參考點序列進行精簡,實現高效的RBF 網格尺度分布控制方法.進一步將提出一種采用人工神經網絡進行非結構網格尺度控制的方法,通過引入相對壁面距離和相對網格尺度,初步確定合理的神經網絡輸入輸出參數,建立人工神經網絡訓練模型,采用商業軟件生成二維圓柱和二維翼型非結構三角形網格作為訓練樣本,通過訓練和學習建立起相對壁面距離和相對網格尺度之間的映射關系,進而實現不同密度的二維圓柱、不同二維翼型在不同遠場大小情況下的網格尺度分布控制.

1 基于背景網格法的尺度控制方法回顧

背景網格法可采用規則的笛卡爾直角結構網格、非規則的非結構網格或者自適應笛卡爾網格[9],各種類型的背景網格各有優缺點,由于規則結構網格求解和插值效率較高而得到廣泛使用.

在背景網格上布置“點源”、“線源”或者“面源”等局部網格分布控制參數,通過求解熱傳導方程,得到穩態解即為全流場的尺度控制參數分布[5].該方法能夠生成分布均勻的網格,且能更靈活地考慮流動局部特征對網格分布的影響.

1.1 背景網格法

以二維問題為例,網格尺度S滿足如式(1)所示的穩態熱傳導方程(泊松方程),以確保尺度源項的作用在物理空間內光滑分布.

式中G為源項,定義如式(2)所示

式中下標“i”,“j”代表背景網格節點,N為熱源總數,Ψn為第n個源的強度因子,函數In及Jn分別定義如下

式中S和f為在各源處的尺度,r為背景網格節點到源的距離,|ln|為線源長度.

文獻[5]還提出了考慮方向性的網格密度控制方法,可以通過改變源的強度函數來實現,具體可參考文獻,本文不再贅述.

因此,采用背景網格法進行網格尺度控制的步驟為:

(1)在矩形背景網格上采用中心差分離散泊松方程;

(2)采用Gauss-Seidel 超松弛迭代求解離散的泊松方程,得到背景網格上的網格尺度分布;

(3)在網格生成過程中,根據當地位置在背景網格中插值得到當地網格尺度,用于控制網格生成過程.

1.2 背景網格法網格生成實例

本節采用二維圓柱、NACA0012 翼型和30P30N多段翼型作為考核算例,對背景網格法進行實例測試.本文算例中采用的三角形網格生成算法為作者發展的基于ANN (artificial neural network)的陣面推進法[22],該方法在傳統陣面推進法的基礎上,通過引入ANN 進行生成模式判斷和新點預測,減少了相交性判斷次數,網格生成效率提高30%,其網格生成步驟可簡要概括為:

(1)將幾何邊界離散成初始陣元;

(2)從最小陣面出發,自動選擇網格模板點,人工神經網絡根據網格模板判斷生成模式并預測新點坐標;

(3)根據生成模式、新點坐標和局部網格尺度生成新網格單元;

(4)判斷新單元是否合適,合適則更新數據結構;

(5)回到步驟2,直至所有面變成非活躍面,整個計算域被網格填滿.

圖1~圖3 給出了3 個算例點源設置示意圖和生成的網格,圖中藍色圓點即為點源所在位置.3 個算例人工設置的點源數量分別為12,24 和44.圖中背景網格僅作為示意,實際背景網格節點數量需要根據邊界網格尺度進行調整,以確保對網格尺度分布場的分辨率和插值精度.比如對于在前后緣網格尺度較小的翼型,背景網格必須足夠密才可以有效反映出空間網格尺度的變化.

圖1 圓柱算例點源設置及網格生成情況Fig.1 Nodal source settings and corresponding triangular mesh over a 2D cylinder

圖2 NACA0012 算例點源設置及網格生成情況Fig.2 Nodal source settings and corresponding triangular mesh over NACA0012 airfoil

圖3 30P30N 算例點源設置及網格生成情況Fig.3 Nodal source settings and corresponding triangular mesh over 30P30N airfoil

在本文三個算例中,背景網格規模分別為51 × 51,301 × 301 和401 × 401.由于是在背景網格上迭代求解泊松方程,因此背景網格的網格數量直接決定了迭代求解的效率.

圖中結果顯示,在人工設置合適的點源參數后,可以在背景網格上得到恰當的網格尺度分布.

要說明的是本文暫未采用線源,實際上線源可以看作按線段排列的具有一定強度分布的點源集合,因此在處理本文的簡單問題時,只采用點源進行網格尺度的控制,用以說明背景網格法的優缺點.

2 基于RBF 方法的尺度控制方法

徑向基函數是一種常用的插值函數,常用于網格變形的插值[14-15,23-24],如圖4 所示即為用NACA-0012 翼型的S 型啟動來模擬魚的游動過程,采用RBF 插值方法將物面的變形插值到空間來生成動網格.RBF 插值方法也可以用于網格尺度的插值控制,在已知邊界網格尺度的情況下,求得插值系數矩陣,可將邊界尺度插值到整個計算域.

圖4 RBF 網格變形方法Fig.4 Mesh deformation controlled by RBF interpolation

2.1 基于貪婪算法的RBF 插值方法

RBF 插值方法是比較成熟的插值方法,對于一個變量場,如網格變形情況下的位移場,或網格尺度控制情況下的網格密度場,可以用RBF 插值公式表示為

其中,N為參考點數目,f(r) 為某待求點的函數值(網格尺度),r為待求點的位置矢量,ri為參考點的位置矢量,φ 為RBF 基函數,‖r-ri‖ 為參考點與待求點之間的歐氏距離.wi為第i個參考點的權重系數.

以網格尺度分布控制為例,權重系數滿足

式中Sp為各個參考點的網格尺度,參考點通常為邊界點,也即邊界網格尺度.基函數的類型包括全域型基函數和緊支型基函數,具體可以參考文獻[14-15,23-24],本文不再贅述.

權重系數的求解通常需要求解系數矩陣的逆,在參考點數目較大時,會導致求解效率較低,因此可以引入貪婪算法對參考點進行精簡.

基于貪婪算法的RBF 方法的步驟可簡要歸納為:

(1)初始參考點集為空集,隨機取一邊界點作為初始參考點;

(2)根據參考點,采用RBF 插值計算其他邊界點處的網格尺度;

(3)由于參考點數目太少,RBF 插值得到的邊界網格尺度與給定值存在一定誤差,找出誤差最大的點;

(4)若誤差最大的點不為已有參考點,則將該點做為新參考點加入參考點序列,否則重新任選一點加入參考點序列;

(5)重復步驟(2)~ (4),直到最大誤差或者參考點數目滿足要求,確定最終參考點序列;

(6)根據最終參考點序列進行RBF 插值,得到空間所有位置的網格尺度分布.

2.2 RBF 方法網格生成實例

本節仍然采用二維圓柱、NACA0012 翼型和30P30N 多段翼型作為考核算例,對RBF 方法進行實例測試.算例中所采用的三角形網格生成算法仍為基于ANN 的陣面推進法,具體可參考文獻[22].

RBF 基函數取為Wendland’s C0,緊支半徑取為計算域范圍大小的1/4.網格尺度控制對尺度插值誤差的要求不高,本文以尺度的最大相對誤差不超過5%為標準,進行參考點序列精簡,表1 給出了采用貪婪算法對參考點進行精簡的結果,結果顯示對于翼型算例,精簡序列后參考點數量至少減少一半,RBF 插值的耗時也減少了近一半,在保證了插值效果的同時提高了插值效率.

表1 RBF 方法參考點的數目及插值耗時Table 1 Number of reference nodes and time consumption on interpolation

圖5~ 圖7 給出了精簡后的參考點的位置及根據精簡參考點進行尺度控制生成的非結構網格.結果顯示RBF 方法能夠根據邊界網格尺度插值得到空間的網格尺度,插值得到的網格尺度分布過渡均勻.由于初始參考點的選擇是隨機的,參考點的分布具有一定的隨機性,但是均能夠保證插值誤差滿足要求.

圖5 圓柱算例精簡后的參考點及網格生成情況Fig.5 Reference nodes corresponding triangular mesh over a 2D cylinder

圖6 NACA0012 算例精簡后的參考點及網格生成情況Fig.6 Reference nodes and corresponding triangular mesh over NACA0012 airfoil

圖7 30P30N 算例精簡后的參考點及網格生成情況Fig.7 Reference nodes and corresponding triangular mesh over 30P30N airfoil

3 基于人工神經網絡的尺度控制方法

基于人工神經網絡的深度學習具有較強的非線性擬合能力,能夠通過現有樣本數據的訓練識別出數據中隱含的非線性映射關系.網格尺度分布實際上是一個關于幾何特征、流場特征的非線性映射,如圖8 所示,幾何特征(物面曲率、狹縫、細小結構等)可以直接影響網格尺度分布,幾何特征也可以決定流動特征(梯度量等),從而間接決定網格分布.不同幾何外形的網格尺度分布不同,同一幾何外形在不同來流條件(流場特征)情況下,網格尺度分布也不同.

圖8 網格分布控制與幾何特征和流場特征的關系Fig.8 Relationship between mesh size control,geometry,and flow features

流場特征與幾何外形、來流條件和邊界條件相關,目前采用機器學習方法進行流場預測是一個熱點研究問題[25-26].而根據流場特征確定網格分布則可以根據網格自適應的相關準則進行,也可以采用機器學習的方法進行[27].本文初步考慮幾何外形對網格分布的影響,采用ANN 建立幾何特征與網格尺度分布之間的關系.

3.1 ANN 輸入輸出參數

Chedid 和Najiar 等[28]曾嘗試用ANN 建立起網格密度與幾何特征之間的關系,如圖9 所示,神經網絡輸入考慮了空間點gi到邊界點的最小距離 δ1,次小距離 δ2及其對應的夾角a*和a**及到夾角對應兩邊的投影距離j1a,j1b,j2a和j2b共8 個輸入參數,輸出該空間點處的網格密度,定義為一定大小區域內節點的數量.該方法能夠較為全面地反映空間點與邊界之間的幾何關系(距離及投影距離),以及空間點對應的邊界處的幾何特征(夾角、曲率),但計算量較大不利于提高尺度控制效率,同時由于無法網格局部加密控制,失去了對網格尺度控制的靈活性,因此優勢并未得到發揮.

圖9 文獻[28]中神經網絡的輸入參數Fig.9 Input parameters for the artificial neural network in Ref.[28]

為提高效率,本文初步選擇最小壁面距離wdist作為輸入參數,網格尺度Sp作為輸出參數.為提高ANN 的泛化性,其輸入輸出參數通常需要進行歸一化操作,本文采用計算域的大小Lr_d,遠場邊界網格尺度Lr_f和物面邊界網格尺度Lr_w對輸入輸出進行分別歸一化.

同時,由于從物面到遠場,網格尺度變化較大,網格尺度相對值可能在103量級以上,為盡量縮小輸入輸出的值域范圍,提高ANN 訓練效果,本文對輸入輸出參數進行開根號.此外,在物面附近網格尺度變化快,需要同時考慮物面參考值和遠場參考值進行歸一化,具體輸入輸出參數形式如表2 所示.

表2 ANN 輸入輸出模型Table 2 Parameter model for artificial neural network

3.2 訓練方法及參數

本文基于Matlab 神經網絡訓練工具設計全連接的人工神經網絡,網絡含有1 個輸入層,1 個隱藏層,1 個輸出層.輸入層含有1 個神經元,輸出層含有1 個神經元,隱藏層的神經元數量為10 個,激活函數采用Sigmoid 函數,損失函數為均方誤差函數,訓練方法采用Levenberg-Marquardt 反向傳播方法,該方法比常規的梯度下降反向傳播算法訓練效率更高,具體可以參考文獻[29].神經網絡的結構如圖10所示.

圖10 基于Matlab 的人工神經網絡訓練工具Fig.10 Artificial neural network training tool based on Matlab

選擇二維圓柱網格和NACA0012 翼型網格作為訓練樣本,如圖11 所示,網格面個數分別為4295,3950,每個網格面對應一組樣本數據點.網格訓練樣本按70%,15%和15%的比例隨機劃分為訓練集、測試集和驗證集,圖12 給出了在翼型三角形網格訓練集上的Loss 值及在驗證集上的預測精度收斂歷程.結果顯示:經過120 次迭代后,Loss 值及預測誤差均下降到了0.001 58 左右.

圖11 網格分布訓練樣本網格Fig.11 Sample grids for ANN training

圖12 訓練Loss 值和精度收斂歷程Fig.12 Convergence of loss and accuracy on sample grids

3.3 ANN 預測結果

基于前述的ANN 訓練結果,分別預測生成了不同密度的二維圓柱網格、NACA0012 密網格、RAE2822 翼型以及三段翼型,在物面附近和遠場均取得了較好的效果,同時還能適應不同遠場大小的情況,如圖13 所示.算例中所采用的三角形網格生成算法仍為基于ANN 的陣面推進法,具體可參考文獻[22].

圖13 ANN 模型各向同性網格預測結果Fig.13 Mesh size controlled by ANN model for isotropic triangular grids

圖13 ANN 模型各向同性網格預測結果(續)Fig.13 Mesh size controlled by ANN model for isotropic triangular grids (continued)

同時,本文還將ANN 應用于各向異性混合網格尺度控制,各向異性四邊形采用層推進逐層推進生成[30-32],而各向同性三角形采用基于ANN 的陣面推進生成[22].層推進的推進方向、多方向推進數量均采用ANN 預測,同時在凹角處考慮局部推進步長,避免網格相交,具體方法細節可以參考文獻[30].為與各向同性網格尺度控制相一致,本文將層推進生成的最后一層網格作為虛擬物面,用于計算最小壁面距離,作為ANN 控制網格尺度的輸入參數.另外,層推進的網格尺度控制仍然采用指定物面第一層網格高度和增長率的方式給定.圖14 給出了NA0012翼型和30P30N 三段翼型的生成結果.由圖中結果可見,網格分布均勻合理,網格質量滿足要求,說明本文發展的ANN 方法可應用于各向異性混合網格的尺度控制.

圖14 ANN 模型各向異性混合網格預測結果Fig.14 Mesh size controlled by ANN model for anisotropic hybrid grids

4 背景網格法、RBF 方法與ANN 方法的比較

本文分別采用3 種方法生成了幾個典型幾何外形的非結構網格,其網格生成質量、生成效率、尺度控制靈活性、自動化程度等情況存在一定差別,本節對3 種方法進行對比分析.

從網格質量的角度,在網格尺度較小的部位,如翼型前后緣、狹縫等,背景網格法要求加密背景網格,以分辨最小網格尺度,而RBF 方法和ANN 方法不依賴于背景網格,因此在尺度較小的部位得到的尺度分布控制效果比采用矩形背景網格要好.

從自動化程度和靈活性角度,采用在背景網格上設置點源控制網格分布,對分布的控制最為靈活,能夠根據幾何特征和流場特征預先設置點源來改變網格分布.而背景網格法需要人工設置點源,本文算例中,在遠場和翼型附近人工設置數十個點源,每個點源人為給定強度、位置等參數,需要一定的經驗,自動化程度較低.而RBF 方法和ANN 方法只需要給定離散的邊界節點,就可以得到空間網格尺度分布,自動化程度相對較高,但是靈活性較低.

從網格生成效率的角度,背景網格法需要在笛卡爾網格上迭代求解泊松方程,在背景網格很密時,求解效率較低.而RBF 方法雖然需要求解矩陣的逆,但在引入貪婪算法之后,參考點數量減少,RBF 方法的效率得到提高.ANN 方法雖然需要求解最小壁面距離,但壁面距離的值并不需要十分精確,因此可以采用一些近似求法,也可以達到更高的效率.

表3 給出了3 種方法在生成網格過程中尺度控制所耗費的時間,由于每種方法生成的網格單元數存在一定差異,因此表中也給出了網格單元數量.

表3 3 種方法控制網格尺度耗時對比Table 3 Efficiency comparison of the three methods

背景網格法耗時基本由背景網格的規模及點源的數量決定.在本文三個算例中,背景網格規模為51 × 51,301 × 301 和401 × 401,點源的數量大致為20 個,因此在采用背景網格法時,翼型算例耗時明顯增加.

由表3 中數據可見,在背景網格數量較少時,背景網格法的效率較高,但是隨著背景網格數量和點源數量的增加,背景網格法效率明顯下降.相較于傳統背景網格法,RBF 方法和ANN 方法的耗時明顯減少,耗時僅為背景網格法耗時的1/10~ 1/5,網格尺度控制效率相應提高了5~ 10 倍.而且,在外形相對復雜的情況下,ANN 方法展現了更好的控制效率.可以預見,在三維復雜外形情況下,ANN 的控制效果和效率會更高.

5 總結與展望

本文提出了一種采用ANN 進行網格尺度分布控制的方法,初步確定了合理的神經網絡輸入輸出參數,基于Matlab 建立人工神經網絡模型,采用商業軟件生成二維圓柱和NACA0012 翼型非結構三角形網格作為網格密度訓練樣本,通過訓練和學習建立起壁面距離和網格尺度的映射關系,實現了不同密度的二維圓柱和不同二維翼型的網格尺度預測.同時,發展了基于RBF 方法的網格尺度控制方法,采用貪婪算法對插值參考點序列進行精簡,將RBF 插值效率提高了一倍.

與傳統背景網格法相比,基于ANN 方法和RBF方法的網格尺度預測效率提高5~10 倍,有助于進一步提高網格生成效率.最后將基于ANN 的網格尺度控制方法拓展應用于各向異性混合網格的尺度控制,得到的網格質量滿足要求,證明了方法的實用性.

展望ANN 方法在尺度分布控制領域的應用前景,在以往的CFD 研究和工程實踐中,已經積累了大量各種類型的網格,其中包含了網格尺度分布信息.對已有的網格數據按照幾何外形進行分類,訓練得到可以處理不同類別幾何外形的神經網絡,在需要對新的幾何外形進行網格尺度控制時,只需選擇已經訓練好的該類別的神經網絡即可高效快速完成尺度控制.進一步還可以考慮流場特征(來流條件),對神經網絡的適用范圍進行細分,使得其能夠同時考慮幾何外形和流場特征進行尺度控制.

可以進一步考慮物面幾何曲率、局部特征尺寸(窄邊、狹縫)等參數,提高神經網絡的泛化性,同時研究通過神經網絡進行流場預測,并根據流場特征進行網格自適應.

猜你喜歡
背景效率方法
“新四化”背景下汽車NVH的發展趨勢
《論持久戰》的寫作背景
當代陜西(2020年14期)2021-01-08 09:30:42
提升朗讀教學效率的幾點思考
甘肅教育(2020年14期)2020-09-11 07:57:42
晚清外語翻譯人才培養的背景
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
跟蹤導練(一)2
“錢”、“事”脫節效率低
中國衛生(2014年11期)2014-11-12 13:11:32
提高講解示范效率的幾點感受
體育師友(2011年2期)2011-03-20 15:29:29
主站蜘蛛池模板: 97视频免费看| 欧美h在线观看| 一区二区三区高清视频国产女人| 久久亚洲天堂| 国产精品成| 丰满人妻久久中文字幕| 国产成人综合亚洲欧美在| 亚洲精品图区| 午夜视频日本| 中文字幕欧美成人免费| 国产a v无码专区亚洲av| 国产激情无码一区二区APP | 精品久久久久成人码免费动漫| A级毛片高清免费视频就| 国产高潮流白浆视频| 尤物成AV人片在线观看| 国产福利微拍精品一区二区| 伊人成人在线| 91视频首页| 亚洲综合中文字幕国产精品欧美| 久久99精品久久久久纯品| 老色鬼欧美精品| 日韩AV手机在线观看蜜芽| 五月婷婷伊人网| 日韩免费毛片视频| 亚洲视频在线观看免费视频| 国内精自线i品一区202| a欧美在线| 精品视频第一页| 中文成人在线视频| 丰满的熟女一区二区三区l| 亚洲床戏一区| 亚洲V日韩V无码一区二区| 国产亚洲视频免费播放| 亚洲欧美日韩色图| 成人精品亚洲| 四虎永久在线| 亚洲免费福利视频| 青草视频在线观看国产| 亚洲网综合| 亚洲欧美日韩成人在线| 成人在线不卡| 国产网友愉拍精品| 天天色综网| 亚欧美国产综合| 国产精品亚洲五月天高清| 野花国产精品入口| 久久人妻系列无码一区| 国产成人精品2021欧美日韩| 亚洲欧洲日产无码AV| 亚洲日韩每日更新| 国产精品美女网站| 五月天综合婷婷| 一级毛片无毒不卡直接观看 | 高清无码一本到东京热| 亚洲男人的天堂久久香蕉| 91精品啪在线观看国产91| 第一区免费在线观看| 国产精品久久久久久久久| 国产精品香蕉| 亚洲一区无码在线| 国产精品久久久久久影院| 亚洲经典在线中文字幕| 亚洲五月激情网| 美女潮喷出白浆在线观看视频| 午夜久久影院| 中文字幕永久在线观看| 国产乱视频网站| 亚洲欧美自拍中文| 久久无码av三级| 国产精品国产三级国产专业不| 人妻精品久久无码区| 日韩区欧美国产区在线观看| 久草视频中文| 国产99精品久久| 国产成人精品亚洲日本对白优播| 成人福利在线观看| 日本免费福利视频| 国产爽爽视频| 日韩黄色精品| 一本一道波多野结衣av黑人在线| 精品三级网站|