• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    海上浮式風(fēng)機(jī)的一體化建模及其整體可靠性

    2022-08-23 06:50:58陳建兵宋玉鵬
    空氣動力學(xué)學(xué)報 2022年4期
    關(guān)鍵詞:結(jié)構(gòu)分析方法

    陳建兵,宋玉鵬

    (1. 同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室,上海 200092;2. 同濟(jì)大學(xué) 土木工程學(xué)院,上海 200092;3. 南京工業(yè)大學(xué) 土木工程學(xué)院,南京 211816)

    0 引言

    深遠(yuǎn)海風(fēng)能具有資源豐富、風(fēng)力穩(wěn)定、環(huán)境友好等諸多優(yōu)點(diǎn),近年來在全球范圍內(nèi)獲得了廣泛關(guān)注和積極發(fā)展[1]。為了降低深遠(yuǎn)海風(fēng)能開發(fā)成本,發(fā)電機(jī)組采用漂浮式基礎(chǔ)結(jié)構(gòu)已成為共識[2]。當(dāng)前,深遠(yuǎn)海風(fēng)能開發(fā)主要面臨機(jī)組服役條件十分惡劣和由于離岸較遠(yuǎn)而造成高昂維修保養(yǎng)費(fèi)用等挑戰(zhàn)[3-4]。此外,與陸上風(fēng)能和近海風(fēng)能開發(fā)相比,深遠(yuǎn)海風(fēng)能開發(fā)中支撐結(jié)構(gòu)的費(fèi)用急劇增長[1]。上述因素使得當(dāng)前浮式風(fēng)機(jī)尚未完全實(shí)現(xiàn)商業(yè)化,限制了深遠(yuǎn)海風(fēng)能的大規(guī)模開發(fā)。為此,應(yīng)當(dāng)采用基于整體可靠度的結(jié)構(gòu)設(shè)計方法,實(shí)現(xiàn)兼顧浮式風(fēng)機(jī)結(jié)構(gòu)安全性與經(jīng)濟(jì)性的定量設(shè)計,為深遠(yuǎn)海風(fēng)能的大規(guī)模開發(fā)提供關(guān)鍵技術(shù)保障。因此,浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析具有重要意義。

    當(dāng)前,復(fù)雜結(jié)構(gòu)的整體可靠度評估通常采用數(shù)值分析方法,主要涉及三個方面的問題:結(jié)構(gòu)響應(yīng)分析的數(shù)值模型、隨機(jī)荷載的建模、結(jié)構(gòu)可靠度分析方法[5]。浮式風(fēng)機(jī)結(jié)構(gòu)的精細(xì)化數(shù)值分析模型是其整體可靠度分析的關(guān)鍵基礎(chǔ)[6]。由于浮式風(fēng)機(jī)結(jié)構(gòu)中同時存在多種復(fù)雜耦合效應(yīng)[7-8],例如槳葉-機(jī)艙-塔筒-浮體-錨索等部件之間的相互作用、剛體運(yùn)動與彈性變形耦合、氣動-結(jié)構(gòu)-機(jī)電-水動耦合、槳葉轉(zhuǎn)動效應(yīng)和幾何非線性效應(yīng)等,使得浮式風(fēng)機(jī)結(jié)構(gòu)的一體化動力學(xué)建模成為一項(xiàng)難題[9]。雖然針對該建模問題當(dāng)前已有較多研究,但是在實(shí)際工程中仍主要采用分離式設(shè)計方法,難以合理考慮上述復(fù)雜的耦合效應(yīng),也難以進(jìn)行結(jié)構(gòu)的整體可靠度分析。一般來說,浮式風(fēng)機(jī)結(jié)構(gòu)的一體化建模主要基于多體動力學(xué)理論[9],采用Lagrange方法或者Kane方法建立結(jié)構(gòu)的整體動力學(xué)方程,并可以進(jìn)一步分為多剛體動力學(xué)模型和多柔性體動力學(xué)模型[10]。多剛體動力學(xué)模型將浮式風(fēng)機(jī)中的所有部件均考慮為剛體,因而無法獲得柔性部件(槳葉、塔筒等)的彈性變形和內(nèi)力響應(yīng),具有較大的局限性[10]。但是,由于多剛體動力學(xué)模型建模簡便、分析效率較高,因此多用于浮式風(fēng)機(jī)結(jié)構(gòu)的概念設(shè)計和初始設(shè)計階段[11]。多柔性體模型能夠同時考慮系統(tǒng)中的剛體運(yùn)動和彈性變形響應(yīng),其中一般采用有限元方法或者模態(tài)疊加方法建立柔性部件的彈性模型[7]。雖然多柔性體模型能夠全面反映系統(tǒng)的力學(xué)特性和復(fù)雜的耦合效應(yīng),但是該類模型建模較為復(fù)雜,同時存在分析效率不高的問題,難以應(yīng)用于需要進(jìn)行大量數(shù)值模擬分析的問題中,例如結(jié)構(gòu)的可靠度評估與優(yōu)化設(shè)計等[12]。

    浮式風(fēng)機(jī)在服役期間的主要環(huán)境荷載為風(fēng)荷載和波浪荷載,且兩種荷載具有相關(guān)性[13]。因此,在浮式風(fēng)機(jī)結(jié)構(gòu)的動力分析與可靠度評估中應(yīng)當(dāng)合理考慮風(fēng)和波浪荷載的聯(lián)合作用[14]。實(shí)測數(shù)據(jù)表明,10 min內(nèi)的風(fēng)速時程和20 min至3~6 h內(nèi)的波浪時程均可以認(rèn)為是平穩(wěn)過程,并可以采用平均風(fēng)速、有效波高、譜峰周期、風(fēng)-浪夾角等物理量描述上述風(fēng)況和海況[13]。為了刻畫結(jié)構(gòu)在服役期內(nèi)所經(jīng)歷的所有工況的概率大小,可基于結(jié)構(gòu)所在場址的長期觀測數(shù)據(jù)建立上述物理量的聯(lián)合概率分布,這也是當(dāng)前海洋結(jié)構(gòu)設(shè)計規(guī)范中所推薦的方法[13]。在建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型時,可采用條件概率分布建模方法[15]、Nataf變換方法[16]和Copula方法[17-18]。條件概率分布建模方法是當(dāng)前海洋結(jié)構(gòu)設(shè)計規(guī)范中推薦使用的方法[13]。該方法首先確定主導(dǎo)隨機(jī)變量,并建立其邊緣概率分布模型,然后根據(jù)主導(dǎo)隨機(jī)變量的取值將數(shù)據(jù)劃分為一系列區(qū)間,進(jìn)而估計主導(dǎo)隨機(jī)變量取不同值條件下的其他隨機(jī)變量的概率分布。但是,當(dāng)隨機(jī)變量的個數(shù)較多時,該方法需要進(jìn)行多次條件概率分布建模。此外,在某些條件下(例如高風(fēng)速條件),觀測數(shù)據(jù)較少,將導(dǎo)致該條件下的概率分布建模結(jié)果誤差較大[18]。Nataf變換方法首先基于等概率變換原則,將原始隨機(jī)變量變換為服從標(biāo)準(zhǔn)高斯分布的隨機(jī)變量,進(jìn)而利用多維高斯分布模型建立變換后的隨機(jī)變量的聯(lián)合概率分布,最后再利用反變換獲得原始隨機(jī)變量的聯(lián)合概率分布[16]。該方法雖然簡便易行,但是當(dāng)隨機(jī)變量為非線性、非高斯相關(guān)時,例如風(fēng)-浪參數(shù)間的相關(guān)結(jié)構(gòu)通常呈現(xiàn)這種特點(diǎn),該方法將引入較大誤差。Copula方法是近年來發(fā)展的多維隨機(jī)變量聯(lián)合概率分布建模的新方法,對于非線性相關(guān)問題十分有效[17]。該方法在聯(lián)合概率分布建模中,將隨機(jī)變量的邊緣概率分布建模與變量間的概率相關(guān)結(jié)構(gòu)建模分開處理,具有簡便、靈活、精確等優(yōu)點(diǎn)。該方法目前已在多個領(lǐng)域獲得廣泛關(guān)注[19],但在海洋工程中尚應(yīng)用較少。

    結(jié)構(gòu)可靠度分析方法主要包括一次和二次可靠度方法、基于跨越過程理論的首次超越破壞可靠度方法、隨機(jī)模擬方法、概率密度演化理論等[5,20]。一次和二次可靠度方法常用于靜力問題的分析,具有較高的效率,然而該方法需要將系統(tǒng)的極限狀態(tài)函數(shù)展開為一階或二階形式,對于強(qiáng)非線性系統(tǒng)(例如浮式風(fēng)機(jī)結(jié)構(gòu))會產(chǎn)生較大誤差[9]?;诳缭竭^程理論的首次超越破壞可靠度方法適用于動力系統(tǒng),但通常引入Poisson假定,具有一定的近似性,尤其對于大失效概率問題其誤差較大[5]。隨機(jī)模擬方法主要包括Monte Carlo模擬及其改進(jìn)方法以及子集模擬等方法,這類方法不受問題的維度、非線性程度等因素的限制,應(yīng)用范圍較廣,但通常具有隨機(jī)收斂和效率不高等問題[20]。為了將這類方法應(yīng)用于浮式風(fēng)機(jī)的可靠度分析中,通常發(fā)展浮式風(fēng)機(jī)結(jié)構(gòu)分析的代理模型以提高效率,但同時也引入了新的誤差[21]。概率密度演化理論給出了隨機(jī)動力系統(tǒng)響應(yīng)的概率密度函數(shù)的演化方程,揭示了系統(tǒng)響應(yīng)隨機(jī)性演化的內(nèi)在物理規(guī)律[5],該方法在分析系統(tǒng)的可靠度時同樣不受問題類型、問題維度和非線性程度等因素的影響,能給出系統(tǒng)響應(yīng)量的概率密度函數(shù)演化過程,同時具有較高的分析效率,因此適用于浮式風(fēng)機(jī)的可靠度分析,但目前在此領(lǐng)域的應(yīng)用研究較少。

    為了分析浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度,首先引入結(jié)合多剛體動力學(xué)理論和有限元方法建立起來的浮式風(fēng)機(jī)結(jié)構(gòu)一體化分析模型(StoDRAFOWT模型)。該模型不僅能同時獲得浮式風(fēng)機(jī)的剛體運(yùn)動和彈性變形響應(yīng),而且具有建模簡便、易于求解等優(yōu)點(diǎn)[9]。然后,基于我國南海某場址長期風(fēng)-浪數(shù)據(jù)并引入Copula方法,建立風(fēng)-浪聯(lián)合作用的多參數(shù)非線性相關(guān)的聯(lián)合概率分布模型,進(jìn)而采用概率密度演化理論分析浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度。由此,建立了風(fēng)-浪聯(lián)合作用下浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠度分析的基本框架。針對風(fēng)機(jī)結(jié)構(gòu)設(shè)計中需要考慮的額定風(fēng)速和極端風(fēng)剪不利工況,采用上述方法進(jìn)行了單柱式浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析。

    1 浮式風(fēng)機(jī)一體化動力學(xué)建模

    1.1 基本建模思想與坐標(biāo)系統(tǒng)

    如前所述,浮式風(fēng)機(jī)的多剛體動力學(xué)模型雖然分析效率較高,但是無法分析結(jié)構(gòu)的彈性變形。而多柔性體動力學(xué)模型雖然能夠全面分析結(jié)構(gòu)的動力響應(yīng),但是建模復(fù)雜、分析效率不高。本文首先介紹將多剛體動力學(xué)模型與有限元方法相結(jié)合的建模策略,不僅能夠全面分析結(jié)構(gòu)響應(yīng),而且建模較為簡便,分析效率較高[9]。已有研究表明[10-11],對于當(dāng)前兆瓦級的浮式風(fēng)機(jī)系統(tǒng),結(jié)構(gòu)的彈性變形對系統(tǒng)剛體運(yùn)動響應(yīng)的影響可以忽略,因此可以首先采用多剛體動力學(xué)模型分析系統(tǒng)的剛體運(yùn)動響應(yīng)。而系統(tǒng)的剛體運(yùn)動對結(jié)構(gòu)的彈性變形具有重要影響。為此,采用Kane方法[22-23]推導(dǎo)了考慮系統(tǒng)剛體運(yùn)動效應(yīng)的結(jié)構(gòu)有限元模型,從而能夠全面獲得結(jié)構(gòu)的動力響應(yīng)。

    圖1為單柱式浮式風(fēng)機(jī)的示意圖和坐標(biāo)系的定義。為描述浮式風(fēng)機(jī)系統(tǒng)的狀態(tài),定義三組笛卡爾坐標(biāo)系,分別為慣性坐標(biāo)系Ox0y0z0、浮體-塔筒的局部坐標(biāo)系Oxyz、 槳葉的局部坐標(biāo)系Ox′y′z′。 其中Ox0y0z0坐標(biāo)系的原點(diǎn)位于海平面位置,x0軸與風(fēng)速方向一致,z0軸豎直向上。在初始時刻浮式風(fēng)機(jī)未發(fā)生位移時,坐標(biāo)系Oxyz與Ox0y0z0重 合。坐標(biāo)系Ox′y′z′的原點(diǎn)位于輪轂中心,x′軸與x軸 平行,z′軸由葉根指向葉尖。

    圖1 單柱式浮式風(fēng)機(jī)示意圖與坐標(biāo)系統(tǒng)Fig. 1 Schematic diagram of a spar-type floating offshore wind turbine and the coordinate system

    1.2 一體化結(jié)構(gòu)分析模型

    為了建立浮式風(fēng)機(jī)系統(tǒng)的多剛體動力學(xué)模型,通?;贚agrange方程得到考慮轉(zhuǎn)子轉(zhuǎn)動和浮體六自由度剛體運(yùn)動的動力方程[10],即包含七個自由度。由于現(xiàn)代風(fēng)力機(jī)通常配備槳距控制系統(tǒng),為此,將變槳距控制方程與上述七自由度動力方程耦合,得到考慮變速、變槳距控制的浮式風(fēng)機(jī)多剛體動力學(xué)模型[9,24]:

    式中,M0和C0分別為浮式風(fēng)機(jī)的質(zhì)量矩陣和阻尼矩陣,其表達(dá)式可參見文獻(xiàn)[24];Irx為轉(zhuǎn)子相對于x′軸的轉(zhuǎn)動慣量;vp和 ωp=(ωx,ωy,ωz)T分 別為浮體在Oxyz坐標(biāo)系中的平動和轉(zhuǎn)動速度; ωb和 φ分別為槳葉的轉(zhuǎn)動速度和轉(zhuǎn)動角度;β為槳葉的槳距角;c0、G、 τd和 τp是與風(fēng)力機(jī)槳距控制和轉(zhuǎn)矩控制有關(guān)的參數(shù); ω0為低轉(zhuǎn)速軸的額定轉(zhuǎn)速;FG為重力在Oxyz坐標(biāo)系下的效應(yīng);MG為 在Oxyz坐標(biāo)系下重力關(guān)于坐標(biāo)系原點(diǎn)的力矩;FD為 外部荷載在Oxyz坐 標(biāo)系下的效應(yīng);MD為在Oxyz坐標(biāo)系下外部荷載關(guān)于坐標(biāo)系原點(diǎn)的力矩;FM為低轉(zhuǎn)速軸上的凈扭矩。

    通過求解方程式(1)即可獲得浮式風(fēng)機(jī)中浮體的六自由度剛體運(yùn)動響應(yīng)、槳葉轉(zhuǎn)速和槳葉的槳距角。為了分析浮式風(fēng)機(jī)中柔性部件的彈性變形,基于Kane方程分別推導(dǎo)得到塔筒單元和槳葉單元的有限元表達(dá)式。

    在坐標(biāo)系Oxyz中,可以得到考慮剛體運(yùn)動效應(yīng)的塔筒單元的動力方程[9]:

    式中:Mt和 μ分別為單元的質(zhì)量矩陣和結(jié)點(diǎn)位移向量;Ct、Kt和Ft分別為單元的剛度矩陣、阻尼矩陣和荷載向量,表達(dá)式為:

    式中:ρ為單元的材料密度;V為單元體積;N?為形函數(shù)矩陣;上標(biāo)“T”表示轉(zhuǎn)置矩陣; ω?為浮體轉(zhuǎn)動速度ωp對 應(yīng)的反對稱矩陣[24];R為坐標(biāo)轉(zhuǎn)換矩陣,將坐標(biāo)系Oxyz中的矢量轉(zhuǎn)換至坐標(biāo)系Ox0y0z0中 ;ro為Oxyz的坐標(biāo)系原點(diǎn)在坐標(biāo)系Ox0y0z0中 的位矢;rk為單元中任意一點(diǎn)在坐標(biāo)系Oxyz中的位矢。

    類似地,在坐標(biāo)系Ox′y′z′中,可以得到考慮剛體運(yùn)動和旋轉(zhuǎn)效應(yīng)的槳葉單元的動力方程[9]:

    式中:Mb為 單元的質(zhì)量矩陣;Cb、Kb和Fb分別為單元的剛度矩陣、阻尼矩陣和荷載向量,其表達(dá)式與式(3)、式(4)、式(5)類似。其中,由于剛體運(yùn)動和旋轉(zhuǎn)效應(yīng)產(chǎn)生的槳葉單元的附加阻尼矩陣、附加剛度矩陣和附加荷載向量的表達(dá)式為:

    其中:T為坐標(biāo)轉(zhuǎn)換矩陣,將坐標(biāo)系Ox′y′z′中的矢量轉(zhuǎn)換至坐標(biāo)系Oxyz中;為 槳葉轉(zhuǎn)動速度 ( ωb,0,0)T對應(yīng)的 反 對 稱 矩 陣; ¨ro為Oxyz的 坐 標(biāo) 系 原 點(diǎn) 在 坐 標(biāo) 系Ox0y0z0中 的加速度;ro′為Ox′y′z′的坐標(biāo)系原點(diǎn)在坐標(biāo)系Oxyz中的位矢;rq為 單元中任意一點(diǎn)在坐標(biāo)系Ox′y′z′中的位矢。

    值得注意的是,槳葉截面為非對稱形狀,并且槳葉截面存在沿翼展方向變化的預(yù)扭轉(zhuǎn)角,同時還存在槳距控制。因此,槳葉在風(fēng)輪平面內(nèi)、外兩個方向的振動相互耦合,在建模過程中應(yīng)當(dāng)合理考慮。此外,槳葉在風(fēng)機(jī)正常運(yùn)轉(zhuǎn)工況下的變形較大,應(yīng)當(dāng)在槳葉的三維梁單元模型中考慮幾何非線性效應(yīng)。這是與塔筒有限元模型的不同之處。

    基于塔筒和槳葉的有限元模型,并結(jié)合前述浮式風(fēng)機(jī)的剛體運(yùn)動響應(yīng),可建立浮式風(fēng)機(jī)柔性構(gòu)件的整體有限元模型,進(jìn)而分析浮式風(fēng)機(jī)結(jié)構(gòu)彈性變形響應(yīng)。

    1.3 荷載計算

    在浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)分析中,主要應(yīng)考慮槳葉和塔筒的氣動荷載、浮體的水動力荷載和錨泊荷載。本文采用葉素動量理論方法[25]計算槳葉的氣動荷載,并考慮了Prandtl葉尖損失因子和高推力系數(shù)條件下的Glauert修正[25]。塔筒的氣動荷載可基于準(zhǔn)定常假定計算[25]。

    在計算大型浮體的水動力荷載時,應(yīng)當(dāng)合理考慮輻射效應(yīng)和繞射效應(yīng)[26]。本文關(guān)注的是單柱式浮式風(fēng)機(jī),浮體的截面尺寸相對較小,因此可以簡化考慮這兩種效應(yīng)。在實(shí)際中可采用相對形式的Morison公式計算單柱式浮體的水動力荷載[10,26],此時,單位長度浮體的水動力荷載為[27]:

    式中:ρw為 海水密度;dp為浮體外徑;Cd和CA分別為拖曳力系數(shù)和附加質(zhì)量系數(shù);vf和af分別為該浮體微元的速度和加速度向量;vw和aw分別為水質(zhì)點(diǎn)的速度和加速度,可采用Airy波浪理論、Stokes波浪理論或其他波浪理論 計 算[13,25]; ∥·∥表示 向 量 的幅值。應(yīng) 指 出,求解結(jié)構(gòu)的動力方程(式1)時,通常將式(13)中的移至動力方程的左邊,由此產(chǎn)生浮體的附加質(zhì)量矩陣。

    在計算浮式風(fēng)機(jī)的錨泊荷載時,需要建立錨索的力學(xué)分析模型,主要包括錨索的靜力模型、擬靜力模型和動力模型[28]。由于靜力和擬靜力模型在某些工況下會低估錨索的疲勞荷載,因此推薦采用錨索的動力模型[28]。本文采用錨索的集中質(zhì)量動力模型[28]計算錨泊荷載。該模型不僅考慮錨索的受拉彈性變形,而且將錨索的動力效應(yīng)和水動力效應(yīng)以及與海底的接觸效應(yīng)等也考慮在內(nèi)。錨索底部固定在海床上,錨索頂部的邊界條件由浮體的位置確定。通過求解錨索的動力方程,可以獲得錨索對浮體的荷載效應(yīng)。

    1.4 模型驗(yàn)證

    為了檢驗(yàn)上述模型的正確性,分別利用該模型和當(dāng)前風(fēng)力機(jī)結(jié)構(gòu)分析的常用軟件FAST,分析風(fēng)-浪聯(lián)合作用下的海上單柱式浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng),并進(jìn)行對比。浮式風(fēng)機(jī)結(jié)構(gòu)參數(shù)采用美國國家可再生能源實(shí)驗(yàn)室(NREL)給出的5 MW海上風(fēng)力發(fā)電機(jī)組參數(shù)[29]和OC3項(xiàng)目中海上單柱式浮式風(fēng)機(jī)基礎(chǔ)參數(shù)[30],關(guān)鍵參數(shù)如表1所示。

    表1 浮式風(fēng)力發(fā)電機(jī)組關(guān)鍵參數(shù)Table 1 Key parameters of the floating offshore wind turbine

    本文僅給出一種工況條件下的分析對比結(jié)果,其他工況條件下的對比結(jié)果見文獻(xiàn)[24]。此處的工況條件為:定常風(fēng)速場,且輪轂高度處的風(fēng)速為15 m/s;規(guī)則波浪,波高為6 m,波浪周期為10 s。分析結(jié)果如圖2所示。結(jié)果對比表明,本文模型的分析結(jié)果與FAST的分析結(jié)果吻合良好,尤其是在響應(yīng)進(jìn)入平穩(wěn)階段以后。因此,可以利用本文模型進(jìn)行浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度分析。

    圖2 風(fēng)-浪作用下浮式風(fēng)機(jī)結(jié)構(gòu)響應(yīng)的對比結(jié)果Fig. 2 Comparisons of the responses of floating offshore wind turbine under wind and wave excitations

    2 風(fēng)-浪聯(lián)合作用的概率模型

    由于浮式風(fēng)機(jī)在服役期間同時受到風(fēng)荷載和波浪荷載作用,因此在進(jìn)行浮式風(fēng)機(jī)的動力響應(yīng)與可靠度分析時,應(yīng)當(dāng)合理考慮風(fēng)-浪聯(lián)合作用。如前所述,當(dāng)前海洋結(jié)構(gòu)的設(shè)計規(guī)范中建議采用風(fēng)-浪參數(shù)的聯(lián)合概率分布模型考慮風(fēng)-浪聯(lián)合作用[13]。本文采用Copula方法建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型。

    在實(shí)際工程中,通常需要考慮平均風(fēng)速Vw、有效波高Hs、 波浪譜峰周期Tp以及風(fēng)-浪夾角等物理量。本文主要關(guān)注額定風(fēng)速條件下發(fā)生極端風(fēng)剪時浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度。根據(jù)風(fēng)力機(jī)結(jié)構(gòu)的設(shè)計規(guī)范,在分析該工況下風(fēng)力機(jī)結(jié)構(gòu)的動力響應(yīng)時,取風(fēng)-浪夾角為0°[31]。因此,需要建立額定風(fēng)速條件下Hs和Tp的條件聯(lián)合概率分布模型。

    在建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型時,應(yīng)當(dāng)基于場址的長期風(fēng)-浪觀測數(shù)據(jù)。本文采用歐洲中尺度天氣預(yù)報中心ERA-Interim數(shù)據(jù)庫[32]提供的風(fēng)-浪參數(shù)再分析數(shù)據(jù),獲取了我國南海場址(21°N,113°E)處1979年至2018年(共計40年)的風(fēng)-浪再分析數(shù)據(jù)。這些數(shù)據(jù)均為30 min內(nèi)觀測數(shù)據(jù)的平均值,時間分辨率為6 h。通過對數(shù)據(jù)預(yù)處理,可得到輪轂處平均風(fēng)速、有效波高和譜峰周期的長期數(shù)據(jù)。

    為了建立額定風(fēng)速條件下Hs和Tp的條件聯(lián)合概率分布模型,在上述風(fēng)-浪數(shù)據(jù)的基礎(chǔ)上選取額定風(fēng)速附近的數(shù)據(jù)進(jìn)行概率建模。本文中,風(fēng)力機(jī)的額定風(fēng)速為11.4 m/s,因此取風(fēng)速在11.3 m/s ~11.5 m/s區(qū)間內(nèi)的風(fēng)-浪數(shù)據(jù)進(jìn)行Hs和Tp的聯(lián)合概率分布建模,共有681個數(shù)據(jù)點(diǎn)。

    根據(jù)Copula聯(lián)合概率分布建模的基本理論[17],對于任意兩個隨機(jī)變量X1和X2,其聯(lián)合概率密度函數(shù)的表達(dá)式為:

    式中:F1(x1) 和F2(x2) 分 別 為X1和X2的邊緣概率分布函數(shù);f1(x1)和f2(x2)為相應(yīng)的邊緣概率密度函數(shù);c12(·,·)為二元Copula密度函數(shù),可從常用的Copula函數(shù)中選用,例如高斯Copula、t Copula、Clayton Copula、Gumbel Copula、Frank Copula等[17]。

    從式(14)中可以看到,在利用Copula方法建立多維隨機(jī)變量的聯(lián)合概率分布模型時,需要分別建立各隨機(jī)變量的邊緣概率分布模型及變量間的相關(guān)結(jié)構(gòu)模型。由于邊緣概率分布模型及變量間的相關(guān)結(jié)構(gòu)模型可以分開處理,因此該方法具有靈活、便捷、精確等優(yōu)點(diǎn)。

    基于上述681個額定風(fēng)速條件下的數(shù)據(jù)點(diǎn),首先采用常用的參數(shù)概率分布模型,建立Hs和Tp的邊緣概率分布模型,并進(jìn)行假設(shè)檢驗(yàn),進(jìn)而基于AIC準(zhǔn)則[19]確定最優(yōu)的邊緣概率分布模型。結(jié)果表明,Hs和Tp的最優(yōu)邊緣分布均為Gamma分布,結(jié)果見表2和圖3。表中, Γ (·)為Gamma函數(shù)。

    圖3 額定風(fēng)速下有效波高和譜峰周期的邊緣概率分布Fig. 3 Marginal probability distributions of Hs and Tp under the rated wind speed condition

    進(jìn)一步,采用常用的Copula函數(shù)建立Hs和Tp的相關(guān)結(jié)構(gòu)模型。備選的Copula函數(shù)除了上述常用的Copula函數(shù)外,還包括按照如下方式構(gòu)建的非對稱Copula函數(shù)[33]:

    在確定備選的Copula函數(shù)后,采用極大似然方法估算Copula函數(shù)中的待定參數(shù),然后根據(jù)AIC準(zhǔn)則選取最優(yōu)Copula函數(shù)并進(jìn)行擬合優(yōu)度檢驗(yàn)[17,19]。對于上述額定風(fēng)速條件下的風(fēng)浪數(shù)據(jù),識別的最優(yōu)結(jié)果為按照式(15)構(gòu)造的非對稱高斯Copula模型,見表2。表中,C(u,v) 表示高斯Copula, Φλ(·)為相關(guān)系數(shù)為 λ(? 1≤λ≤1) 的二維標(biāo)準(zhǔn)正態(tài)分布函數(shù), Φ(·)為一維標(biāo)準(zhǔn)正態(tài)分布函數(shù), Φ?1(·) 為 Φ(·)的反函數(shù)。進(jìn)而,基于式(14),可以得到額定風(fēng)速條件下Hs和Tp的聯(lián)合概率密度函數(shù),如圖4所示。

    圖4 額定風(fēng)速下有效波高和譜峰周期的聯(lián)合概率密度函數(shù)Fig. 4 Joint probability density function of Hs and Tp under the rated wind speed condition

    表2 波浪參數(shù)概率分布建模結(jié)果Table 2 Probabilistic modeling results of wave parameters

    應(yīng)指出,本文采用Copula方法建立了二維隨機(jī)環(huán)境變量的條件聯(lián)合概率分布模型。對于高維隨機(jī)環(huán)境變量,可采用C-vine Copula方法便捷地建立其聯(lián)合概率分布模型[18-19]。

    3 浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠性分析

    3.1 概率密度演化理論簡介

    本文采用概率密度演化理論進(jìn)行額定風(fēng)速下發(fā)生極端風(fēng)剪時浮式風(fēng)機(jī)的可靠度分析。為此,先簡要介紹概率密度演化理論。

    不失一般性,對于多自由度動力系統(tǒng),考慮系統(tǒng)的隨機(jī)因素后,動力方程可寫為:

    式中:Y(t)為 系統(tǒng)的狀態(tài)響應(yīng)量;Y0為初始時刻的系統(tǒng)狀態(tài); Θ=(Θ1,...,Θs)為s維隨機(jī)向量,刻畫了系統(tǒng)中結(jié)構(gòu)參數(shù)和外部激勵的隨機(jī)性。 Θ的聯(lián)合概率密度函數(shù)記為pΘ(θ)。

    式(16)的解可寫為Y(t)=H(Θ,t),該系統(tǒng)中其他物理量的響應(yīng)Z(t)=(Z1(t),...,Zm(t))均 可通過Y(t)和(t)獲得?;诟怕适睾阍?,李杰和陳建兵[5]建立了廣義概率密度演化方程。當(dāng)只關(guān)注一個物理響應(yīng)量Z(t)時,廣義概率密度演化方程具有如下形式[5]:

    且具有如下的初值條件:

    式 中:pZΘ(z,θ,t)為 (Θ,Z(t))的 聯(lián) 合 概 率 密 度 函 數(shù);δ(·)為 Dirac函數(shù);z0為Z(t)的初值。

    因此,Z(t)的概率密度函數(shù)為:

    式中, ?Θ為 隨機(jī)向量 Θ的分布域。

    為了求解廣義概率密度演化方程,已發(fā)展了點(diǎn)演化和群演化兩類數(shù)值求解方法。本文采用點(diǎn)演化方法,具體的求解步驟可參考文獻(xiàn)[5]。

    基于概率密度演化理論發(fā)展了兩種結(jié)構(gòu)首次超越破壞可靠度分析方法,即,基于吸收邊界條件的可靠度分析方法和基于極值分布的可靠度分析方法[5]。本文采用基于極值分布的可靠度分析方法。首先構(gòu)造隨機(jī)過程Z(Θ,t)的等價極值事件:

    式中,T為結(jié)構(gòu)響應(yīng)的分析時間。

    進(jìn)而構(gòu)造如下的虛擬隨機(jī)過程[5]:

    由于Z(Θ,t) 的 極值分布等價于t=1時 刻Zvir(Θ,t)的概率分布,因此可先通過求解廣義概率密度演化方程獲得Zvir(Θ,t)的 概率密度演化過程,然后獲得Z(Θ,t)的極值分布,即:

    由此,可得到結(jié)構(gòu)的可靠度為:

    式中,b為Z(Θ,t)的失效閾值。

    3.2 算例

    針對風(fēng)力機(jī)結(jié)構(gòu)設(shè)計規(guī)范[31]中定義的額定風(fēng)速條件下發(fā)生極端風(fēng)剪的工況,進(jìn)行了單柱式浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析。此處僅考慮豎直方向的正極端風(fēng)剪,此時風(fēng)剖面為[31]:

    式中:Vhub為輪轂高度處的平均風(fēng)速;z和zhub分別為高度和輪轂高度;α=0.2為風(fēng)剖面指數(shù)律的冪指數(shù);T0=12s為 極端風(fēng)剪的持續(xù)時間;t為 時間;D為風(fēng)輪直徑; β=6.4 ; σ1為 脈 動 風(fēng) 速 標(biāo) 準(zhǔn) 差 的 代 表 值,σ1=Iref(0.75Vhub+5.6), 其中Iref為 湍 流 強(qiáng) 度; Λ1為湍流積分尺度; σ1和 Λ1的取值可根據(jù)規(guī)范[31]確定。當(dāng)Vhub=15m/s、Iref=0.12 時,5 MW風(fēng)力機(jī)在t=T0/2時的極端風(fēng)剖面和正常風(fēng)剖面對比圖見圖5。

    圖5 正常風(fēng)剖面和極端風(fēng)剖面對比Fig. 5 Comparison of normal and extreme wind profiles

    在分析結(jié)構(gòu)的可靠度時,應(yīng)當(dāng)充分、合理考慮系統(tǒng)中的隨機(jī)源。本文采用的概率密度演化理論在分析結(jié)構(gòu)可靠度時能夠同時全面考慮結(jié)構(gòu)自身和外部激勵的隨機(jī)性[5]。作為示例,本文不僅將額定風(fēng)速條件下的有效波高Hs和 譜峰周期Tp作為隨機(jī)變量,而且考慮了脈動風(fēng)速、波面高程以及塔筒阻尼比 ξT的隨機(jī)性。本文基于譜表達(dá)方法并分別采用波數(shù)-頻率聯(lián)合功率譜[34-35]和JONSWAP譜[13]模擬隨機(jī)脈動風(fēng)場和隨機(jī)波浪場。為了降低隨機(jī)場模擬中的隨機(jī)變量(相位)個數(shù),采用演化相位譜方法[36],由此分別引入4個和1個相位演化時間,分別記為Te1~Te4和Te5。相位演化時間為隨機(jī)變量,基于實(shí)測風(fēng)速數(shù)據(jù)和波浪數(shù)據(jù)可以識別其概率分布[36]。相位演化時間和塔筒阻尼比的概率分布信息如表3所示[9]。Hs和Tp的概率分布信息由第2節(jié)方法確定。

    表3 隨機(jī)變量的概率分布信息Table 3 Probability distributions of random variables

    在求解廣義概率密度演化方程時,首先需要將概率空間剖分為Nq個子域,并在每個子域內(nèi)選取一個代表點(diǎn)[5]。然后對每個代表點(diǎn)進(jìn)行確定性動力響應(yīng)分析以獲得概率密度演化方程(偏微分方程)中的時變系數(shù),進(jìn)而求解概率密度演化方程,并獲得響應(yīng)量的概率密度演化結(jié)果[5]。本文共考慮8個隨機(jī)變量,采用基于點(diǎn)集GF-偏差最小的點(diǎn)集優(yōu)選策略[37]在八維概率空間中選取了329個代表點(diǎn)。圖6中的紅點(diǎn)展示的是Hs和Tp的代表點(diǎn),圖中的灰線為Hs和Tp聯(lián)合概率密度函數(shù)的等值線。

    圖6 額定風(fēng)速下隨機(jī)動力分析中的波浪參數(shù)代表點(diǎn)Fig. 6 Representative points of wave parameters for the stochastic analysis under rated wind speed condition

    在選定代表點(diǎn)之后,需要模擬每種情況下的隨機(jī)風(fēng)場和波浪場,并分析浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)。作為示例,圖7給出了Hs=2m、Tp=8.4s工況下浮式風(fēng)機(jī)結(jié)構(gòu)動力響應(yīng)的典型樣本。同時分析了不考慮極端風(fēng)剪時的結(jié)構(gòu)響應(yīng)并進(jìn)行對比。動力響應(yīng)的時間分析長度為1500 s,并從第500 s開始施加極端風(fēng)剪條件。為了避免初始瞬態(tài)響應(yīng)的影響,僅展示后1200 s的動力激勵和結(jié)構(gòu)動力響應(yīng)結(jié)果。圖中,“EWS”表示極端風(fēng)剪。

    圖7 浮式風(fēng)機(jī)結(jié)構(gòu)在特定工況下的響應(yīng)Fig. 7 Responses of the floating offshore wind turbine under specified load conditions

    從圖7(a)可以看出,在施加極端風(fēng)剪條件后,z=30m 處和z=150m處的風(fēng)速差異明顯增大,由此導(dǎo)致塔筒和槳葉的位移響應(yīng)發(fā)生變化(如圖7d、圖7f所示)。從圖7(d)還可以看出,考慮極端風(fēng)剪后,該樣本中塔頂位移響應(yīng)的極值有所增大,因此會影響結(jié)構(gòu)的整體可靠度結(jié)果。圖7(c)和圖7(f)顯示,大約在700~750 s時間段內(nèi),結(jié)構(gòu)響應(yīng)有突然的“跌落”,這是由于啟動了變槳控制,降低了氣動荷載(如圖7(j)所示)。圖7(g)~圖7(j)表明,極端風(fēng)剪對浮式風(fēng)機(jī)的縱蕩響應(yīng)、縱搖響應(yīng)、槳葉轉(zhuǎn)速、槳葉的槳距角變化均影響較小。

    采用概率密度演化理論可以獲得結(jié)構(gòu)響應(yīng)極值的概率分布。以塔頂順風(fēng)向位移為例,通過求解廣義概率密度演化方程,得到其響應(yīng)極值的概率密度函數(shù)如圖8所示。在給定的失效閾值條件下,當(dāng)塔頂位移響應(yīng)超過該閾值時,即認(rèn)為發(fā)生失效。由此可以得到不同失效閾值條件下的失效概率,如圖9所示。

    圖8 塔頂順風(fēng)向位移響應(yīng)極值的概率密度函數(shù)Fig. 8 Probability density function of the extreme value of the fore-aft displacement response at the tower top

    圖9 失效概率隨失效閾值變化關(guān)系Fig. 9 Relationship between failure probability and failure threshold

    作為示例,本文在分析浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度時,僅考慮了塔頂位移。在實(shí)際中,應(yīng)考慮多個指標(biāo),例如塔頂加速度、塔底應(yīng)力以及浮式風(fēng)機(jī)的縱搖響應(yīng)等。概率密度演化理論能夠同時考慮多個失效指標(biāo)[5]。當(dāng)同時采用多個失效指標(biāo)分析時,可靠度結(jié)果可能發(fā)生變化。

    4 結(jié)論

    浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠性分析是進(jìn)行基于可靠度優(yōu)化設(shè)計的關(guān)鍵基礎(chǔ),對于大規(guī)模深遠(yuǎn)海風(fēng)能的開發(fā)具有重要意義。本文針對海上單柱式浮式風(fēng)機(jī)結(jié)構(gòu),首先介紹了一體化耦合動力學(xué)分析模型,然后建立了環(huán)境參數(shù)的聯(lián)合概率分布模型以合理考慮風(fēng)-浪聯(lián)合作用,進(jìn)而分析浮式風(fēng)機(jī)結(jié)構(gòu)在極端工況條件下的動力可靠度。主要研究結(jié)論如下:

    1)對于當(dāng)前兆瓦級浮式風(fēng)機(jī),柔性部件的彈性變形對浮式風(fēng)機(jī)的剛體運(yùn)動響應(yīng)幾乎沒有影響,因此在建模過程中可以忽略這一因素,從而使建模方法與分析過程更為簡便。本文基于這一思想,將多剛體動力學(xué)理論與有限元方法結(jié)合起來考察了浮式風(fēng)機(jī)結(jié)構(gòu)的一體化模型,不僅能夠彌補(bǔ)一般多剛體動力學(xué)模型的不足,而且利用Kane方法便捷地考慮了系統(tǒng)中剛體運(yùn)動對柔性部件變形的影響。雖然從理論上看,全耦合多柔性體模型(例如FAST)的建模方法考慮得更為全面,但是從實(shí)際工程應(yīng)用的角度來看,本文所提的方法更為便捷。兩種模型的對比分析結(jié)果表明了本文模型具有良好的精度。

    2)在浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)分析中,應(yīng)當(dāng)合理考慮風(fēng)-浪聯(lián)合作用?;趫鲋返拈L期風(fēng)-浪觀測數(shù)據(jù)并采用Copula方法,能夠精確、便捷地建立風(fēng)-浪多參數(shù)非線性相關(guān)的聯(lián)合概率分布模型,從而合理地考慮風(fēng)和波浪對結(jié)構(gòu)的荷載效應(yīng)。

    3)基于概率密度演化理論的可靠度分析方法適用于海上大型浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠性分析。該方法不僅能夠同時考慮外部激勵和結(jié)構(gòu)自身的隨機(jī)性,而且能夠獲得結(jié)構(gòu)響應(yīng)的概率密度演化過程,同時具有很高的分析效率。

    本文研究建立了海上大型浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠性分析的基本框架,為進(jìn)行基于整體可靠度的結(jié)構(gòu)優(yōu)化設(shè)計奠定了基礎(chǔ)。由于我國東南海域的風(fēng)電場常年遭受臺風(fēng)侵襲[38],因此在未來工作中需要進(jìn)一步分析在臺風(fēng)等極端環(huán)境條件下的浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠度,并開展基于整體可靠度的結(jié)構(gòu)優(yōu)化設(shè)計等研究。

    猜你喜歡
    結(jié)構(gòu)分析方法
    《形而上學(xué)》△卷的結(jié)構(gòu)和位置
    隱蔽失效適航要求符合性驗(yàn)證分析
    論結(jié)構(gòu)
    中華詩詞(2019年7期)2019-11-25 01:43:04
    電力系統(tǒng)不平衡分析
    電子制作(2018年18期)2018-11-14 01:48:24
    電力系統(tǒng)及其自動化發(fā)展趨勢分析
    可能是方法不對
    論《日出》的結(jié)構(gòu)
    用對方法才能瘦
    Coco薇(2016年2期)2016-03-22 02:42:52
    四大方法 教你不再“坐以待病”!
    Coco薇(2015年1期)2015-08-13 02:47:34
    捕魚
    精品久久久久久久久亚洲| www.色视频.com| 麻豆国产97在线/欧美| 欧美成人a在线观看| 男人狂女人下面高潮的视频| 中文字幕av在线有码专区| 亚洲第一区二区三区不卡| 26uuu在线亚洲综合色| 97在线视频观看| 性色avwww在线观看| 午夜福利成人在线免费观看| 欧美激情国产日韩精品一区| 久久99蜜桃精品久久| 国产亚洲午夜精品一区二区久久 | 色播亚洲综合网| 国产成人aa在线观看| 男女那种视频在线观看| 国内精品美女久久久久久| 精品一区二区免费观看| 午夜老司机福利剧场| 蜜桃久久精品国产亚洲av| 人妻制服诱惑在线中文字幕| 男女啪啪激烈高潮av片| 激情 狠狠 欧美| 搡老乐熟女国产| 成人性生交大片免费视频hd| 少妇的逼水好多| 欧美日韩一区二区视频在线观看视频在线 | 人体艺术视频欧美日本| 国产白丝娇喘喷水9色精品| 国产成人一区二区在线| 免费无遮挡裸体视频| 黄片无遮挡物在线观看| 亚洲av男天堂| 视频中文字幕在线观看| 一个人观看的视频www高清免费观看| 最后的刺客免费高清国语| 久久久久精品性色| 亚洲欧美日韩卡通动漫| 精品一区在线观看国产| 青春草亚洲视频在线观看| 中国国产av一级| 久久久精品94久久精品| 深爱激情五月婷婷| 91在线精品国自产拍蜜月| 欧美zozozo另类| 亚洲国产最新在线播放| 国产亚洲av嫩草精品影院| 国产淫片久久久久久久久| 久99久视频精品免费| 三级毛片av免费| 日韩欧美一区视频在线观看 | 高清欧美精品videossex| 中文字幕人妻熟人妻熟丝袜美| 久久久久久久久久人人人人人人| 国产精品国产三级国产av玫瑰| 午夜老司机福利剧场| 亚洲天堂国产精品一区在线| 成人毛片60女人毛片免费| 国产男人的电影天堂91| 国产成人aa在线观看| 中文字幕av在线有码专区| 亚洲av国产av综合av卡| 日日干狠狠操夜夜爽| 亚洲av电影不卡..在线观看| 精品国产露脸久久av麻豆 | 欧美极品一区二区三区四区| 国产男女超爽视频在线观看| 国产男人的电影天堂91| 国产精品国产三级国产专区5o| 免费av观看视频| 2021天堂中文幕一二区在线观| 18禁在线播放成人免费| 熟妇人妻不卡中文字幕| av.在线天堂| 国产精品人妻久久久影院| 国产在视频线精品| 午夜精品一区二区三区免费看| 老司机影院毛片| 国产 亚洲一区二区三区 | 日韩成人av中文字幕在线观看| 男插女下体视频免费在线播放| 好男人视频免费观看在线| 亚洲性久久影院| 久久亚洲国产成人精品v| 欧美日本视频| 内射极品少妇av片p| 赤兔流量卡办理| 一级爰片在线观看| 国产成人午夜福利电影在线观看| 大香蕉久久网| 91久久精品电影网| 国产 一区精品| 成人国产麻豆网| 国产亚洲午夜精品一区二区久久 | 又粗又硬又长又爽又黄的视频| 国产视频首页在线观看| 最近视频中文字幕2019在线8| av在线播放精品| 国产在视频线精品| av在线老鸭窝| 国产成人精品久久久久久| 亚洲一级一片aⅴ在线观看| 最近手机中文字幕大全| 久久久久久久久中文| 免费高清在线观看视频在线观看| 在线观看免费高清a一片| 亚洲av在线观看美女高潮| 九九在线视频观看精品| 免费不卡的大黄色大毛片视频在线观看 | 在线播放无遮挡| av播播在线观看一区| 最近中文字幕高清免费大全6| 中文字幕制服av| 老司机影院毛片| 亚洲av男天堂| 少妇熟女欧美另类| 欧美不卡视频在线免费观看| 国产精品国产三级国产av玫瑰| 久久久久网色| 国产成人一区二区在线| 国产三级在线视频| 午夜福利在线观看吧| 99热这里只有是精品50| 午夜福利在线在线| 亚洲婷婷狠狠爱综合网| 成人一区二区视频在线观看| 久久久午夜欧美精品| 69人妻影院| 特大巨黑吊av在线直播| 国产熟女欧美一区二区| 国产精品日韩av在线免费观看| 久久99热6这里只有精品| 看免费成人av毛片| 国产男人的电影天堂91| 久久久久性生活片| 黄色日韩在线| 激情五月婷婷亚洲| 熟妇人妻不卡中文字幕| 欧美精品一区二区大全| 女人十人毛片免费观看3o分钟| 欧美zozozo另类| 直男gayav资源| 国产成人精品久久久久久| 亚洲国产欧美在线一区| 久久午夜福利片| 一区二区三区乱码不卡18| 18禁裸乳无遮挡免费网站照片| 亚洲精品456在线播放app| 国产激情偷乱视频一区二区| 国产日韩欧美在线精品| 欧美激情国产日韩精品一区| 亚洲第一区二区三区不卡| 国产熟女欧美一区二区| 黄色日韩在线| 久久99蜜桃精品久久| 欧美一区二区亚洲| 免费黄频网站在线观看国产| 在线观看av片永久免费下载| 久久综合国产亚洲精品| 黄色欧美视频在线观看| 久久99精品国语久久久| 日本一二三区视频观看| 联通29元200g的流量卡| 久久这里只有精品中国| 精品一区二区三区视频在线| 色综合色国产| 国产男人的电影天堂91| 精品久久久久久成人av| 久久精品人妻少妇| 久久97久久精品| 最近最新中文字幕大全电影3| 22中文网久久字幕| 91在线精品国自产拍蜜月| 少妇人妻精品综合一区二区| 亚洲欧美成人综合另类久久久| 亚洲成人中文字幕在线播放| 成年女人看的毛片在线观看| 日本-黄色视频高清免费观看| 亚洲欧美精品专区久久| 亚洲18禁久久av| 少妇猛男粗大的猛烈进出视频 | 最后的刺客免费高清国语| 欧美 日韩 精品 国产| 97超碰精品成人国产| 日本免费在线观看一区| 婷婷色综合大香蕉| 国产男人的电影天堂91| av国产久精品久网站免费入址| 男人和女人高潮做爰伦理| 亚洲一区高清亚洲精品| 80岁老熟妇乱子伦牲交| 男女啪啪激烈高潮av片| 日韩强制内射视频| 成人无遮挡网站| 国产一区亚洲一区在线观看| 一级二级三级毛片免费看| 少妇熟女欧美另类| 欧美一级a爱片免费观看看| 大香蕉久久网| 大话2 男鬼变身卡| 色综合亚洲欧美另类图片| 18禁在线播放成人免费| 51国产日韩欧美| 国产日韩欧美在线精品| 夫妻性生交免费视频一级片| 亚洲精品,欧美精品| 国产欧美另类精品又又久久亚洲欧美| 舔av片在线| 亚洲欧美日韩卡通动漫| 人妻少妇偷人精品九色| 成人亚洲精品av一区二区| 亚洲欧美成人精品一区二区| 日本熟妇午夜| 日本黄大片高清| 99热网站在线观看| 午夜免费激情av| 观看美女的网站| 麻豆国产97在线/欧美| 亚洲欧美中文字幕日韩二区| 精品国内亚洲2022精品成人| 人妻夜夜爽99麻豆av| 有码 亚洲区| 国产人妻一区二区三区在| 欧美97在线视频| 免费观看av网站的网址| 日韩不卡一区二区三区视频在线| 深爱激情五月婷婷| 久久精品国产亚洲网站| 国内精品一区二区在线观看| 91久久精品国产一区二区三区| 久久久久久久大尺度免费视频| 免费观看av网站的网址| 在线观看免费高清a一片| 精品久久久久久久人妻蜜臀av| 亚洲怡红院男人天堂| 亚洲人与动物交配视频| 日本黄大片高清| 69av精品久久久久久| 久久这里有精品视频免费| 成人国产麻豆网| 久久精品国产亚洲网站| 国产爱豆传媒在线观看| 一本一本综合久久| 波野结衣二区三区在线| 久久6这里有精品| 秋霞伦理黄片| 中文精品一卡2卡3卡4更新| 如何舔出高潮| 亚洲自拍偷在线| 精品一区二区三区视频在线| 亚洲国产高清在线一区二区三| 免费观看性生交大片5| a级一级毛片免费在线观看| 亚洲伊人久久精品综合| 亚洲久久久久久中文字幕| 91午夜精品亚洲一区二区三区| 最后的刺客免费高清国语| 久久精品国产亚洲av涩爱| 日本猛色少妇xxxxx猛交久久| 精品99又大又爽又粗少妇毛片| 国产精品综合久久久久久久免费| 日本欧美国产在线视频| 国产男女超爽视频在线观看| 日韩av在线大香蕉| 天天一区二区日本电影三级| 欧美性感艳星| 国产成人午夜福利电影在线观看| 国产亚洲午夜精品一区二区久久 | 夜夜爽夜夜爽视频| 久久久久久久久大av| 亚洲一区高清亚洲精品| 亚洲真实伦在线观看| 午夜免费男女啪啪视频观看| 最近最新中文字幕免费大全7| 亚洲精品456在线播放app| 日韩人妻高清精品专区| 最近手机中文字幕大全| 最近中文字幕2019免费版| 日本av手机在线免费观看| 国产伦精品一区二区三区四那| 久久久久久伊人网av| 超碰97精品在线观看| 国产女主播在线喷水免费视频网站 | 中文字幕久久专区| 久久久久久久久久黄片| 精品久久久久久久久久久久久| 亚洲乱码一区二区免费版| 成人性生交大片免费视频hd| 两个人的视频大全免费| 好男人在线观看高清免费视频| h日本视频在线播放| 淫秽高清视频在线观看| 国产探花在线观看一区二区| 国产男人的电影天堂91| 欧美一级a爱片免费观看看| 菩萨蛮人人尽说江南好唐韦庄| 亚洲精品456在线播放app| 亚洲精品日韩在线中文字幕| 国产黄色视频一区二区在线观看| 中文天堂在线官网| 成人av在线播放网站| 精品国产露脸久久av麻豆 | 中文资源天堂在线| 国产视频首页在线观看| 精品国内亚洲2022精品成人| av又黄又爽大尺度在线免费看| 亚洲图色成人| 99久久精品热视频| 亚洲不卡免费看| 深夜a级毛片| 国产乱来视频区| 日韩人妻高清精品专区| 欧美xxxx性猛交bbbb| 最近中文字幕2019免费版| 国国产精品蜜臀av免费| 欧美3d第一页| 少妇熟女aⅴ在线视频| 久久久久精品久久久久真实原创| 麻豆成人午夜福利视频| 亚洲第一区二区三区不卡| 好男人视频免费观看在线| 国产成人精品福利久久| 一二三四中文在线观看免费高清| 国产精品.久久久| 免费观看性生交大片5| 日日摸夜夜添夜夜添av毛片| 欧美一区二区亚洲| 搞女人的毛片| av在线天堂中文字幕| 午夜免费男女啪啪视频观看| 男女边吃奶边做爰视频| 成人一区二区视频在线观看| 一级片'在线观看视频| 国产一区有黄有色的免费视频 | 国产一级毛片在线| 国产亚洲5aaaaa淫片| 免费黄色在线免费观看| 蜜臀久久99精品久久宅男| 成人av在线播放网站| 亚洲av日韩在线播放| 毛片女人毛片| 国产午夜福利久久久久久| 精品久久久精品久久久| 九九在线视频观看精品| 插阴视频在线观看视频| 全区人妻精品视频| 男插女下体视频免费在线播放| 国产老妇女一区| 亚洲高清免费不卡视频| 插逼视频在线观看| 2021少妇久久久久久久久久久| 一级毛片黄色毛片免费观看视频| 亚洲最大成人手机在线| 夫妻性生交免费视频一级片| 色综合站精品国产| 免费观看a级毛片全部| 女人十人毛片免费观看3o分钟| 18+在线观看网站| 国产精品一区二区三区四区久久| 熟妇人妻不卡中文字幕| 亚洲精品国产av蜜桃| 秋霞伦理黄片| 一级毛片 在线播放| 伦理电影大哥的女人| 免费观看性生交大片5| av在线亚洲专区| 搡老乐熟女国产| 国产精品不卡视频一区二区| 一个人免费在线观看电影| av在线亚洲专区| 一级毛片久久久久久久久女| 黄色配什么色好看| 国产成人a∨麻豆精品| 一个人看的www免费观看视频| 国产精品久久久久久av不卡| 一级毛片黄色毛片免费观看视频| 在现免费观看毛片| ponron亚洲| 欧美xxxx性猛交bbbb| 亚洲av中文av极速乱| 韩国av在线不卡| 只有这里有精品99| 99久国产av精品| 精品久久久久久久久亚洲| 欧美极品一区二区三区四区| 美女国产视频在线观看| 男女那种视频在线观看| 麻豆国产97在线/欧美| 天堂av国产一区二区熟女人妻| 精品久久国产蜜桃| 中国美白少妇内射xxxbb| 免费看美女性在线毛片视频| 天天一区二区日本电影三级| 18禁在线无遮挡免费观看视频| 国产高清不卡午夜福利| 午夜精品在线福利| 日韩人妻高清精品专区| 啦啦啦啦在线视频资源| 男人舔女人下体高潮全视频| 亚洲精品一区蜜桃| 亚洲av免费在线观看| 久久久久久久亚洲中文字幕| 日韩av免费高清视频| 欧美97在线视频| 男女边吃奶边做爰视频| 亚洲,欧美,日韩| av在线播放精品| 一个人看的www免费观看视频| 97热精品久久久久久| 免费av不卡在线播放| 日本一二三区视频观看| 久久热精品热| 又爽又黄无遮挡网站| 26uuu在线亚洲综合色| 日本三级黄在线观看| 免费电影在线观看免费观看| 成人欧美大片| 国产综合精华液| 不卡视频在线观看欧美| videossex国产| 欧美xxⅹ黑人| 国产中年淑女户外野战色| 国产av在哪里看| 女人被狂操c到高潮| 3wmmmm亚洲av在线观看| 国产亚洲av片在线观看秒播厂 | 欧美成人一区二区免费高清观看| 国产精品久久久久久久电影| 午夜福利视频1000在线观看| 日韩人妻高清精品专区| 一级av片app| 草草在线视频免费看| 丝袜喷水一区| 国产视频首页在线观看| 成人高潮视频无遮挡免费网站| 午夜福利视频精品| 国产午夜福利久久久久久| 亚洲成人久久爱视频| 免费看a级黄色片| 99久国产av精品| 色综合亚洲欧美另类图片| 午夜精品一区二区三区免费看| 97在线视频观看| 国产成人精品一,二区| 日韩大片免费观看网站| 欧美变态另类bdsm刘玥| 成人午夜高清在线视频| 国产精品久久久久久久电影| 伊人久久精品亚洲午夜| 国精品久久久久久国模美| 在线a可以看的网站| 大香蕉97超碰在线| 国产av在哪里看| 午夜精品国产一区二区电影 | 日本猛色少妇xxxxx猛交久久| 亚洲欧美成人精品一区二区| 国产成人a∨麻豆精品| 午夜福利在线观看免费完整高清在| 国产精品美女特级片免费视频播放器| 亚洲熟女精品中文字幕| 国产乱人偷精品视频| xxx大片免费视频| 国产又色又爽无遮挡免| 少妇猛男粗大的猛烈进出视频 | 一级片'在线观看视频| 欧美日本视频| 国产精品久久视频播放| 国产精品综合久久久久久久免费| 蜜桃亚洲精品一区二区三区| 亚洲一级一片aⅴ在线观看| 精品久久久久久久人妻蜜臀av| 搡老妇女老女人老熟妇| 国产熟女欧美一区二区| 免费在线观看成人毛片| 男女下面进入的视频免费午夜| 中文在线观看免费www的网站| 精品久久国产蜜桃| 久久久久国产网址| 午夜福利视频1000在线观看| av在线观看视频网站免费| 精品久久久久久成人av| 久久人人爽人人爽人人片va| 亚洲性久久影院| av国产久精品久网站免费入址| 一区二区三区高清视频在线| 蜜桃亚洲精品一区二区三区| 爱豆传媒免费全集在线观看| 丝袜喷水一区| 97超视频在线观看视频| 又爽又黄a免费视频| 国产又色又爽无遮挡免| 亚洲国产精品sss在线观看| 97精品久久久久久久久久精品| 我的老师免费观看完整版| 18禁裸乳无遮挡免费网站照片| 欧美激情久久久久久爽电影| 国语对白做爰xxxⅹ性视频网站| 一区二区三区高清视频在线| 18禁动态无遮挡网站| 国产一区亚洲一区在线观看| 久久久久久久久久成人| 春色校园在线视频观看| 亚洲乱码一区二区免费版| 尤物成人国产欧美一区二区三区| av一本久久久久| 国内精品一区二区在线观看| 成年人午夜在线观看视频 | 一二三四中文在线观看免费高清| av在线蜜桃| 欧美3d第一页| 亚州av有码| 91精品一卡2卡3卡4卡| 午夜福利视频精品| 在线观看av片永久免费下载| 亚洲久久久久久中文字幕| 日韩av免费高清视频| 成人高潮视频无遮挡免费网站| 亚洲av免费高清在线观看| 亚洲自拍偷在线| av在线老鸭窝| 一个人免费在线观看电影| 免费观看在线日韩| 国产亚洲一区二区精品| 亚洲一级一片aⅴ在线观看| 欧美日本视频| 男插女下体视频免费在线播放| 国产亚洲5aaaaa淫片| 中文字幕人妻熟人妻熟丝袜美| or卡值多少钱| 国产av码专区亚洲av| 建设人人有责人人尽责人人享有的 | 直男gayav资源| a级毛片免费高清观看在线播放| 国产黄片美女视频| 在线免费观看不下载黄p国产| 国产伦精品一区二区三区视频9| 成人美女网站在线观看视频| 国产精品一区二区三区四区久久| 久久精品熟女亚洲av麻豆精品 | 女人被狂操c到高潮| 亚洲自拍偷在线| 国产91av在线免费观看| 亚洲国产欧美人成| 黄色日韩在线| 麻豆精品久久久久久蜜桃| 亚洲怡红院男人天堂| 18+在线观看网站| 最后的刺客免费高清国语| 2021天堂中文幕一二区在线观| 国产 一区 欧美 日韩| 舔av片在线| 一级毛片黄色毛片免费观看视频| 日韩 亚洲 欧美在线| 婷婷色麻豆天堂久久| 日韩大片免费观看网站| 欧美xxxx性猛交bbbb| 午夜日本视频在线| 免费看av在线观看网站| 热99在线观看视频| 又爽又黄a免费视频| 国产精品av视频在线免费观看| 亚洲精品乱码久久久v下载方式| av播播在线观看一区| 国产精品蜜桃在线观看| 婷婷六月久久综合丁香| 国产大屁股一区二区在线视频| 国产精品一区二区三区四区久久| 国产亚洲av片在线观看秒播厂 | 插阴视频在线观看视频| 国产精品无大码| 大又大粗又爽又黄少妇毛片口| 啦啦啦中文免费视频观看日本| 国产精品综合久久久久久久免费| 99久久人妻综合| 欧美性感艳星| 我的老师免费观看完整版| 内地一区二区视频在线| 噜噜噜噜噜久久久久久91| 高清午夜精品一区二区三区| 免费看日本二区| 国产亚洲午夜精品一区二区久久 | av免费在线看不卡| 亚洲性久久影院| 精品一区二区三区视频在线| 国产伦精品一区二区三区视频9| 日韩成人伦理影院| 国产精品熟女久久久久浪| 少妇人妻一区二区三区视频| 中文字幕久久专区| 美女内射精品一级片tv| 国产亚洲5aaaaa淫片| 午夜精品一区二区三区免费看| 少妇猛男粗大的猛烈进出视频 | 少妇裸体淫交视频免费看高清| 两个人的视频大全免费| av线在线观看网站| 欧美一区二区亚洲| 成人毛片a级毛片在线播放| 久久精品国产亚洲av天美| 深爱激情五月婷婷| 国产高潮美女av| av黄色大香蕉| 91在线精品国自产拍蜜月| 男女边摸边吃奶| 免费黄网站久久成人精品| 久久久精品94久久精品| 欧美最新免费一区二区三区| 麻豆精品久久久久久蜜桃| 又粗又硬又长又爽又黄的视频| 夫妻性生交免费视频一级片| 亚洲精品中文字幕在线视频 | 国产 一区 欧美 日韩| 久久6这里有精品| 国产亚洲5aaaaa淫片|