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

基于譜元法和北方蒼鷹算法的振動損傷檢測

2024-01-01 00:00:00周運來姚峰白春玉李凱翔朱勝陽
華東交通大學學報 2024年5期

摘要:【目的】為了解決結構損傷檢測中有限元法建模不精確和計算成本高的問題,提出了一種結合譜元法與北方蒼鷹優化算法(NGO)的結構健康監測(SHM)技術。【方法】首先,采用譜元法建立結構的頻響函數,并應用于結構的損傷定位與損傷檢測目標函數的構造,將損傷檢測劃分為兩階段問題,以降低算法優化維度和檢測復雜性。其次,引入北方蒼鷹優化算法(NGO),對目標函數進行優化求解。最后,以平面桁架結構和ASCE Benchmark結構為研究對象,利用NGO、粒子群優化(PSO)和蟻獅優化(ALO)算法對其各種損傷工況進行損傷檢測性能對比。【結果】結果表明,在低維度和簡單結構中,NGO,PSO和ALO算法均表現出良好的求解能力;但在高維度和大型復雜結構中,NGO相較于PSO和ALO算法具有更高的損傷檢測能力和魯棒性。【結論】改進后的方法提高了損傷檢測數值建模的精度。

關鍵詞:譜元法;損傷檢測;NGO;ASCE Benchmark Structure

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

本文引用格式:周運來,姚峰,白春玉,等. 基于譜元法和北方蒼鷹算法的振動損傷檢測[J]. 華東交通大學學報,2024,41(5):29-38.

Vibration Damage Detection Based on the Spectral Element Method and Northern Goshawk Optimizer Algorithm

Zhou Yunlai1, Yao Feng1, Bai Chunyu2, Li Kaixiang2, Zhu Shengyang3

(1. State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace Engineering, Xi′an Jiaotong University, Xi′an 710049, China; 2. National Key Laboratory of Strength and Structural Integrity, Aircraft Strength Research Institute of China, Xi′an 710065, China; 3. State Key Laboratory of Rail Transit Vehicle System, Southwest Jiaotong

University, Chengdu 611756, China)

Abstract: 【Objective】 To address the issues of imprecise modeling and high computational cost in structural damage detection using the finite element method (FEM), this study proposes a structural health monitoring (SHM) technique that combines the spectral element method (SEM) with the Northern Goshawk Optimization (NGO) algorithm.【Method】Firstly, the spectral element method was used to establish the frequency response function of the structure, which was then applied to construct the objective function for damage localization and detection. This approach divided the damage detection problem into two stages, reducing the optimization dimension and complexity. Secondly, NGO algorithm was introduced to optimize and solve the objective function. Finally, planar truss structure and ASCE Benchmark Structure were used as case studies to compare the damage detection performance of NGO, Particle Swarm Optimization (PSO), and Ant Lion Optimization (ALO) algorithms under various damage cases. 【Result】The results show that for low-dimensional and simple structures, NGO, PSO, and ALO algorithms all exhibit good solving capabilities. However, for high-dimensional and large complex structures, NGO demonstrates superior damage detection capability and robustness compared to PSO and ALO. 【Conclusion】 The improved method enhances the accuracy of numerical modeling in damage detection.

Key words: spectral element method; damage detection; NGO; ASCE Benchmark Structure

Citation format: ZHOU Y L, YAO F, BAI C Y, et al. Vibration damage detection based on the spectral element method and Northern Goshawk Optimizer algorithm[J]. Journal of East China Jiaotong University, 2024, 41(5): 29-38.

【研究意義】在結構的使用壽命期間,其性能會因各種損傷逐步衰退,主要源于材料退化、結構疲勞、超載及環境因素。這些損傷可能改變結構材料的性質或幾何特性。為確保結構的穩定性、維護性和可靠性,早期準確定位和診斷損傷尤為重要。簡而言之,準確識別和定位結構損傷對于保持結構完整性和安全性至關重要,有助于防止漸進式故障并提高安全性能[1-2]。

