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

考慮RNA運行作用的近海風電結構主動質量阻尼器振動控制研究

2023-12-18 08:59:32白久林李晨輝王宇航
振動與沖擊 2023年23期
關鍵詞:風速結構模型

白久林, 李晨輝, 王宇航

(1.重慶大學 山地城鎮建設與新技術教育部重點實驗室,重慶 400045;2.重慶大學 土木工程學院,重慶 400045; 3.西北電力設計院有限公司,西安 710075)

發展風電產業是我國“十四五”期間實現能源低碳轉型并推動形成綠色發展方式的重要途徑。近年來,為了達成更高的新能源供能比例,我國風電裝機規模不斷擴大,塔架不斷升高,風輪尺寸不斷增加。然而,如此迅猛的發展中存在兩個重要問題:① 風電市場擴張迅速,對風電項目的開發“降本”提出了更高的要求,開發新型的高性能風電結構體系是實現“降本”目標的有效選擇。② 海上風電塔支撐結構常采用單樁基礎和鋼塔筒,是典型的高柔結構,其基準周期長,在環境荷載與風輪-機艙組件(rotor-nacelle assembly, RNA)的運行作用激勵下,振動問題突出,倒塔事故頻發[1-2],風電結構體系“增效”需求迫切。

質量阻尼器是結構振動控制的常用手段,主要包括被動、半主動、主動三種形式,其通過附加質量產生的慣性作用為主體結構提供相對反向的控制力,同時振動能量由阻尼元件逐漸耗散。采用質量阻尼器的有益效果體現為:① 可降低風電結構內力,使塔架更容易滿足設計要求,節省塔筒材料、運輸、施工等建設成本;② 可緩解結構的疲勞損傷,延長使用年限。因此,質量阻尼器可解決“降本”與“增效”兩方面問題,為電能輸出提供保障,具有實際意義[3-4]。

目前,學者們采用不同的質量阻尼器方案對風電塔進行了控制研究。在被動阻尼器方面,Murtagh等[5]考慮塔架和旋轉葉片的自振特性,建立了風電塔耦合模型,初步探究了調諧質量阻尼器(tuned mass damper, TMD)對風電塔的風致振動控制。Lackner等[6]采用自主開發的FAST-SC平臺,研究了TMD對不同基礎形式風電結構的控制作用。Li等[7]研究了擺式TMD對臺風作用下風電結構動力響應的控制作用。Stewart等[8-9]使用FAST-SC平臺研究了風浪方向錯位下TMD對海上風電結構的控制效果。Sun等[10]考慮風浪錯位和渦流致振,提出一種三維擺式調諧質量阻尼器(3D-PTMD)來抑制近海風電結構多方向的動力響應,并先后研究了3D-PTMD在疲勞損傷控制[11]及冰激振動控制[12]方面的表現。李萬潤等[13]提出一種雙向TMD,開展了風電結構風致響應減振的研究。在半主動阻尼器的研究中,線性二次型調節(linear quadratic regulator, LQR)算法是一種有效的控制策略:Sarkar等[14-15]與Rezaee等[16]分別基于LQR算法考察了半主動調諧液柱阻尼器(SA-TLCD)與半主動調諧質量阻尼器(SA-TMD)對風電塔動力響應的控制作用。

白久林等[17-18]的研究結果表明,RNA運行作用下風電塔頂前后向響應含量復雜;在氣動阻尼的影響下塔架一階振動受到限制。然而,在已有研究中,風電結構普遍被考慮為傳統的高聳結構針對塔架一階頻率進行質量阻尼器的調諧與設計。Ghassempour等[19]、Chen等[20-21]與Xie等[22]分別通過不同的耦合模型(GH Bladed,FAST,Simpack)得出了類似的結論:風機運行狀態下,針對塔架一階頻率調諧的TMD對前后向響應的控制效果會受到負面影響??紤]到風電塔前后向為主要受載向,塔頂加速度是工程中的重要監測指標,加速度超限可能會導致風機故障甚至局部脫網,因此發展一種控制效率更高的阻尼器具有重要的應用價值。

