• <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陳建兵宋玉鵬
    關(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è)計(jì)方法,實(shí)現(xiàn)兼顧浮式風(fēng)機(jī)結(jié)構(gòu)安全性與經(jīng)濟(jì)性的定量設(shè)計(jì),為深遠(yuǎn)海風(fēng)能的大規(guī)模開發(fā)提供關(guān)鍵技術(shù)保障。因此,浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析具有重要意義。

    當(dāng)前,復(fù)雜結(jié)構(gòu)的整體可靠度評估通常采用數(shù)值分析方法,主要涉及三個(gè)方面的問題:結(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)中同時(shí)存在多種復(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è)計(jì)方法,難以合理考慮上述復(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è)計(jì)和初始設(shè)計(jì)階段[11]。多柔性體模型能夠同時(shí)考慮系統(tǒng)中的剛體運(yùn)動和彈性變形響應(yīng),其中一般采用有限元方法或者模態(tài)疊加方法建立柔性部件的彈性模型[7]。雖然多柔性體模型能夠全面反映系統(tǒng)的力學(xué)特性和復(fù)雜的耦合效應(yīng),但是該類模型建模較為復(fù)雜,同時(shí)存在分析效率不高的問題,難以應(yīng)用于需要進(jìn)行大量數(shù)值模擬分析的問題中,例如結(jié)構(gòu)的可靠度評估與優(yōu)化設(shè)計(jì)等[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)速時(shí)程和20 min至3~6 h內(nèi)的波浪時(shí)程均可以認(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è)計(jì)規(guī)范中所推薦的方法[13]。在建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型時(shí),可采用條件概率分布建模方法[15]、Nataf變換方法[16]和Copula方法[17-18]。條件概率分布建模方法是當(dāng)前海洋結(jié)構(gòu)設(shè)計(jì)規(guī)范中推薦使用的方法[13]。該方法首先確定主導(dǎo)隨機(jī)變量,并建立其邊緣概率分布模型,然后根據(jù)主導(dǎo)隨機(jī)變量的取值將數(shù)據(jù)劃分為一系列區(qū)間,進(jìn)而估計(jì)主導(dǎo)隨機(jī)變量取不同值條件下的其他隨機(jī)變量的概率分布。但是,當(dāng)隨機(jī)變量的個(gè)數(shù)較多時(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)時(shí),例如風(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)。該方法目前已在多個(gè)領(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)分析的代理模型以提高效率,但同時(shí)也引入了新的誤差[21]。概率密度演化理論給出了隨機(jī)動力系統(tǒng)響應(yīng)的概率密度函數(shù)的演化方程,揭示了系統(tǒng)響應(yīng)隨機(jī)性演化的內(nèi)在物理規(guī)律[5],該方法在分析系統(tǒng)的可靠度時(shí)同樣不受問題類型、問題維度和非線性程度等因素的影響,能給出系統(tǒng)響應(yīng)量的概率密度函數(shù)演化過程,同時(shí)具有較高的分析效率,因此適用于浮式風(fēng)機(jī)的可靠度分析,但目前在此領(lǐng)域的應(yīng)用研究較少。

    為了分析浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度,首先引入結(jié)合多剛體動力學(xué)理論和有限元方法建立起來的浮式風(fēng)機(jī)結(jié)構(gòu)一體化分析模型(StoDRAFOWT模型)。該模型不僅能同時(shí)獲得浮式風(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è)計(jì)中需要考慮的額定風(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軸豎直向上。在初始時(shí)刻浮式風(fēng)機(jī)未發(fā)生位移時(shí),坐標(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],即包含七個(gè)自由度。由于現(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)角,同時(shí)還存在槳距控制。因此,槳葉在風(fēng)輪平面內(nèi)、外兩個(gè)方向的振動相互耦合,在建模過程中應(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 荷載計(jì)算

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

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

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

    在計(jì)算浮式風(fēng)機(jī)的錨泊荷載時(shí),需要建立錨索的力學(xué)分析模型,主要包括錨索的靜力模型、擬靜力模型和動力模型[28]。由于靜力和擬靜力模型在某些工況下會低估錨索的疲勞荷載,因此推薦采用錨索的動力模型[28]。本文采用錨索的集中質(zhì)量動力模型[28]計(jì)算錨泊荷載。該模型不僅考慮錨索的受拉彈性變形,而且將錨索的動力效應(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ī)在服役期間同時(shí)受到風(fēng)荷載和波浪荷載作用,因此在進(jìn)行浮式風(fēng)機(jī)的動力響應(yīng)與可靠度分析時(shí),應(yīng)當(dāng)合理考慮風(fēng)-浪聯(lián)合作用。如前所述,當(dāng)前海洋結(jié)構(gòu)的設(shè)計(jì)規(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)剪時(shí)浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度。根據(jù)風(fēng)力機(jī)結(jié)構(gòu)的設(shè)計(jì)規(guī)范,在分析該工況下風(fēng)力機(jī)結(jié)構(gòu)的動力響應(yīng)時(shí),取風(fēng)-浪夾角為0°[31]。因此,需要建立額定風(fēng)速條件下Hs和Tp的條件聯(lián)合概率分布模型。

    在建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型時(shí),應(yīng)當(dāng)基于場址的長期風(fēng)-浪觀測數(shù)據(jù)。本文采用歐洲中尺度天氣預(yù)報(bào)中心ERA-Interim數(shù)據(jù)庫[32]提供的風(fēng)-浪參數(shù)再分析數(shù)據(jù),獲取了我國南海場址(21°N,113°E)處1979年至2018年(共計(jì)40年)的風(fēng)-浪再分析數(shù)據(jù)。這些數(shù)據(jù)均為30 min內(nèi)觀測數(shù)據(jù)的平均值,時(shí)間分辨率為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個(gè)數(shù)據(jù)點(diǎn)。

    根據(jù)Copula聯(lián)合概率分布建模的基本理論[17],對于任意兩個(gè)隨機(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)合概率分布模型時(shí),需要分別建立各隨機(jī)變量的邊緣概率分布模型及變量間的相關(guān)結(jié)構(gòu)模型。由于邊緣概率分布模型及變量間的相關(guān)結(jié)構(gòu)模型可以分開處理,因此該方法具有靈活、便捷、精確等優(yōu)點(diǎn)。

    基于上述681個(gè)額定風(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)剪時(shí)浮式風(fēng)機(jī)的可靠度分析。為此,先簡要介紹概率密度演化理論。

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

    式中:Y(t)為 系統(tǒng)的狀態(tài)響應(yīng)量;Y0為初始時(shí)刻的系統(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)注一個(gè)物理響應(yīng)量Z(t)時(shí),廣義概率密度演化方程具有如下形式[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)的等價(jià)極值事件:

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

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

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

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

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

    3.2 算例

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

    式中:Vhub為輪轂高度處的平均風(fēng)速;z和zhub分別為高度和輪轂高度;α=0.2為風(fēng)剖面指數(shù)律的冪指數(shù);T0=12s為 極端風(fēng)剪的持續(xù)時(shí)間;t為 時(shí)間;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 時(shí),5 MW風(fēng)力機(jī)在t=T0/2時(shí)的極端風(fēng)剖面和正常風(fēng)剖面對比圖見圖5。

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

    在分析結(jié)構(gòu)的可靠度時(shí),應(yīng)當(dāng)充分、合理考慮系統(tǒng)中的隨機(jī)源。本文采用的概率密度演化理論在分析結(jié)構(gòu)可靠度時(shí)能夠同時(shí)全面考慮結(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ī)變量(相位)個(gè)數(shù),采用演化相位譜方法[36],由此分別引入4個(gè)和1個(gè)相位演化時(shí)間,分別記為Te1~Te4和Te5。相位演化時(shí)間為隨機(jī)變量,基于實(shí)測風(fēng)速數(shù)據(jù)和波浪數(shù)據(jù)可以識別其概率分布[36]。相位演化時(shí)間和塔筒阻尼比的概率分布信息如表3所示[9]。Hs和Tp的概率分布信息由第2節(jié)方法確定。

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

    在求解廣義概率密度演化方程時(shí),首先需要將概率空間剖分為Nq個(gè)子域,并在每個(gè)子域內(nèi)選取一個(gè)代表點(diǎn)[5]。然后對每個(gè)代表點(diǎn)進(jìn)行確定性動力響應(yīng)分析以獲得概率密度演化方程(偏微分方程)中的時(shí)變系數(shù),進(jìn)而求解概率密度演化方程,并獲得響應(yīng)量的概率密度演化結(jié)果[5]。本文共考慮8個(gè)隨機(jī)變量,采用基于點(diǎn)集GF-偏差最小的點(diǎn)集優(yōu)選策略[37]在八維概率空間中選取了329個(gè)代表點(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)的典型樣本。同時(shí)分析了不考慮極端風(fēng)剪時(shí)的結(jié)構(gòu)響應(yīng)并進(jìn)行對比。動力響應(yīng)的時(shí)間分析長度為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時(shí)間段內(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)超過該閾值時(shí),即認(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í),僅考慮了塔頂位移。在實(shí)際中,應(yīng)考慮多個(gè)指標(biāo),例如塔頂加速度、塔底應(yīng)力以及浮式風(fēng)機(jī)的縱搖響應(yīng)等。概率密度演化理論能夠同時(shí)考慮多個(gè)失效指標(biāo)[5]。當(dāng)同時(shí)采用多個(gè)失效指標(biāo)分析時(shí),可靠度結(jié)果可能發(fā)生變化。

    4 結(jié)論

    浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠性分析是進(jìn)行基于可靠度優(yōu)化設(shè)計(jì)的關(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)合作用。基于場址的長期風(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)的整體可靠性分析。該方法不僅能夠同時(shí)考慮外部激勵和結(jié)構(gòu)自身的隨機(jī)性,而且能夠獲得結(jié)構(gòu)響應(yīng)的概率密度演化過程,同時(shí)具有很高的分析效率。

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

    猜你喜歡
    結(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国产在线视频色| 国产成人精品无人区| 国产欧美日韩一区二区三| 纯流量卡能插随身wifi吗| 亚洲中文日韩欧美视频| 欧美乱妇无乱码| 亚洲av第一区精品v没综合| 欧美激情 高清一区二区三区| 久久欧美精品欧美久久欧美| 亚洲国产精品合色在线| 午夜免费鲁丝| 日本在线视频免费播放| 十分钟在线观看高清视频www| 这个男人来自地球电影免费观看| 久久精品人人爽人人爽视色| 久久国产乱子伦精品免费另类| 亚洲三区欧美一区| 无限看片的www在线观看| 久久久久精品国产欧美久久久| 无人区码免费观看不卡| 国产精品久久视频播放| 亚洲中文av在线| 天堂动漫精品| 午夜福利一区二区在线看| 欧美亚洲日本最大视频资源| 久久婷婷成人综合色麻豆| 国产精品98久久久久久宅男小说| 可以在线观看毛片的网站| 国产av一区在线观看免费| 老司机午夜福利在线观看视频| 欧美人与性动交α欧美精品济南到| 免费人成视频x8x8入口观看| 日本免费a在线| 国产精品 国内视频| 国产亚洲av高清不卡| 天天一区二区日本电影三级 | avwww免费| 村上凉子中文字幕在线| 动漫黄色视频在线观看| 亚洲专区国产一区二区| 99国产精品一区二区三区| 色综合婷婷激情| 亚洲情色 制服丝袜| 免费久久久久久久精品成人欧美视频| 色av中文字幕| 很黄的视频免费| 99在线视频只有这里精品首页| 免费看十八禁软件| av超薄肉色丝袜交足视频| 俄罗斯特黄特色一大片| 大陆偷拍与自拍| 久久久水蜜桃国产精品网| 岛国视频午夜一区免费看| 欧美老熟妇乱子伦牲交| 午夜福利欧美成人| 亚洲自拍偷在线| 好男人电影高清在线观看| 老熟妇乱子伦视频在线观看| av电影中文网址| 亚洲av片天天在线观看| 午夜成年电影在线免费观看| 女人精品久久久久毛片| 97超级碰碰碰精品色视频在线观看| 99在线视频只有这里精品首页| 欧美成人性av电影在线观看| 国产一区二区三区综合在线观看| 99re在线观看精品视频| 999久久久国产精品视频| 亚洲人成网站在线播放欧美日韩| 一级黄色大片毛片| 97人妻天天添夜夜摸| 久久精品亚洲精品国产色婷小说| 久久人人爽av亚洲精品天堂| 在线播放国产精品三级| 97超级碰碰碰精品色视频在线观看| 99久久精品国产亚洲精品| 两人在一起打扑克的视频| 久久久久久亚洲精品国产蜜桃av| 久久精品国产99精品国产亚洲性色 | 999久久久国产精品视频| 一区二区三区高清视频在线| 国产视频一区二区在线看| 日韩精品中文字幕看吧| 久久久久国产精品人妻aⅴ院| 好看av亚洲va欧美ⅴa在| 日日爽夜夜爽网站| 91精品国产国语对白视频| 男男h啪啪无遮挡| 久久久久久久久中文| 看黄色毛片网站| 99久久99久久久精品蜜桃| netflix在线观看网站| 18禁黄网站禁片午夜丰满| 99精品久久久久人妻精品| 51午夜福利影视在线观看| 欧美激情久久久久久爽电影 | 黑人操中国人逼视频| 黄色丝袜av网址大全| www.熟女人妻精品国产| 亚洲一区高清亚洲精品| 一区在线观看完整版| 97超级碰碰碰精品色视频在线观看| 三级毛片av免费| 人人妻人人澡人人看| aaaaa片日本免费| 亚洲av电影不卡..在线观看| 色尼玛亚洲综合影院| 黄色a级毛片大全视频| 欧美久久黑人一区二区| 99国产精品99久久久久| av电影中文网址| 国产精品免费视频内射| 欧美成人午夜精品| 妹子高潮喷水视频| 亚洲专区中文字幕在线| 神马国产精品三级电影在线观看 | 亚洲成人国产一区在线观看| av欧美777| 男人舔女人下体高潮全视频| 久久国产精品男人的天堂亚洲| 亚洲精品在线观看二区| 村上凉子中文字幕在线| 成人18禁在线播放| 一级作爱视频免费观看| 中文字幕人妻熟女乱码| 首页视频小说图片口味搜索| 熟女少妇亚洲综合色aaa.| 最近最新中文字幕大全电影3 | 一夜夜www| 午夜亚洲福利在线播放| 大型av网站在线播放| 亚洲aⅴ乱码一区二区在线播放 | 一级毛片女人18水好多| 亚洲在线自拍视频| 成人欧美大片| 久久国产亚洲av麻豆专区| 看黄色毛片网站| 成人精品一区二区免费| 咕卡用的链子| 午夜福利影视在线免费观看| 91九色精品人成在线观看| 咕卡用的链子| 免费在线观看视频国产中文字幕亚洲| 欧美乱码精品一区二区三区| 两个人看的免费小视频| 午夜免费鲁丝| 亚洲午夜理论影院| 成人永久免费在线观看视频| 757午夜福利合集在线观看| 亚洲色图av天堂| 国产精品国产高清国产av| 在线观看66精品国产| 国产乱人伦免费视频| 每晚都被弄得嗷嗷叫到高潮| 18禁观看日本| 国产真人三级小视频在线观看| 免费av毛片视频| 女性生殖器流出的白浆| 精品免费久久久久久久清纯| 亚洲五月天丁香| 午夜福利视频1000在线观看 | 欧美 亚洲 国产 日韩一| 国产亚洲精品一区二区www| 国产成人免费无遮挡视频| 国产亚洲精品综合一区在线观看 | 亚洲少妇的诱惑av| 99国产极品粉嫩在线观看| 亚洲自拍偷在线| 他把我摸到了高潮在线观看| 岛国视频午夜一区免费看| 制服丝袜大香蕉在线| a在线观看视频网站| tocl精华| 国产精品久久电影中文字幕| 级片在线观看| 国产成人av激情在线播放| 国产1区2区3区精品| 日韩精品中文字幕看吧| 在线观看66精品国产| 亚洲av成人不卡在线观看播放网| 国产野战对白在线观看| 欧美午夜高清在线| 成人亚洲精品一区在线观看| 人成视频在线观看免费观看| 可以在线观看的亚洲视频| 亚洲第一电影网av| 成人欧美大片| 在线十欧美十亚洲十日本专区| 国产亚洲精品久久久久久毛片| 欧美老熟妇乱子伦牲交| 亚洲色图av天堂| 老司机靠b影院| 在线视频色国产色| 婷婷六月久久综合丁香| 亚洲熟妇熟女久久| 性色av乱码一区二区三区2| 69精品国产乱码久久久| 最新在线观看一区二区三区| 久久伊人香网站| 黄色毛片三级朝国网站| 两性午夜刺激爽爽歪歪视频在线观看 | 欧美大码av| 精品国产乱子伦一区二区三区| 久久精品亚洲熟妇少妇任你| 久久久久久免费高清国产稀缺| 神马国产精品三级电影在线观看 | 老鸭窝网址在线观看| 精品欧美国产一区二区三| 最新在线观看一区二区三区| 欧美中文综合在线视频| 一级毛片女人18水好多| 超碰成人久久| 一边摸一边做爽爽视频免费| av天堂在线播放| 深夜精品福利| 欧美激情高清一区二区三区| 电影成人av| 久久久精品国产亚洲av高清涩受| 国产麻豆69| 成人亚洲精品av一区二区| 国产精品一区二区免费欧美| 亚洲国产精品合色在线| 乱人伦中国视频| 午夜福利一区二区在线看| 欧美日韩亚洲综合一区二区三区_| 日韩大码丰满熟妇| 欧美中文日本在线观看视频| 黑人欧美特级aaaaaa片| 韩国av一区二区三区四区| 久久国产精品影院| 在线视频色国产色| 久久香蕉国产精品| 老鸭窝网址在线观看| 国产国语露脸激情在线看| 多毛熟女@视频| 黄色毛片三级朝国网站| 午夜两性在线视频| 99久久久亚洲精品蜜臀av| 久久狼人影院| 美女高潮到喷水免费观看| 国产精品久久视频播放| 久久久国产成人精品二区| 亚洲 欧美一区二区三区| 午夜激情av网站| 免费看a级黄色片| 日本精品一区二区三区蜜桃| 满18在线观看网站| 国产国语露脸激情在线看| 黑人欧美特级aaaaaa片| 久久精品亚洲精品国产色婷小说| 亚洲人成网站在线播放欧美日韩| 欧美在线黄色| 国产精品精品国产色婷婷| 国产熟女午夜一区二区三区| 日韩一卡2卡3卡4卡2021年| 一进一出抽搐动态| 久99久视频精品免费| 国产色视频综合| 又大又爽又粗| 亚洲一区高清亚洲精品| 亚洲精品国产色婷婷电影| 又黄又粗又硬又大视频| 国产亚洲精品久久久久久毛片| 国产精品久久久久久亚洲av鲁大| 日韩欧美三级三区| 一个人观看的视频www高清免费观看 | 一边摸一边做爽爽视频免费| 人人妻,人人澡人人爽秒播| 久久精品国产亚洲av香蕉五月| 日韩视频一区二区在线观看| 男人操女人黄网站| 日韩欧美国产在线观看| 欧美日韩精品网址| 大型黄色视频在线免费观看| 亚洲精品久久成人aⅴ小说| 久久久精品欧美日韩精品| 国产欧美日韩综合在线一区二区| 亚洲男人的天堂狠狠| 色婷婷久久久亚洲欧美| 成人国产一区最新在线观看| 国产一卡二卡三卡精品| 色精品久久人妻99蜜桃| 少妇的丰满在线观看| 99国产精品99久久久久| 一级a爱视频在线免费观看| 午夜免费成人在线视频| 日韩欧美三级三区| 免费人成视频x8x8入口观看| 日韩一卡2卡3卡4卡2021年| 精品第一国产精品| 亚洲色图av天堂| 韩国精品一区二区三区| 国产成年人精品一区二区| 亚洲国产精品成人综合色| 自拍欧美九色日韩亚洲蝌蚪91| 在线观看舔阴道视频| 91字幕亚洲| 国产黄a三级三级三级人| 麻豆国产av国片精品| 淫妇啪啪啪对白视频| 精品国产超薄肉色丝袜足j| 在线观看一区二区三区| 极品教师在线免费播放| 久久久国产成人免费| 国产精品一区二区免费欧美| 麻豆久久精品国产亚洲av| 亚洲成人国产一区在线观看| a级毛片在线看网站| 国产精品久久电影中文字幕| 国产免费男女视频| 精品久久久久久,| 免费女性裸体啪啪无遮挡网站| 久久久久久亚洲精品国产蜜桃av| 麻豆一二三区av精品| 身体一侧抽搐| 黄色视频不卡| 欧美色欧美亚洲另类二区 | 国产一区二区三区综合在线观看| 在线观看免费日韩欧美大片| 久久热在线av| 国产99白浆流出| 欧美激情极品国产一区二区三区| 无人区码免费观看不卡| 手机成人av网站| 一区福利在线观看| 欧美成人性av电影在线观看| 亚洲精品国产一区二区精华液| 成人永久免费在线观看视频| 涩涩av久久男人的天堂| 99国产精品一区二区蜜桃av| 日日干狠狠操夜夜爽| 老司机福利观看| 不卡av一区二区三区| 国产不卡一卡二| 国产av精品麻豆| 自拍欧美九色日韩亚洲蝌蚪91| 9热在线视频观看99| 男女做爰动态图高潮gif福利片 | 69av精品久久久久久| 97超级碰碰碰精品色视频在线观看| 中文字幕精品免费在线观看视频| 此物有八面人人有两片| 成人亚洲精品一区在线观看| 97超级碰碰碰精品色视频在线观看| netflix在线观看网站| 亚洲精品在线观看二区| 一二三四社区在线视频社区8| 久久久国产成人免费| 国产亚洲精品综合一区在线观看 | 女性生殖器流出的白浆| 男人操女人黄网站| 精品久久久久久久人妻蜜臀av | 国产三级黄色录像| 热99re8久久精品国产| 99国产精品一区二区三区| 一进一出好大好爽视频| 欧美成狂野欧美在线观看| 中文字幕av电影在线播放| 波多野结衣高清无吗| 制服丝袜大香蕉在线| 法律面前人人平等表现在哪些方面| 变态另类成人亚洲欧美熟女 | 国产单亲对白刺激| 美女国产高潮福利片在线看| 午夜免费观看网址| 国产亚洲av嫩草精品影院| 啦啦啦韩国在线观看视频| 国产野战对白在线观看| 亚洲国产精品成人综合色| 国产日韩一区二区三区精品不卡| 亚洲精品国产区一区二| 午夜福利免费观看在线| 极品人妻少妇av视频| 超碰成人久久| 色综合亚洲欧美另类图片| 一区在线观看完整版| 国产成人免费无遮挡视频| 亚洲国产毛片av蜜桃av| 国产成人精品久久二区二区91| 日本三级黄在线观看| 999久久久精品免费观看国产| 日韩三级视频一区二区三区| 97超级碰碰碰精品色视频在线观看| www日本在线高清视频| 午夜a级毛片| 欧美乱妇无乱码| 麻豆一二三区av精品| 村上凉子中文字幕在线| 免费人成视频x8x8入口观看| 亚洲在线自拍视频| 欧美日本亚洲视频在线播放| 亚洲av熟女| www.熟女人妻精品国产| 成人国产一区最新在线观看| 免费看美女性在线毛片视频| 亚洲国产中文字幕在线视频| 国产一区二区三区视频了| 久久久久久国产a免费观看| 一进一出抽搐动态| 熟妇人妻久久中文字幕3abv| 国产精品久久久人人做人人爽| 18禁黄网站禁片午夜丰满| 又黄又粗又硬又大视频| 亚洲精品久久成人aⅴ小说| 两性夫妻黄色片| 丝袜美足系列| 在线观看免费视频网站a站| 高清毛片免费观看视频网站| 50天的宝宝边吃奶边哭怎么回事| 日韩免费av在线播放| 国产精品av久久久久免费| 在线播放国产精品三级| 神马国产精品三级电影在线观看 | 午夜福利影视在线免费观看| 久久国产乱子伦精品免费另类| 久久人妻熟女aⅴ| 日韩欧美国产在线观看| 99久久久亚洲精品蜜臀av| 香蕉丝袜av| 欧美一级a爱片免费观看看 | 久久久久国产一级毛片高清牌| 国产精品98久久久久久宅男小说| 老汉色av国产亚洲站长工具| 91成人精品电影| 宅男免费午夜| 免费在线观看黄色视频的| 人妻丰满熟妇av一区二区三区| 国产精品一区二区精品视频观看| 国产99白浆流出| 色在线成人网| www.自偷自拍.com| 在线观看免费视频日本深夜| 亚洲va日本ⅴa欧美va伊人久久| 一区二区三区国产精品乱码| 激情视频va一区二区三区| 国产精品二区激情视频| 极品人妻少妇av视频| 欧美乱色亚洲激情| 亚洲片人在线观看| 久久久国产成人免费| 搡老熟女国产l中国老女人| 亚洲 国产 在线| 琪琪午夜伦伦电影理论片6080| 亚洲成国产人片在线观看| 欧美成人午夜精品| 人人妻人人爽人人添夜夜欢视频| 日韩大尺度精品在线看网址 | 国产精品av久久久久免费| 免费久久久久久久精品成人欧美视频| 曰老女人黄片| 老汉色av国产亚洲站长工具| 久久这里只有精品19| 亚洲五月色婷婷综合| 天天添夜夜摸| 亚洲第一电影网av| 亚洲成人久久性| 极品教师在线免费播放| 老司机午夜福利在线观看视频| 亚洲一区中文字幕在线| 99精品欧美一区二区三区四区| 免费不卡黄色视频| 精品不卡国产一区二区三区| 90打野战视频偷拍视频| 欧美另类亚洲清纯唯美| 久久欧美精品欧美久久欧美| 亚洲一码二码三码区别大吗| 啦啦啦韩国在线观看视频| 男女下面进入的视频免费午夜 | 少妇被粗大的猛进出69影院| 国产极品粉嫩免费观看在线| 夜夜看夜夜爽夜夜摸| 国产一区二区三区在线臀色熟女| 国产高清videossex| 波多野结衣av一区二区av| 欧美色视频一区免费| 午夜激情av网站| 少妇裸体淫交视频免费看高清 | 999精品在线视频| 中文字幕另类日韩欧美亚洲嫩草| 级片在线观看| 欧美日本亚洲视频在线播放| 波多野结衣一区麻豆| av免费在线观看网站| 亚洲aⅴ乱码一区二区在线播放 | 淫秽高清视频在线观看| 久久久久精品国产欧美久久久| 成人三级黄色视频| 免费在线观看黄色视频的| 亚洲人成网站在线播放欧美日韩| 亚洲av熟女| 国产精品电影一区二区三区| 国产蜜桃级精品一区二区三区| 女人被狂操c到高潮| 不卡一级毛片| 国产成人精品久久二区二区91| 人妻丰满熟妇av一区二区三区| avwww免费| 日本撒尿小便嘘嘘汇集6| 不卡一级毛片| 免费看十八禁软件| 99国产精品一区二区三区| 久久香蕉国产精品| 国产高清videossex| 欧美日本亚洲视频在线播放| 别揉我奶头~嗯~啊~动态视频| 又紧又爽又黄一区二区| 少妇裸体淫交视频免费看高清 | 欧美国产精品va在线观看不卡| 亚洲精品粉嫩美女一区| 国产片内射在线| 国产亚洲精品第一综合不卡| 99香蕉大伊视频| 精品国产亚洲在线| 国产av精品麻豆| 日日夜夜操网爽| 亚洲成人国产一区在线观看| 18禁黄网站禁片午夜丰满| 中文字幕av电影在线播放| 国产蜜桃级精品一区二区三区| 亚洲欧美激情综合另类| 女人被躁到高潮嗷嗷叫费观| 视频在线观看一区二区三区| 免费一级毛片在线播放高清视频 | 欧美日韩乱码在线| 国产精品一区二区免费欧美| 国产免费av片在线观看野外av| 三级毛片av免费| 999久久久国产精品视频| 99在线视频只有这里精品首页| 少妇 在线观看| 中文字幕另类日韩欧美亚洲嫩草| 怎么达到女性高潮| 一进一出好大好爽视频| svipshipincom国产片| 国产成人系列免费观看| 女人精品久久久久毛片| av天堂在线播放| 亚洲人成伊人成综合网2020| 国产一区二区三区视频了| 又大又爽又粗| 亚洲三区欧美一区| 久久天躁狠狠躁夜夜2o2o| 欧美绝顶高潮抽搐喷水| 成在线人永久免费视频| 欧美最黄视频在线播放免费| 国产私拍福利视频在线观看| 久久伊人香网站| 美女免费视频网站| 久久这里只有精品19| 丝袜在线中文字幕| 老鸭窝网址在线观看| 国产亚洲精品一区二区www| 成在线人永久免费视频| 夜夜躁狠狠躁天天躁| 亚洲熟妇中文字幕五十中出| 伦理电影免费视频| 99re在线观看精品视频| 国产麻豆69| 男人舔女人的私密视频| 在线观看免费午夜福利视频| 在线观看免费视频日本深夜| 电影成人av| 精品少妇一区二区三区视频日本电影| 中亚洲国语对白在线视频| 国产亚洲av高清不卡| 禁无遮挡网站| 亚洲精品在线美女| 久9热在线精品视频| 亚洲中文字幕一区二区三区有码在线看 | 亚洲国产精品999在线| av片东京热男人的天堂| 此物有八面人人有两片| 久久久精品国产亚洲av高清涩受| av中文乱码字幕在线| 国产精品亚洲av一区麻豆| 啦啦啦免费观看视频1| 亚洲avbb在线观看| 亚洲视频免费观看视频| 久久热在线av| 老汉色av国产亚洲站长工具| 黄色毛片三级朝国网站| 露出奶头的视频| 亚洲一区二区三区不卡视频| 日本免费a在线| 女性被躁到高潮视频| 母亲3免费完整高清在线观看| 亚洲色图av天堂| 高清毛片免费观看视频网站| 999精品在线视频| 99在线视频只有这里精品首页| 12—13女人毛片做爰片一| 久久这里只有精品19| 欧美av亚洲av综合av国产av| 国产亚洲欧美精品永久| 欧美黑人精品巨大| 级片在线观看| 亚洲熟妇中文字幕五十中出| 丝袜美足系列| 免费看美女性在线毛片视频| 老司机深夜福利视频在线观看| cao死你这个sao货| 久热爱精品视频在线9| 欧美成狂野欧美在线观看| 欧美一区二区精品小视频在线| 精品无人区乱码1区二区| 天天躁狠狠躁夜夜躁狠狠躁| 精品久久久久久成人av| 免费高清视频大片| 黄片播放在线免费| 美女高潮到喷水免费观看| 精品久久久精品久久久| 欧美在线一区亚洲| 亚洲人成网站在线播放欧美日韩|