【研究進展】過去數十年中,基于多種元啟發式算法對振動分析數據進行處理以評估結構損傷的方法,已經受到了學術界和工業界的廣泛關注[3-4]。Majumdar等[5]采用有限元法計算的固有頻率與實際測量結果之間的均方根值作為目標函數,使用蟻群優化算法進行結構損傷檢測,并驗證了平面和空間桁架結構。Mohan等[6]對懸臂梁和平面框架結構進行損傷檢測,采用基于結構頻域響應的差異和基于固有頻率的差異兩種目標函數,并用粒子群優化算法最小化目標函數,結果顯示基于頻域響應的方法精度更高。Nabavi等[7]提出結合巨型犰狳優化算法(GOA)與結構時域響應的損傷檢測方法,通過結構損傷加速度和解析模型的加速度定義優化目標函數,用GOA算法解決損傷識別問題,確定結構損傷位置和程度。李成等[8]分別針對兩層鋼架以及三層框架結構,以實測模態參數和計算模態參數之間的誤差作為目標函數,采用人工魚群算法進行優化求解,從而實現結構的損傷檢測。Seyedpoor等[9]采用兩階段方法識別桁架結構損傷:第1階段使用柔度矩陣損傷指數識別疑似損壞單元,第2階段用差分進化算法預測損傷程度。Khatir等[10]提出改進頻響函數(FRF)指標用于復雜結構損傷識別,第一階段檢測和定位單個及多個損傷,之后用人工大猩猩部隊優化算法(GTO)、野狗優化算法(DOA)、非洲禿鷲優化算法(AVOA)和基于梯度的優化算法(GBO)等求解損傷程度。

【關鍵問題】結構建模在損傷檢測中至關重要。之前的研究中,均采用有限元法進行損傷檢測數值建模。而有限元法需要細化網格以獲得準確結果,這增加了結構自由度和計算成本[11]。譜元法結合了有限元法、動力學剛度法和譜分析法,具有高精度、少自由度、低計算成本等優點,能夠有效求解頻域問題[12-13]。北方蒼鷹優化(northern goshawk optimization, NGO)算法是一種啟發式優化算法,于2022年由Dehghani等提出[14]。該算法通過模擬蒼鷹的狩獵行為,包括獵物識別與攻擊、追逐逃生兩個階段,來尋找優化問題的解,具有結構簡單、性能優良等特點,已在多個領域得到了應用[15-16]。

【創新特色】本文結合譜元法與NGO進行損傷檢測。首先,利用譜元法構造結構的頻響函數,彌補了有限元法在高頻段精度不高的缺點,并利用頻響函數建立損傷定位模型以及損傷檢測目標函數。其次,引入NGO,對目標函數進行優化,并與PSO以及ALO算法進行比較。最后,分別以平面桁架結構以及ASCE Benchmark結構為研究對象,對其各種損傷工況進行檢測。為了驗證NGO算法的魯棒性并模擬真實測量數據,在數值計算的頻響曲線中加入5%的高斯白噪聲。

1 損傷檢測原理

1.1 譜元法求解頻響函數

1.1.1 桿的軸向振動

桿單元的運動方程為

[ρA?2ux,t?t2-EA?2ux,t?x2=0] (1)

式中:[ux,t]為桿的縱向位移;ρ,A和E分別為桿的密度,截面面積和彈性模量。

桿單元的動力學剛度矩陣為

[SRω=EALkLLsinkLLcoskLL-1-1coskLL] (2)

式中:L為長度;

[kL=ωρE] (3)

則桿單元頻域上節點力與節點位移的關系為

[SRωU1U2=F1F2] (4)

式中:U1,U2均為節點位移;F1,F2均為節點軸向力。

1.1.2 Euler梁的橫向振動

Euler梁的振動方程為

[EI?4w?x4+ρA?2w?t2=0] (5)