與TMD不同,主動質量阻尼器(active mass damper, AMD)[23]結構體系,如圖1所示。通過主動的驅動元件直接與受控結構相連,通過一定反饋算法實時輸出驅動力,具有無需調諧的特點,可更好地適應風電結構復雜的動力特性,其對風電結構的控制研究目前幾乎空白。本文基于LQR算法,探究了考慮RNA運行作用下,AMD對風電結構前后向響應的控制作用,從時域和頻域分別考察了其減振性能與控制特點,旨在為風電結構減振提供參考,為風電塔的結構安全與供能穩定提供支撐與保障。

1 近海風電結構分析模型

1.1 NREL 5 MW 單樁風電塔

美國可再生能源實驗室(National Renewable Energy Laboratory, NREL) 5 MW風機模型[24]是目前最常用的風機模型,其開發綜合了Multibrid M5000、REpower 5M等實際風機數據,可真實地反應風機各項動力特性,部分數據如表1所示。本文以最常見的單樁風電塔為對象開展系列分析,如圖2所示。

表1 NREL 5 MW單樁風電塔參數

圖2 NREL 5 MW 單樁風電塔

風電結構的顯著特點在于塔筒頂部用于捕獲風能的RNA。在風力發電的過程中,為了風能的最大利用,風機常采用變速-變槳控制策略進行發電量的調節。根據貝茲理論[25],風輪捕獲功率P的表達式為

(1)

式中:ρa為空氣密度,取1.225 kg/m3;v為風速;R為風輪半徑;功率系數Cp為尖速比λ與槳距角β的函數。其中尖速比λ為葉片尖端線速度(ωr為角速度)與風速v的比值

(2)

可知,固定風速時,功率系數Cp決定了風輪捕獲功率。實際運行中,Cp取決于風輪轉速與槳距角的實時調節:當實時風速處于切入風速與額定風速之間時,葉片槳距角維持0°,通過變速控制使風機始終運行在最佳尖速比λopt附近,實現功率系數Cp最大化;當風速超過額定風速時,通過槳距角調節(變槳)控制功率系數Cp,從而使發電機以額定功率穩定輸出。

對于NREL 5 MW風機,其運行過程分為三種狀態:風速在3.0~11.4 m/s為切入狀態,在11.4~25.0 m/s為額定狀態,超過25.0 m/s為停機狀態,如圖3所示。

圖3 變速-變槳控制策略

1.2 單自由度模型

由于細長、柔性的風電塔在風荷載作用下彎曲變形形式單一,振動位移沿塔高的分布基本符合一階振型,可在考慮RNA運行作用對體系參數影響的基礎上,將風電塔簡化為單自由度(single degree of freedom, SDOF)運動模型[26-28],如圖4所示。

圖4 風電結構單自由度模型示意圖

風電塔SDOF體系的各項參數需考慮RNA運行作用的影響,結合圖3,選取兩特定風速條件7 m/s (切入狀態)與18 m/s (額定狀態),使用OpenFAST軟件[29]分別計算不同狀態下耦合模型前后向自由衰減曲線,如圖5所示。質量、頻率及阻尼比可具體考慮為:

圖5 OpenFAST衰減測試

(1) 質量:本文聚焦塔頂前后向加速度的控制方案,因此將塔葉片和機艙簡化為集中質量,總計350 000 kg;塔筒作為支撐結構,主要提供體系整體剛度與自身阻尼,其質量在重點研究塔頂響應的等效模型中可以忽略。

(2) 頻率:由于塔架在環境荷載下運動行為主要受結構一階振型控制,因此等效模型僅考慮前后向一階模態。對應主頻率根據圖5的衰減曲線相鄰振峰間距可計算前后向一階頻率fs=0.278 Hz,與表1一致。

(3) 阻尼:風電塔前后向阻尼不可采用塔架結構阻尼 (約0.005~0.010) ,而應考慮包括風輪旋轉產生氣動阻尼的整體阻尼,通過動力衰減法計算切入狀態整體阻尼為0.059,額定狀態整體阻尼為0.076,如圖5所示。

1.3 氣動荷載

SDOF體系需要輸入RNA運行作用下的氣動荷載。對于塔架前后向,主要考慮風輪旋轉時的氣動推力Taero,由風場條件、風機尺寸、風輪轉速、槳距角等參量共同決定。

風況的考慮包含平均風與湍流風。平均風速可采用指數律風剖面模型以考慮風輪表面的風剪切效應,表示為