Euler梁單元的動力學剛度矩陣[SBω]為

[SBω=EIL3αλL-ξγLλLβL2-γLμL2-ξ-γLα-λLγLμL2-λLβL2] (6)

式中:L為長度;I為截面慣性矩,

[α=CSh+SChkL3/Δ ,ξ=S+ShkL3/Δβ=-CSh+SChkL/Δ ,μ=-S+ShkL/Δγ=-C+ChkL2/Δ ,λ=SShkL2/ΔC=coskL , " S=sinkL ,Ch=cosh kL ,Sh=sinh kLΔ=1-CCh , " " "k=ωρAEI1/4] 則Euler梁頻域上節點力與節點位移的關系為

[SBων1θ1ν2θ2=V1M1V2M2] (7)

式中:[v1],[θ1]分別為左節點上的位移和轉角;[v2],[θ2]分別為右節點上的位移和轉角;[V1],[M1]分別為左節點上的剪力和彎矩;[V2],[M2]分別為右節點上的剪力和彎矩。

1.1.3 軸的扭轉振動

軸的扭轉振動方程為

[?2φx,t?x2-ρG?2φx,t?t2=0] (8)

式中:[φx,t]為桿的縱向位移;ρ和G分別為桿的密度和剪切模量。

軸扭轉的動力學剛度矩陣為

[STω=kTGJsinkTLcoskTL-1-1coskTL] (9)

式中:L為長度;J為截面極慣性矩;

[kT=ωρG] (10)

則頻域上節點扭矩與節點扭角的關系為

[STωΦ1Φ2=T1T2] (11)

式中:Φ1,Φ2均為節點扭角;T1,T2均為節點扭矩。

1.1.4 頻響函數

以空間梁單元為例,將式(4),式(7),式(11)整合,得到空間梁單元的節點力與節點位移之間的關系。將其利用有限元剛度陣組合方法,得到空間結構的譜形式方程為

[SωUω=Fω] (12)

[Uω=Sω-1Fω] (13)

式中:[Sω]為空間結構動力學剛度矩陣;[Uω]為節點的頻域響應向量;[Fω]為頻域的節點力向量。

由式(13)可得出,空間結構的頻響函數為

[Hω=Sω-1] (14)

1.2 損傷定義

賦予結構中每個單元一個損傷系數a,則損傷的動力學剛度矩陣可以表示為

[SωDe=1-aSω] (15)

式中:[Sω]為未損傷的動力學剛度矩陣;[SωDe]為損傷的動力學剛度矩陣;a為損傷系數,在0~1取值,a為0時表示該單元沒有損傷,a為1時表示該單元完全損傷。

1.3 損傷定位原理

通過譜元模型模擬了未損傷結構的頻響函數為

[HωA=Sω-1A] (16)

損傷結構的頻響函數為

[HωT=Sω-1T] (17)

因損傷引起的剛度變化ΔS[ω]為

[ΔSω=Hω-1A-Hω-1T] (18)

式(18)兩邊乘以[Hω]T得到

[HωT[ΔSω]=HωTHω-1A-HωTHω-1T] (19)

β值是根據自由度、元素位置以及全局FRF矩陣的第j行計算得到的,用于表示損傷位置,其定義為

[βj,i=H(ω)T1nH(ω)A-1(:,i)-I(j,i)] (20)

1.4 損傷檢測目標函數

將結構的損傷檢測問題轉化為一個優化問題,使損傷結構計算的頻響函數與實際的頻響函數之間的誤差最小,從而得到損傷參數。基于頻響函數的誤差函數如下

[Eqa=j=1mi=1nHcomputeωi,a-Hactualωi] (21)

式中:m為考慮響應的節點總數;n為計算時采用的頻率總數;[Hcomputeωi,a]為計算的頻響函數;[Hactualωi]為實際的頻響函數。

2 北方蒼鷹算法(NGO)

北方蒼鷹的狩獵策略可以分為兩個階段:獵物識別階段,追擊與逃逸階段。

2.1 第1階段(獵物識別階段)

在捕獵的第1階段,北方蒼鷹隨機挑選一個獵物,隨后快速向其發起攻擊。由于在搜索空間中對獵物的選擇是隨機的,因此這一階段增加了NGO算法的勘探能力。這個階段對搜索空間進行全局搜索,目的是確定最優區域。在這一階段,北方蒼鷹進行獵物選擇和攻擊的行為,用式(22)~式(24)描述為

[Pi=Xk,i=1,2,…,N,k=1,2,…,i-1,i+1,…,N] (22)

[xnew,P1i,j=xi,j+rpi,j-sxi,j,FPilt;Fixi,j+rxi,j-pi,j,FPi≥Fi] (23)

[Xi=Xnew,P1i,Fnew,P1ilt;FiXi,Fnew,P1i≥Fi] (24)

式中:[Pi]為第i只北方蒼鷹所選獵物的位置;[FPi]為對應的目標函數值;k為[1,N]范圍內且不等于i的隨機整數;[xnew,P1i,j]為第i只北方蒼鷹的新解在第j維的值;[Fnew,P1i]為其對應的目標函數值;r為[0,1]范圍內的隨機數;s的值為1或2。

2.2 第2階段(追擊及逃逸)

在北方蒼鷹攻擊獵物后,獵物會試圖逃跑。因此,在追逐獵物的收尾過程中,北方蒼鷹需要繼續追逐獵物。由于北方蒼鷹的追擊速度很高,幾乎可以在任何情況下追逐獵物,并最終捕獲獵物。對這種行為的模擬提高了算法對搜索空間的局部搜索能力。假設這種狩獵活動接近于一個半徑為R的攻擊位置。在第2階段中,用式(25)~式(27)描述為

[xnew,P2i,j=xi,j+R2r-1xi,j] (25)

[R=0.021-tT] (26)

[Xi=Xnew,P2i,Fnew,P2ilt;FiXi,Fnew,P2i≥Fi] (27)

式中:t為當前迭代次數;T為最大迭代次數;R為攻擊半徑,隨迭代次數增加而縮小;[xnew,P2i,j]是第i只北方蒼鷹在新位置的第j維解的值;[Fnew,P2i]為其對應的目標函數值。

3 數值算例

本節將討論基于譜元法和NGO算法的損傷檢測性能評估。本研究以平面桁架結構和ASCE Benchmark結構為研究對象,并與ALO以及PSO算法進行損傷檢測性能比較,為防止偶然性,每種算法進行10次求解,并取平均值。另外,在實際測量結構頻響的實驗中,往往會產生誤差,從而導致頻響曲線產生噪聲。為了模擬實際測量結果,在數值計算的頻響曲線中加入5%的高斯白噪聲。

3.1 平面桁架結構損傷檢測

圖1為平面桁架結構,將平面桁架損傷分為3種類型,如表1所示。

在節點6處施加豎直方向單位力,利用譜元法分別計算節點2,節點3,節點5垂直自由度的頻響,將其作為實際測量數據。圖2為平面桁架在節點3垂直自由度處未損傷與3種損傷工況的頻響曲線。

3.1.1 平面桁架結構損傷工況1

假設頻響是在450 rad/s的測試頻率下測量的。利用全局頻響矩陣的第4行,式(20)中向量[β4,i]為