(3)

式中:Uz為不同高度z處平均風速;Uhub為輪轂高度zhub=90 m處的平均風速;α為粗糙度指數,IEC-61400-1[30]建議取0.14。

湍流風采用Von Karman譜進行計算,其中順風向湍流功率譜表示為

(4)

式中:σu為順風向的湍流強度;f為風湍流頻率;L為積分尺度參數。L可取為

L=2.45·min(60 m,zhub)

(5)

根據IEC-61400-1選擇B級湍流特性及正常湍流模型(normal turbulence model, NTM)進行計算。NTM中順風向風速分量的標準差代表值σu表達式為

σu=I15(0.75Uhub+b)

(6)

式中:I15為15 m/s風速時湍流強度的期望值,按照B級湍流特性取0.14;b=5.6 m/s。

根據IEC-61400-1,風場的空間相干性考慮為

(7)

式中:r為風場中兩點i與j之間的距離;a為衰減系數,取12;b為偏移系數,可按下式取值

(8)

基于式(3)~式(8),本文采用Turbsim[31]模擬風輪平面的湍流風場。以輪轂為中心,采用覆蓋風輪掃掠面積的風速網格進行風場的模擬,網格上設有31×31個風速點,覆蓋區域尺寸為145 m×145 m,如圖6所示。

圖6 葉素翼型氣動矢量圖

本文聚焦的塔頂前后向響應主要來自風輪的氣動推力Taero,可通過葉素動量理論求取。計算過程中,分別用動量方程和葉素理論建立氣動推力Taero的表達式,使二者相等來進行風速誘導因子參數的迭代計算,最后將誘導因子帶回葉素理論方程得到計算結果,其中氣動推力Taero表達式分別為

(9)

式中:B為葉片數;ω為葉片旋轉角速度;c為弦長;φ為入流角;r為葉素控制體的半徑;v0為風輪上游無窮遠處風速;a和a′分別為軸向和切向誘導因子;Cn和Ct分別為葉素翼型的法向力系數和切向力系數。

可以看出,氣動荷載的計算與葉片的旋轉速度ω和槳距角β(入流角φ與風速局部攻角α的差,見圖6)密切相關,即風電結構動力響應會受到變速、變槳控制的直接影響。本文重點研究前后向響應的控制,因此僅介紹氣動推力Taero。本文在計算氣動推力Taero過程中考慮了普朗特葉尖損失修正、葛朗沃特修正,具體可參考文獻[25]。

考慮兩特定風速條件7 m/s(切入狀態)與18 m/s(額定狀態),采用OpenFAST軟件中的AeroDyn模塊[32]分別計算切入狀態與額定狀態下的氣動推力Taero,如圖7與圖8所示。

(a) 氣動推力時程

(a) 氣動推力時程

可以看出,切入狀態下變速控制發揮作用,氣動推力Taero受到風輪轉速的顯著影響;而額定狀態下以變槳控制為主,氣動推力Taero與槳距角的動態趨勢保持一致。另外,由于工作狀態的不同,不同風速下氣動推力Taero時程曲線的波動方式存在差異,這會導致不同的結構動力響應。從圖7(b)、圖8(b)中可知,3P激勵是氣動推力Taero的重要組成部分。由于空氣動力產生的超低頻荷載會使塔頂發生緩慢的運動;而3P的頻率較高,會顯著影響塔頂加速度。

1.4 對比結果

在相同的氣動推力Taero的條件下,單自由度模型與OpenFAST耦合模型的塔頂加速度頻譜圖對比如圖9所示。從圖9可知,二者的頻譜特征吻合度高,塔架前后向一階頻峰(后簡稱為一階峰)、風輪旋轉3P頻峰(后簡稱為3P峰)的頻譜特征均可對應,等效模型可以較好地反映塔頂前后向響應的頻域特性。

(a)

然而,SDOF模型的3P頻峰偏低,這是由于集中質量忽略了柔性葉片的動力放大作用??紤]到需要重點關注的主頻峰與3P峰在等效模型中(尤其是切入狀態)可以有較好的體現,不會影響振動控制參數研究的規律,故此等效模型可以進行下一步應用。需要說明的是,由于氣動阻尼由風輪旋轉產生,SDOF模型無法精細化考慮風輪的運動特性,阻尼比未按照衰減測試結果嚴格取值,取0.04~0.05之間對應效果最佳。

2 主動質量阻尼器

2.1 控制原理

主動質量阻尼器(active mass damper, AMD)的振動控制應用已經較為成熟,由TMD進一步發展而來,同樣依靠質量體的慣性運動提供控制力。二者相比,TMD的控制性能受到連接剛度(調諧)的限制,減振頻帶較窄。而AMD的特點在于沒有剛度、阻尼元件,不依靠調諧對結構施加相對反向的控制力,直接通過作動器將質量與被控結構連為一體,根據響應反饋進行實時控制。

AMD系統內一般設有傳感器、智能算法與作動器,通過傳感器系統測量結構的干擾和(或)反應,反饋至算法系統并進行最優控制力計算,以此驅動作動器并推動AMD進行慣性質量運動,最終達到減振控制的目的。在風電結構的振動控制領域中,AMD的相關研究幾乎處于空白,考慮到風電結構響應較復雜的頻域特征,無需調諧的AMD具有更高的應用價值。

進一步,受阻尼器控制的塔頂集中質量運動方程表示為

(10)

式中:Taero為風輪表面的氣動荷載;FDamp為阻尼器對結構提供的控制力;ms、cs、ks分別代表結構的整體質量、阻尼系數及剛度系數,可表示為

(11)

cs=2msωsξs

(12)

式中:fs為風電塔的一階頻率;ξs為風電塔的整體阻尼比。

2.2 LQR控制算法

TMD的參數確定一般采用定點理論[33]根據結構一階頻率確定,而AMD控制力需在一定智能算法和反饋機制的基礎上進行實時計算。

LQR是工程界常用的控制理論之一,其對象是現代控制理論中以狀態空間形式給出的線性系統,而目標函數為對象狀態和控制輸入的二次型函數。LQR最優設計的原理是設計出狀態反饋控制器G使得二次型目標函數J取最小值。狀態反饋控制器G又稱為反饋增益矩陣,由權矩陣Q和R決定,故Q、R矩陣的選擇是使用LQR進行系統控制的關鍵。

根據式(10),使用狀態空間方程描述風電塔-AMD等效模型以方便建模及觀測,受控系統可表示為

(13)

(14)

(15)

(16)

系統的二次型性能指標J表示為

(17)

式中:Q為半正定對稱權矩陣;R為正定對稱權矩陣,二者用以平衡系統響應與控制力輸入的權重。

為使性能指標J最小,需設計反饋增益矩陣G以獲得最優驅動力

U=-GZ

(18)

式中,反饋增益矩陣G可進一步表示為

G=R-1BTP

(19)

式中,P為萊卡提(Riccati)矩陣函數,是Riccati方程的解

PA+ATP-PBR-1BTP+Q=0

(20)

由此推算反饋增益矩陣G,將式(18)代入式(13),結構-AMD控制系統的狀態方程可進一步表示為

(21)

可知,對比TMD,LQR控制相當于通過反饋增益矩陣G改變系統矩陣A,此時(A-BG)為等效系統矩陣。另外可知,反饋增益矩陣G只由系統自身性質(A,B,Q,R)決定,不受外荷載DT的影響。

在現代控制理論中,權矩陣Q與R通常需要考慮狀態空間方程的特點,根據經驗進行取值;對于本文的控制問題來說,權矩陣Q與R通常取為

(22)

R=β2

(23)

式中:K與M分別為AMD-結構體系的剛度矩陣與質量矩陣;β1與β2分別為權系數,需要根據實際結構進行確定。一般來說,Q越大,受控結構反應越小,控制效果越好;R越大,則控制輸入越小,控制效果越差。

進一步地,根據上述表達式采用SIMULINK平臺建立AMD-風電塔體系模型,系統矩陣A、驅動力輸入矩陣B和荷載輸入矩陣D分別在MATLAB的workspace中定義,反饋增益矩陣G通過lqr函數進行求取,并通過sim函數賦值給SIMULINK模型。

3 減振分析