[β4,i=[0, 0, 0, 0, 0, 0, -0.549 "1, 0.378 "1, -0.235 "1, " " " " " " " " " " " " " " " 0, 0, 0, 2.205 "7, -1.728 "4, -0.165 "4, 0, 0, 0]]

在[β4,i]中,除了第7~9,13~15自由度外,其余全為0。則可以判斷出節點3和節點5所對應的單元有損傷,即疑似損傷單元為單元6。

在完成損傷定位之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取30,求解結果如圖3所示。結果表明,3種算法在損傷工況1中均表現良好。

3.1.2 平面桁架結構損傷工況2

損傷工況2考慮兩個單元損傷,利用全局頻響矩陣的第4行,此時式(20)中向量[β4,i]為

[β4,i=[0, 0, 0, 2.712 "8, 1.486 "3, 0.986 "4, -2.052 "0, " " " " " " " "- 0.874 "3, 0.714 "9, " " " " " " " 0, 0, 0, 0, 0, 0, 0.012 "2, " " " " " " " "- 1.059 "1,-0.051 "5]]

在[β4,i]中,除了第4~9,16~18自由度外,其余全為0。則可以判斷出節點2,節點3以及節點6所對應的單元有損傷,即疑似損傷單元為單元2和單元7。

損傷定位完成之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取50,求解結果如圖4所示。結果表明,3種算法在損傷工況2中均有較強的求解能力。

3.1.3 平面桁架結構損傷工況3

損傷工況3考慮兩個單元損傷,利用全局頻響矩陣的第4行,此時式(20)中向量[β4,i]為

[β4,i=[0, 0, 0, -1.976 "9, 0.619 "5, 0.023 "8, 1.210 "0, " " " " " " " " 0.446 "5,-0.233 "5,0, 0, 0, -1.402 "1, " " " " " " " "- 0.523 "3, -0.354 "5, 0, 0, 0]]

在[β4,i]中,除了第4~9,13~15自由度外,其余全為0。則可以判斷出節點2,節點3以及節點5所對應的單元有損傷,即疑似損傷單元為單元2,單元5以及單元6。

損傷定位完成之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取50,求解結果如圖5所示。結果表明,NGO和PSO在損傷工況3中表現良好,ALO算法求解能力略差。

3.2 ASCE Benchmark Structure損傷檢測

1998年,Black和Ventura提出了Benchmark結構,并建造于加拿大哥倫比亞大學地震工程研究實驗室內。圖6為Benchmark Structure的三維示意圖,表2為其結構參數。由45個節點和116個單元組成。為了研究結構的損傷識別,UBC實驗室針對Benchmark Structure定義了不同損傷工況,主要通過改變剛度對損傷進行模擬。參考其設置方法,本文設置了4種損傷工況,分別模擬了框架中各構件的不同程度損傷,如表3所示。

在節點17的y方向施加單位力,利用譜元法分別計算節點21,節點35,節點33y方向的頻響。將其作為實際測量數據。圖7為Benchmark結構在節點33y方向未損傷與4種損傷工況的頻響曲線。

3.2.1 Benchmark結構損傷工況1

損傷工況1僅考慮一個單元損傷。利用基于頻響函數的損傷定位方法,可以定位出,節點19和節點28所對應的單元疑似損傷,即單元103疑似損傷。

損傷定位完成之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取50,求解結果如圖8所示。結果表明,NGO和ALO算法在損傷工況1中求解結果優于PSO算法。

3.2.2 Benchmark結構損傷工況2

損傷工況2為第2層斜支撐全部移除。利用基于頻響函數的損傷定位方法,可以定位出單元93~單元100疑似損傷。

單元損傷定位完成之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取200,求解結果如圖9所示。結果表明,NGO算法和ALO算法在損傷工況2中,求解能力遠強于PSO。

3.2.3 Benchmark結構損傷工況3

損傷工況3考慮兩個單元受損。利用基于頻響函數的損傷定位方法,可以定位出單元9以及單元58疑似損傷。

損傷定位完成之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取200,求解結果如圖10所示。結果表明,NGO和PSO算法在損傷工況4中,求解能力優于ALO算法。但從求解迭代曲線中可以看出,NGO算法收斂速度比PSO算法快。

3.2.4 Benchmark結構損傷工況4

損傷工況4考慮3個單元受損。利用基于頻響函數的損傷定位方法,可以定位出單元3,單元12以及單元24疑似損傷。

損傷定位完成之后,分別利用NGO,PSO,ALO算法對式(21)進行優化求解。最大迭代次數取300,求解結果如圖11所示。結果表明,NGO算法在損傷工況4中,求解能力優于PSO以及ALO算法。

4 結論

提出了一種基于譜元法和NGO的方法,利用譜元法求解結構的頻響曲線,并以頻響函數為指標,采用了兩階段損傷檢測方法,減少了求解的維度,更適合用于較復雜結構的損傷檢測。以平面桁架和ASCE Benchmark Structure為例,在MATLAB中,將NGO算法與ALO以及PSO算法進行損傷檢測性能對比。可以得到如下結論。

1) 譜元法可以作為一種更高效準確的損傷檢測建模方法來代替有限元法。

2) 在低維度以及較簡單結構中,NGO,PSO以及ALO算法具有良好的求解能力。

3) 對于較高維度以及復雜結構來說,NGO相較于PSO與ALO具有更高的損傷檢測能力以及魯棒性。

參考文獻:

[1] " 巫文君, 董利強, 任暉, 等. 基于改進MPSCO算法的框架結構損傷檢測研究[J]. 地震工程與工程振動, 2017, 37(2): 108-116.

WU W J, DONG L Q, REN H, et al. Damage detection of frame structures based on improved multi-particle swarm coevolution optimization algorithm[J]. Earthquake Engineering and Engineering Dynamics, 2017, 37(2): 108-116.

[2] " MAHDAVI S H, XU C. Time-domain structural damage identification using ensemble bagged trees and evolutionary optimization algorithms[J]. Structural Control and Health Monitoring, 2023: 6321012.

[3] " MORADI S, RAZI P, FATAHI L. On the application of bees algorithm to the problem of crack detection of beam-type structures[J]. Computers amp; Structures, 2011, 89(23/24): 2169-2175.

[4] " 陳震, 朱軍華, 余嶺. 一種基于改進PSO算法的結構損傷識別方法[J]. 振動與沖擊, 2012, 31(5): 17-20.

CHEN Z, ZHU J H, YU L. An improved PSO algorithm for structure damage identification[J]. Journal of Vibration and Shock, 2012, 31(5): 17-20.

[5] " MAJUMDAR A, MAITI D K, MAITY D. Damage assessment of truss structures from changes in natural frequencies using ant colony optimization[J]. Applied Mathematics and Computation, 2012, 218(19): 9759-9772.

[6] " MOHAN S C, MAITI D K, MAITY D. Structural damage assessment using FRF employing particle swarm optimization[J]. Applied Mathematics and Computation, 2013, 219(20): 10387-10400.

[7] " NABAVI S, GHOLAMPOUR S, HAJI M S. Damage detection in frame elements using Grasshopper Optimization Algorithm (GOA) and time-domain responses of the structure[J]. Evolving Systems, 2022, 13(2): 307-318.

[8] " 李成, 余嶺. 基于人工魚群算法的結構模型修正與損傷檢測[J]. 振動與沖擊, 2014, 33(2): 112-116.

LI C, YU L. Structural model updating and damage detection based on artificial fish swarm algorithm[J]. Journal of Vibration and Shock, 2014, 33(2): 112-116.

[9] " SEYEDPOOR S M, MONTAZER M A. Damage identification method for truss structures using a flexibility-based damage probability index and differential evolution algorithm[J]. Inverse Problems in Science and Engineering, 2016, 24(8): 1303-1322.

[10] KHATIR S, TIACHACHT S, THANH C L, et al. A robust FRF damage indicator combined with optimization techniques for damage assessment in complex truss structures[J]. Case Studies in Construction Materials, 2022, 17: e01197.

[11] LEE U. Spectral element method in structural dynamics[M]. Hoboken: John Wiley amp; Sons, 2009.

[12] 鄂林仲陽, 杜強, 李上明. 基于譜元法的空間剛架動力學特性分析[J]. 計算力學學報, 2016, 33(5): 802-806.

E L Z Y, DU Q, LI S M. Dynamics characteristic analysis of space frame structure based on spectral element method[J]. Chinese Journal of Computational Mechanics, 2016, 33(5): 802-806.

[13] KIM T, LEE U. Vibration analysis of thin plate structures subjected to a moving force using frequency-domain spectral element method[J]. Shock and Vibration, 2018: 1908508.

[14] DEHGHANI M, HUBáLOVSKY ?, TROJOVSKY P. Northern Goshawk Optimization: a new swarm-based algorithm for solving optimization problems[J]. IEEE Access, 2021, 9: 162059-162080.

[15] EL-DABAH M A, EL-SEHIEMY R A, HASANIEN H M, et al. Photovoltaic model parameters identification using Northern Goshawk Optimization Algorithm[J]. Energy, 2023, 262: 125522.

[16] NING M, CHEN X, LIN Y, et al. Revealing the hot deformation behavior of AZ42 Mg alloy by using 3D hot processing map based on a novel NGO-ANN model[J]. Journal of Materials Research and Technology, 2023, 27: 2292-2310.

通信作者:周運來(1986—),男,博士,研究員,博士生導師。研究方向為結構優化、無損監測。E-mail: yunlai.zhou@xjtu.edu.cn。

主站蜘蛛池模板: 亚洲AⅤ波多系列中文字幕| 亚洲综合18p| 国产一区成人| 国产成人精品午夜视频'| 国产久草视频| 国产女主播一区| 欧美精品二区| 97视频免费在线观看| 欧美黑人欧美精品刺激| 97视频免费在线观看| 中文字幕精品一区二区三区视频| 亚洲视频四区| 久无码久无码av无码| 5555国产在线观看| 免费人成视网站在线不卡| 国产黑人在线| 欧美一级高清视频在线播放| 免费无码AV片在线观看中文| 3344在线观看无码| 国产精品流白浆在线观看| 99re在线免费视频| 九九线精品视频在线观看| 国产成人综合久久精品下载| 久久久久国产精品熟女影院| 国产无码在线调教| 久久性妇女精品免费| 国产99视频在线| 五月天综合网亚洲综合天堂网| 国产成人精品2021欧美日韩| 久久国产精品麻豆系列| 国产国拍精品视频免费看| 97久久免费视频| 在线观看国产精品日本不卡网| 精品国产Av电影无码久久久| 国产资源免费观看| 国产亚洲精品91| 日韩av手机在线| 国产一级二级三级毛片| 亚洲精品视频在线观看视频| 国产欧美性爱网| 亚洲一本大道在线| 精品亚洲麻豆1区2区3区| 国产精品99久久久| 国产福利在线免费| julia中文字幕久久亚洲| 欧美性色综合网| 亚洲免费成人网| 国产一区免费在线观看| 国产精选自拍| 日本少妇又色又爽又高潮| 日韩专区第一页| 日韩免费毛片| 怡红院美国分院一区二区| 尤物视频一区| 欧美日韩资源| 无码视频国产精品一区二区| 国国产a国产片免费麻豆| 亚洲第一天堂无码专区| 成人小视频网| 午夜久久影院| 国产成人1024精品下载| 亚洲色图综合在线| 巨熟乳波霸若妻中文观看免费| 色欲色欲久久综合网| 国产区福利小视频在线观看尤物| 2021国产在线视频| 福利一区在线| 激情影院内射美女| 白浆视频在线观看| 一区二区三区国产精品视频| 久久婷婷五月综合色一区二区| 青青操国产| 国产哺乳奶水91在线播放| 无码AV高清毛片中国一级毛片| 色综合色国产热无码一| 67194亚洲无码| 亚洲精品视频在线观看视频| 国产精鲁鲁网在线视频| 97国产在线观看| 亚洲综合片| 久久夜夜视频| 制服丝袜国产精品|