對AMD進行參數設計并設計一組相同質量比μ=0.05[34-35]的TMD作為對照組以探究AMD的減振性能與控制特點。其中,質量比μ為阻尼器質量md與機艙、葉片等風機頂部集中質量ms的比值;AMD的權系數為經驗取值,TMD的頻率比γ及阻尼比ξd通過定點理論計算得出,設計參數如表2所示。

表2 阻尼器參數

3.1 控制效果

綜合以上模型與荷載理論,進行AMD減振控制計算,切入與額定狀態下的時程曲線分別如圖10與圖11所示。從圖10和11可知,RNA運行作用下,AMD可以抑制風電塔頂加速度響應;振動響應的幅度與峰值都獲得了較好地控制,AMD的控制作用得到了驗證。

圖10 切入狀態加速度響應對比

圖11 額定狀態加速度響應時程對比

3.2 控制對比

通過對比AMD與TMD的減振響應,可從時域和頻域兩方面考察二者對于加速度響應的控制作用。為方便量化二者對時域響應的控制性能,引入加速度減振率R

(24)

式中:anone與adamp分別泛指未受控與受控結構的加速度指標,如加速度峰值與加速度均方根。據此,可將減振率R分別進一步定義為峰值減振率Rp與均方根減振率Rr。

切入與額定狀態下,TMD與AMD的減振性能分別如圖12、圖13所示。整體來看,AMD對于響應的控制作用都優于TMD:在切入狀態下,TMD幾乎無法有效控制塔頂響應,而AMD的峰值與均方根減振率分別達到9.05%與13.22%;在額定狀態下,TMD的控制效果較切入時有了改善但其減振效率依舊低于AMD,此時AMD的峰值與均方根減振率分別達到11.13%與17.07%。對比AMD在風機不同工作狀態下的表現,可以看出AMD在額定狀態時的減振性能略高于切入狀態。

圖12 切入狀態加速度時程對比

圖13 額定狀態加速度時程對比

切入與額定狀態下,TMD與AMD的減振頻譜圖分別如圖14、圖15所示。從圖14和圖15可知,TMD與AMD的主要控制頻段均位于塔架一階頻率附近,不同之處在于AMD的減振頻帶更寬,甚至會對3P頻帶起到一定控制作用。究其根本,從式(21)中可知,AMD實質上改善了風電結構的系統矩陣A,由于該矩陣主要信息為結構質量、一階頻率及阻尼比,缺乏風輪轉速的信息,因此在控制時也更側重對塔一階頻率及相鄰頻段的控制。另外,AMD可以影響3P頻帶的原因在于LQR算法的反饋增益會根據響應反饋控制力。

(a) AMD

(a) AMD

前文提到,AMD在額定狀態的減振性能略高于切入狀態,可以在頻譜圖中得到解釋:切入狀態下,在低風速的變速控制中,塔頂響應頻譜的3P峰更加突出,塔一階振動在氣動阻尼的影響下較為微弱,AMD發揮的空間比較有限;而在額定狀態,塔一階峰較切入狀態更突出,因此也會受到更好的抑制。這一特點是由AMD的反饋機制決定的。

3.3 性能曲線

RNA運行作用下,風電結構在不同狀態中具有不同的動力響應特點,而AMD在不同風機狀態下也有不同的表現??紤]到風機實際運行中,風速是隨機不可控的,因此探究AMD在全運行風況下的控制性能是有必要的。本文考慮了3~25 m/s的平均風速,每組平均風速中考慮了10組隨機湍流風。通過批處理程序計算質量比μ=0.05的AMD與TMD在每種風速條件下的均方根減振率Rr,如圖16所示。

圖16 TMD與AMD的性能曲線對比

從圖16可知,AMD的減振性能隨風速的增加先降低后提升,拐點風速會低于額定風速(11.4 m/s),這是由于當平均風速接近額定風速時,湍流風的瞬時風速會在某些時刻達到額定風速使風機進入額定狀態;結合前文的分析,此工況下AMD的減振性能會有所提升,故出現拐點。相較之下,TMD的減振性能隨風速的增加而持續上升。對比可知,AMD在任何風況下的控制性能都優于TMD ,減振率Rr高出超10%。

切入狀態減振性能隨風速增加而降低的原因在于:當風速較低時,塔一階峰與3P峰的距離較近(塔一階頻率保持0.278 Hz,3 m/s的3P頻率約為0.35 Hz,7 m/s的3P頻率約為0.42 Hz),AMD在控制塔一階峰的同時更容易輻射到3P頻段的幅值;隨著風速增加,兩峰值距離越來越遠,AMD側重于塔一階峰的抑制,對3P峰的影響降低。額定狀態下,AMD減振性能隨風速增加而增加的原因在于氣動荷載越來越小,塔頂響應中的氣動成分越來越低,塔一階振動越來越劇烈,AMD的發揮空間逐漸變大??傮w來說,作為風電結構RNA運行狀態下的控制方案,AMD具有較可觀的減振性能。

4 結 論

本文考慮風機RNA運行的變速、變槳控制,基于單自由度等效風電塔模型,研究了AMD對塔頂前后向響應的控制作用,通過與TMD的對比考察了AMD的減振性能與控制特點,得到以下結論:

(1) 通過考慮RNA運行作用下的氣動阻尼與氣動推力,單自由度等效模型可以較好地模擬風電塔前后向的響應頻譜特征。

(2) 基于LQR控制算法的AMD可以較好地抑制塔頂前后向加速度響應,其控制性能顯著超過TMD。

(3) AMD的本質在于改善受控風電塔的動力特性,其主要控制頻段依舊為塔一階頻率及相鄰頻段,但由于主動反饋機制,AMD的控制作用會輻射到3P頻段。

(4) 由于3P轉頻的時變性,AMD的減振性能隨風速的增加先降低后提升,拐點位于額定風速(11.4 m/s)處,在風機的全運行風況內都有較好的控制性能。

猜你喜歡
風速結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
3D打印中的模型分割與打包
基于GARCH的短時風速預測方法
主站蜘蛛池模板: 日韩在线网址| 国产三级成人| 亚洲精品男人天堂| 欧美三级自拍| 四虎国产精品永久一区| 天堂成人在线视频| 一本一道波多野结衣一区二区| 色哟哟国产精品| 国产一级视频久久| julia中文字幕久久亚洲| 女人爽到高潮免费视频大全| 亚洲第一视频免费在线| 特级精品毛片免费观看| 国产成在线观看免费视频 | 亚洲日本一本dvd高清| 亚洲国产高清精品线久久| 暴力调教一区二区三区| 亚洲av无码成人专区| 中文国产成人精品久久一| 69av在线| 色综合五月| 欧美不卡视频在线| 男女精品视频| 老司机午夜精品网站在线观看 | 国产欧美日韩在线一区| 国产系列在线| 一本综合久久| 亚洲午夜综合网| 亚洲天堂视频在线播放| 久久综合伊人77777| 曰韩免费无码AV一区二区| 亚洲欧美在线综合一区二区三区| 国产久草视频| 伊大人香蕉久久网欧美| 成人福利视频网| 精品视频第一页| 国产91无毒不卡在线观看| 亚洲有无码中文网| 色网站在线免费观看| 国产福利观看| 久久伊人久久亚洲综合| 少妇精品网站| 丁香六月综合网| 国产一级毛片网站| 国产成人午夜福利免费无码r| 免费无码AV片在线观看中文| 91综合色区亚洲熟妇p| 暴力调教一区二区三区| 亚洲中文字幕97久久精品少妇| 国产香蕉国产精品偷在线观看| 久久网综合| 久久免费观看视频| 日韩久久精品无码aV| 亚洲高清中文字幕| 毛片免费在线视频| 久久精品嫩草研究院| 香蕉国产精品视频| 欧美性久久久久| 欧美啪啪视频免码| 99久久精品国产自免费| 国产欧美日韩在线一区| 在线看片中文字幕| 91福利一区二区三区| 亚洲成人77777| 亚洲色欲色欲www网| AV熟女乱| 欧美黄色网站在线看| 伊人蕉久影院| 日韩欧美高清视频| 日本黄色a视频| 美美女高清毛片视频免费观看| 亚洲一区二区三区麻豆| 国产精品福利导航| 国产精品内射视频| 午夜视频免费一区二区在线看| 精品福利一区二区免费视频| 亚洲最猛黑人xxxx黑人猛交| 亚洲视频在线青青| 亚洲无码视频一区二区三区| 亚洲人成网站在线观看播放不卡| 国产免费黄| 中文字幕首页系列人妻|