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

    基于頻率控制的多約束風(fēng)電塔優(yōu)化方法

    2016-12-08 03:01:14東,東,
    關(guān)鍵詞:塔架固有頻率風(fēng)電

    程 耿 東, 徐 向 東, 劉 曉 峰

    ( 1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室, 遼寧 大連 116024;2.北京金風(fēng)科創(chuàng)風(fēng)電設(shè)備有限公司, 北京 100176 )

    ?

    基于頻率控制的多約束風(fēng)電塔優(yōu)化方法

    程 耿 東*1, 徐 向 東1, 劉 曉 峰2

    ( 1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室, 遼寧 大連 116024;2.北京金風(fēng)科創(chuàng)風(fēng)電設(shè)備有限公司, 北京 100176 )

    提出了基于頻率控制的多約束單管型風(fēng)電塔優(yōu)化方法.塔架簡(jiǎn)化成懸臂梁結(jié)構(gòu),其橫截面參數(shù)作為設(shè)計(jì)變量,以最小化材料體積為目標(biāo)函數(shù),按照將塔架設(shè)計(jì)成剛-剛或剛-柔或柔-柔不同類(lèi)型的要求設(shè)定塔架的固有頻率約束,采用專(zhuān)業(yè)軟件Bladed計(jì)算風(fēng)荷載,按照風(fēng)電塔規(guī)范考慮強(qiáng)度、穩(wěn)定性和疲勞等約束,這使得優(yōu)化結(jié)構(gòu)更符合實(shí)際設(shè)計(jì).考慮到采用Bladed荷載計(jì)算工作量很大,整個(gè)優(yōu)化過(guò)程分為幾個(gè)階段,在每個(gè)階段的開(kāi)始,以前一個(gè)階段的優(yōu)化設(shè)計(jì)作為初始設(shè)計(jì),并重新計(jì)算結(jié)構(gòu)荷載,在該階段內(nèi)于固定荷載下用移動(dòng)漸近線法(MMA)求解優(yōu)化問(wèn)題改進(jìn)設(shè)計(jì),所需的固有頻率、強(qiáng)度及疲勞約束靈敏度采用解析法獲得.對(duì)一現(xiàn)有塔架進(jìn)行優(yōu)化以說(shuō)明方法的有效性.根據(jù)塔架固有頻率和風(fēng)機(jī)工作轉(zhuǎn)速之間的關(guān)系,發(fā)展了高風(fēng)電塔的分類(lèi).在此基礎(chǔ)上,結(jié)合提出的優(yōu)化方法,可以幫助設(shè)計(jì)者判定在指定高度和機(jī)型下哪種類(lèi)型塔架更合適,為塔架概念設(shè)計(jì)提供有價(jià)值的參考.

    風(fēng)電塔;優(yōu)化設(shè)計(jì);多約束

    0 引 言

    隨著世界經(jīng)濟(jì)的快速發(fā)展,能源危機(jī)和環(huán)境污染問(wèn)題日趨嚴(yán)峻,風(fēng)能作為一種潛力巨大的清潔可再生能源,正面臨前所未有的發(fā)展機(jī)遇.進(jìn)入21世紀(jì)以來(lái),風(fēng)電產(chǎn)業(yè)以極快的速度發(fā)展.全球風(fēng)能理事會(huì)于2016年2月公布的數(shù)據(jù)統(tǒng)計(jì)[1]顯示,截至2015年,全球累計(jì)風(fēng)電裝機(jī)容量已達(dá)432.419 GW,這表明風(fēng)電行業(yè)已成為新世紀(jì)極具活力的新能源產(chǎn)業(yè).據(jù)估計(jì),到2020年風(fēng)能將占全球能源市場(chǎng)的5%[2].在丹麥,超過(guò)20%的電能來(lái)自風(fēng)力發(fā)電[3].近十年,中國(guó)的風(fēng)電產(chǎn)業(yè)和技術(shù)突飛猛進(jìn),2015年以30.5 GW的新增風(fēng)電裝機(jī)容量位居全球第一,占全球新增容量的48.4%[1].為了獲得更強(qiáng)更穩(wěn)定的風(fēng)能,不斷提高風(fēng)電塔的高度成為風(fēng)電行業(yè)的共識(shí).然而,隨著高度的增加,一方面風(fēng)電塔會(huì)變得更柔,在循環(huán)荷載下更容易發(fā)生共振和疲勞破壞,另一方面塔架的造價(jià)會(huì)大大提高.所以,高塔的設(shè)計(jì)優(yōu)化非常有意義.

    風(fēng)電塔設(shè)計(jì)涉及結(jié)構(gòu)力學(xué)、結(jié)構(gòu)動(dòng)力學(xué)和空氣動(dòng)力學(xué)等多學(xué)科.不少學(xué)者在風(fēng)電塔設(shè)計(jì)優(yōu)化方面做了研究.Negm等[4]對(duì)典型的圓筒式塔架以每段塔架的長(zhǎng)度、半徑以及壁厚為設(shè)計(jì)變量,考慮強(qiáng)度、最大位移、質(zhì)量和避免共振等約束的情況下,提出并比較了最小質(zhì)量、最大剛度、最大化剛度與質(zhì)量比值、指定固有頻率范圍以及最大化固有頻率等5種優(yōu)化目標(biāo),得出以最大化固有頻率為目標(biāo)函數(shù)效果最好.Uys等[5]對(duì)一幢45 m高的加肋圓筒型鋼塔進(jìn)行了以最小化成本為目標(biāo)的優(yōu)化設(shè)計(jì).其中以每段塔筒的厚度、環(huán)肋的尺寸和數(shù)量等為設(shè)計(jì)變量,依據(jù)材料費(fèi)用及其制造成本構(gòu)造目標(biāo)函數(shù),根據(jù)歐洲規(guī)范[6]計(jì)算沿塔的荷載,在考慮塔筒局部屈曲以及環(huán)肋屈曲的約束下進(jìn)行優(yōu)化.Yildirim等[7]依據(jù)ASCE[8]等設(shè)計(jì)規(guī)范計(jì)算設(shè)計(jì)地震荷載和風(fēng)荷載,考慮強(qiáng)度、疲勞、屈曲和自然頻率的約束,采用遺傳算法,對(duì)某1.5 MW風(fēng)機(jī)的鋼塔進(jìn)行了最輕化設(shè)計(jì).Zwick等[9]針對(duì)某一完整高度的海上格構(gòu)式風(fēng)電塔提出了一種迭代優(yōu)化方法,采用DNV-RP-C203[10]規(guī)范中提供的方法,分析在極限荷載狀態(tài)和疲勞極限狀態(tài)下結(jié)構(gòu)上關(guān)鍵位置的應(yīng)力,對(duì)其進(jìn)行減重設(shè)計(jì),體現(xiàn)了格構(gòu)式輕量化設(shè)計(jì)的潛力.Nicholson等[11]對(duì)風(fēng)電塔和基礎(chǔ)組合系統(tǒng)提出了以最小化材料成本和塔頂位移最小為雙目標(biāo)的優(yōu)化問(wèn)題,采用多目標(biāo)遺傳算法等多種優(yōu)化方法在Isight軟件自動(dòng)求解.Haghi等[12]也對(duì)塔架與基礎(chǔ)組合系統(tǒng)進(jìn)行了研究,提出了一種綜合空氣動(dòng)力學(xué)、流體力學(xué)、結(jié)構(gòu)力學(xué)和土力學(xué)等多學(xué)科的優(yōu)化方法,他們用此方法對(duì)一現(xiàn)有的SWT-3.6-107海上風(fēng)機(jī)的塔架和基礎(chǔ)進(jìn)行了優(yōu)化設(shè)計(jì),結(jié)果使支撐結(jié)構(gòu)的整體質(zhì)量減少了12.1%.

    在這些研究中,優(yōu)化過(guò)程使用的荷載雖然是根據(jù)規(guī)范計(jì)算,但往往只根據(jù)給定的一種風(fēng)速并在整個(gè)優(yōu)化過(guò)程中是固定的.但實(shí)際上,風(fēng)荷載有很多工況并且是一個(gè)隨機(jī)過(guò)程,塔架的荷載是與塔架的尺寸及動(dòng)力性能有關(guān)的,在優(yōu)化迭代過(guò)程中,隨著結(jié)構(gòu)的變化,荷載也會(huì)發(fā)生改變.為此本文提出基于頻率控制的多約束風(fēng)電塔優(yōu)化方法.

    1 多約束風(fēng)電塔優(yōu)化問(wèn)題

    1.1 優(yōu)化列式

    風(fēng)力發(fā)電機(jī)系統(tǒng)很復(fù)雜,包括葉片、輪轂、機(jī)艙、發(fā)電機(jī)、塔架、基礎(chǔ)等部分.本文中將塔架簡(jiǎn)化成非均勻的歐拉-伯努利梁,忽略塔內(nèi)平臺(tái)等非承重結(jié)構(gòu),塔頂結(jié)構(gòu)(如風(fēng)輪、機(jī)艙等)用集中質(zhì)量代替.塔架坐標(biāo)系如圖1所示[13],坐標(biāo)原點(diǎn)位于塔架中軸與基礎(chǔ)上表面的交點(diǎn),整個(gè)坐標(biāo)系不隨機(jī)艙轉(zhuǎn)動(dòng).本文提出優(yōu)化問(wèn)題是使風(fēng)電塔具有指定的固有頻率,并在滿足強(qiáng)度、穩(wěn)定性、疲勞以及橫截面尺寸邊界約束下,得到體積最小的設(shè)計(jì).若將風(fēng)電塔分成N個(gè)梁?jiǎn)卧?,則此問(wèn)題的優(yōu)化列式如下:

    find X=(x1x2…xn)T

    (1a)

    (1b)

    s.t.

    (1c)

    (1d)

    (1e)

    (1f)

    (1g)

    (1h)

    D-1≤0

    (1i)

    Xlower≤X≤Xupper

    (1j)

    式中:X是設(shè)計(jì)變量向量,V是塔架體積,Ae和Le分別代表第e個(gè)單元的橫截面面積和長(zhǎng)度.ωj和φj分別表示塔架第j階固有角頻率和相應(yīng)的特征向量.廣義特征值問(wèn)題(1c)中的K和M分別為對(duì)稱(chēng)的剛度陣和質(zhì)量陣.式(1d)保證各相應(yīng)的特征向量關(guān)于質(zhì)量陣正交歸一化.Xupper和Xlower分別為設(shè)計(jì)變量的上下限.約束條件(1c)~(1i)將在1.4節(jié)討論.

    Xf水平;Zf沿塔架軸方向垂直向上;Yf水平,使Xf、Yf、Zf符合右手定則

    圖1 風(fēng)電塔坐標(biāo)系

    Fig.1 Wind turbine tower coordinate system

    1.2 風(fēng)荷載計(jì)算

    風(fēng)荷載對(duì)風(fēng)電塔優(yōu)化設(shè)計(jì)非常重要,直接影響優(yōu)化結(jié)果的有效性.在許多現(xiàn)有的文章中,風(fēng)荷載簡(jiǎn)化成作用在塔頂?shù)囊粋€(gè)集中力和彎矩,或者是根據(jù)設(shè)計(jì)規(guī)范將風(fēng)荷載看成是沿塔架高度方向的分布式荷載,但均未考慮振動(dòng)的影響,而且只根據(jù)給定的一種風(fēng)速并在整個(gè)優(yōu)化過(guò)程中是固定的.實(shí)際上,風(fēng)荷載的方向和風(fēng)速是不斷變化的,具有極強(qiáng)的隨機(jī)性,因此專(zhuān)業(yè)設(shè)計(jì)時(shí),需要考慮大量不同工況并模擬風(fēng)的非平穩(wěn)隨機(jī)特性.此外,結(jié)構(gòu)荷載會(huì)受到由于塔架與風(fēng)機(jī)耦合振動(dòng)引起的動(dòng)力放大的影響,并且隨著結(jié)構(gòu)尺寸和頻率的變化,其荷載也是有差異的.Bladed是一個(gè)計(jì)算整機(jī)荷載的專(zhuān)業(yè)軟件,它通過(guò)建立正常湍流、極端風(fēng)速、極限持續(xù)陣風(fēng)、極端方向變化等多種風(fēng)模型,考慮塔影效應(yīng)、溫度等因素的影響,對(duì)正常發(fā)電、啟動(dòng)、停機(jī)、故障等多種設(shè)計(jì)荷載工況計(jì)算風(fēng)力發(fā)電機(jī)系統(tǒng)的內(nèi)部荷載(塔架荷載、葉片荷載和輪轂荷載等).

    本文使用的荷載通過(guò)Bladed軟件計(jì)算得到,考慮了700多個(gè)設(shè)計(jì)荷載工況,從中挑選出塔架的極限荷載,其中包括每個(gè)截面的剪力(Fx、Fy、Fz、Fxy)和彎矩(Mx、My、Mz、Mxy),以及用于疲勞計(jì)算的64×64的Markov矩陣.矩陣中的每個(gè)元素分別表示64種均值荷載下64種不同幅值的循環(huán)次數(shù).

    1.3 設(shè)計(jì)變量

    風(fēng)電塔橫截面參數(shù)作為設(shè)計(jì)變量X,依據(jù)不同的截面形狀,如圖2所示,定義相應(yīng)的設(shè)計(jì)變量.例如,對(duì)于圓環(huán)形截面塔架,每個(gè)截面的外半徑R=0.5D和厚度t作為設(shè)計(jì)變量;如果塔架截面形狀為圓角正方形,直邊的長(zhǎng)度l、圓角的半徑R和厚度t可作為設(shè)計(jì)變量;如果采用正六邊形截面的設(shè)計(jì)方案,設(shè)計(jì)變量可以是邊長(zhǎng)l和厚度t.這些設(shè)計(jì)變量都需滿足上下限約束,即Xlower≤X≤Xupper.

    圖2 不同截面形狀的設(shè)計(jì)變量

    本文選取圓環(huán)形截面塔架進(jìn)行優(yōu)化設(shè)計(jì),其橫截面面積A和慣性矩I可依據(jù)設(shè)計(jì)變量R和t寫(xiě)成

    A=π(R2-(R-t)2)

    (2)

    (3)

    根據(jù)每個(gè)截面的面積A、慣性矩I以及荷載,可以得到塔架的固有頻率以及塔殼的強(qiáng)度、穩(wěn)定性和疲勞等相關(guān)信息.

    1.4 設(shè)計(jì)約束及靈敏度分析

    1.4.1 塔架固有頻率與相應(yīng)的靈敏度分析 風(fēng)電塔的固有頻率會(huì)影響結(jié)構(gòu)動(dòng)力響應(yīng)和外部風(fēng)荷載,特別是一階頻率和二階頻率.當(dāng)這些低階頻率落在風(fēng)機(jī)工作轉(zhuǎn)速范圍時(shí),結(jié)構(gòu)荷載會(huì)大大增加,并且可能引發(fā)共振.因此,本文將塔架固有頻率控制作為結(jié)構(gòu)設(shè)計(jì)優(yōu)化時(shí)最重要的約束,即列式(1e)、(1f).

    (4)

    對(duì)于單特征值λj,若所有的設(shè)計(jì)變量同時(shí)變化,那么其線性增量Δλj為

    Δλj=λjΔX

    (5)

    式中:ΔX=(Δx1Δx2… Δxn)T,為設(shè)計(jì)變量增量向量;λj是特征值λj關(guān)于設(shè)計(jì)變量的梯度向量,即

    (6)

    1.4.2 靜強(qiáng)度 風(fēng)電塔的靜強(qiáng)度約束(1g)可以用第四強(qiáng)度理論來(lái)校核,即每個(gè)截面的等效應(yīng)力應(yīng)滿足

    (7)

    式中:fy,k和fy,d分別為材料強(qiáng)度的標(biāo)準(zhǔn)值和設(shè)計(jì)值,γm為材料分項(xiàng)系數(shù).等效應(yīng)力σv可由以下公式求得:

    (8)

    其中A為截面面積,Wxy和Wt分別為抗彎截面系數(shù)和抗扭截面系數(shù).若塔架橫截面是以外半徑R和厚度t為設(shè)計(jì)變量的圓環(huán)形截面,那么Wxy和Wt分別為

    (9)

    (10)

    為了簡(jiǎn)化公式,去掉等效應(yīng)力公式中的根號(hào),則靜強(qiáng)度約束條件可寫(xiě)為

    (11)

    其中σallowable=fy,d為材料強(qiáng)度設(shè)計(jì)值,那么關(guān)于設(shè)計(jì)變量R和t的靈敏度易于求得,即

    (12)

    (13)

    1.4.3 穩(wěn)定性約束 對(duì)于穩(wěn)定性約束(1h),根據(jù)設(shè)計(jì)規(guī)范EN 1993-1-6-2007[16],塔殼的正應(yīng)力設(shè)計(jì)值σx,Ed和實(shí)際軸向失穩(wěn)臨界應(yīng)力σx,Rd應(yīng)滿足

    (14)

    σx,Rd由以下公式計(jì)算:

    (15)

    這里,γm為材料分項(xiàng)系數(shù),根據(jù)設(shè)計(jì)規(guī)范[13]取值.折減系數(shù)χx根據(jù)以下公式計(jì)算:

    (16)

    上式中的相關(guān)參數(shù)可查閱規(guī)范EN 1993-1-6-2007[16]計(jì)算.

    由于以上計(jì)算塔架穩(wěn)定性的解析公式過(guò)于復(fù)雜,雖然可以推導(dǎo)解析的靈敏度,但無(wú)論從編程還是計(jì)算量考慮,本文采用差分法計(jì)算穩(wěn)定性關(guān)于設(shè)計(jì)變量的靈敏度.

    1.4.4 疲勞 對(duì)于疲勞損傷約束(1i),根據(jù)Palmgren-Miner累積損傷理論[17]以及材料的S-N曲線,鋼塔的疲勞損傷按下面公式計(jì)算:

    (17)

    其中D表示損傷,當(dāng)D>1 時(shí),表示損傷破壞.本文中用于計(jì)算疲勞損傷的數(shù)據(jù)來(lái)源于Bladed軟件計(jì)算得到的64×64的Markov矩陣(見(jiàn)1.2),所以這里Nf=64.其他參數(shù)的意義如下:Δσe為應(yīng)力變化幅值,這里僅考慮主要荷載My所產(chǎn)生的應(yīng)力,即Δσe=ΔMy/Wy;ne為應(yīng)力變化幅值Δσe對(duì)應(yīng)的循環(huán)次數(shù);m為材料S-N曲線的斜率;N0為材料S-N曲線上根據(jù)Δσe對(duì)應(yīng)的循環(huán)次數(shù);Δσ0為材料S-N曲線上根據(jù)Δσe對(duì)應(yīng)的應(yīng)力幅值.

    (18)

    1.5 優(yōu)化過(guò)程

    風(fēng)電塔優(yōu)化設(shè)計(jì)問(wèn)題(1)通過(guò)迭代求解,如圖3所示.第一步,初始化設(shè)計(jì)變量xi,i=1,2,…,n,并指定材料、幾何參數(shù)、有限元模型參數(shù),以及約束條件參數(shù).例如,按照將塔架設(shè)計(jì)成剛-剛或剛-柔或柔-柔不同類(lèi)型的要求指定塔架固有頻率上下限、強(qiáng)度和穩(wěn)定性的臨界應(yīng)力,以及設(shè)計(jì)變量上下限等.整個(gè)優(yōu)化過(guò)程包括3個(gè)循環(huán),即荷載循環(huán)、外循序和內(nèi)循環(huán),如圖3所示.

    圖3 迭代求解流程圖

    1.5.1 荷載循環(huán) 采用Bladed計(jì)算荷載是十分耗時(shí)的,所以僅當(dāng)結(jié)構(gòu)變化足夠大時(shí)才應(yīng)用Bladed 重新計(jì)算荷載.在優(yōu)化過(guò)程中,首先應(yīng)用Bladed計(jì)算初始設(shè)計(jì)的荷載,用于結(jié)構(gòu)強(qiáng)度、穩(wěn)定性及疲勞分析.然后,在隨后的優(yōu)化迭代過(guò)程中,結(jié)構(gòu)荷載均保持不變.在每一次迭代結(jié)束后更新設(shè)計(jì)變量,即X=X+ΔX,并且檢驗(yàn)是否需要重新計(jì)算荷載,其判定條件為在該荷載循環(huán)中先前所有迭代步設(shè)計(jì)變量變化的總和是否大于預(yù)先定義的值δ,即

    sum ∑Δx1,…,∑Δxn≥δ

    (19)

    隨著設(shè)計(jì)變量的不斷更新,在沒(méi)有收斂之前如果滿足了判定條件(19),那么需應(yīng)用更新后的設(shè)計(jì)重新計(jì)算結(jié)構(gòu)荷載,如圖3所示.通過(guò)這種方法,整個(gè)優(yōu)化過(guò)程被分成幾個(gè)荷載循環(huán),在每個(gè)荷載循環(huán)內(nèi)使用固定的結(jié)構(gòu)荷載.

    1.5.2 外循環(huán) 如圖3所示,外循環(huán)中的第一步是通過(guò)求解廣義特征值問(wèn)題(1c)得到塔架固有頻率,以及依據(jù)前一迭代步更新的設(shè)計(jì)變量和荷載計(jì)算強(qiáng)度、穩(wěn)定性和疲勞損傷.值得注意的是,由于這里并未出現(xiàn)重特征值的情況,故本文僅考慮特征值問(wèn)題(1c)單特征值的情況,那么相應(yīng)的特征向量φj,j=1,2,…,J也是唯一的.

    外循環(huán)中的第二步,是計(jì)算塔架頻率、強(qiáng)度、穩(wěn)定性及疲勞對(duì)設(shè)計(jì)變量的靈敏度,具體方法見(jiàn)1.1~1.4節(jié).外循環(huán)的第三步是一個(gè)求解子優(yōu)化問(wèn)題的內(nèi)循環(huán)(如圖3所示),其中將用到這些靈敏度.

    在外循環(huán)的第四步中,首先更新設(shè)計(jì)變量,即X=X+ΔX.然后檢驗(yàn)設(shè)計(jì)變量是否收斂,若設(shè)計(jì)變量增量向量ΔX的范數(shù)小于預(yù)先給定的一個(gè)小量ε,則表示收斂,那么優(yōu)化結(jié)束并得到最優(yōu)設(shè)計(jì)塔架(如圖3所示).否則,更新后的設(shè)計(jì)變量作為下一個(gè)外循環(huán)的初始設(shè)計(jì)繼續(xù)迭代優(yōu)化.要注意的是,在開(kāi)始下一個(gè)外循環(huán)開(kāi)始前需檢驗(yàn)是否需要重新計(jì)算結(jié)構(gòu)荷載,見(jiàn)1.5.1節(jié)中介紹.

    1.5.3 內(nèi)循環(huán) 根據(jù)1.1~1.4節(jié)中塔架頻率、強(qiáng)度、穩(wěn)定性及疲勞關(guān)于設(shè)計(jì)變量的靈敏度分析,在外循環(huán)中建立一個(gè)內(nèi)循環(huán)來(lái)求得最優(yōu)的設(shè)計(jì)變量增量Δxi,i=1,2,…,n,即求解下面的子優(yōu)化問(wèn)題:

    (20a)

    (20c)

    e=1,2,…,N

    e=1,2,…,N

    (20h)

    xlower≤xi+Δxi≤xupper;i=1,2,…,n

    (20i)

    注意,在這個(gè)子優(yōu)化問(wèn)題中只有設(shè)計(jì)變量增量Δxi,i=1,2,…,n是待優(yōu)化的.設(shè)計(jì)變量xi,i=1,2,…,n,塔架固有頻率ωj,j=1,2,以及相應(yīng)的靈敏度等在這個(gè)內(nèi)循環(huán)中是不變的.本文采用移動(dòng)漸近線法(MMA)[18]求解子優(yōu)化問(wèn)題(20).

    2 風(fēng)電塔設(shè)計(jì)的分類(lèi)

    風(fēng)電塔會(huì)長(zhǎng)期受到風(fēng)輪的轉(zhuǎn)動(dòng)引起的振動(dòng)荷載,所以避免塔架在循環(huán)荷載下發(fā)生共振是非常必要的.但是,對(duì)于高塔,在風(fēng)電系統(tǒng)啟動(dòng)過(guò)程中,通過(guò)采取適當(dāng)?shù)拇胧ㄖ鲃?dòng)控制,應(yīng)該允許風(fēng)機(jī)與結(jié)構(gòu)發(fā)生短暫的共振.Hu等[19-20]對(duì)某97 m高的風(fēng)電塔經(jīng)過(guò)兩年的檢測(cè)發(fā)現(xiàn)存在這樣的共振.這樣的風(fēng)電塔已經(jīng)不是傳統(tǒng)的剛性塔架.根據(jù)塔架固有頻率與風(fēng)機(jī)工作頻率范圍的關(guān)系,風(fēng)電塔可以分成剛-剛、剛-柔、柔-柔等設(shè)計(jì)類(lèi)型.為了更清晰地描述不同的塔架設(shè)計(jì),本文給出一個(gè)風(fēng)電塔設(shè)計(jì)的Campbell圖(圖4).在圖中,橫坐標(biāo)為風(fēng)輪的轉(zhuǎn)速,縱坐標(biāo)表示塔架的頻率,兩條水平線分別代表風(fēng)電塔的一階和二階頻率,而兩條斜線則分別表示風(fēng)輪轉(zhuǎn)速頻率和葉片通過(guò)塔架的頻率(對(duì)于三葉片風(fēng)機(jī)等于3倍風(fēng)輪轉(zhuǎn)速頻率).定義符號(hào)fij(i=1,3;j=1,2)分別表示兩條水平線和兩條斜線的4個(gè)交點(diǎn),如圖4所示.它們表示第j階塔架頻率ωj和i倍風(fēng)輪轉(zhuǎn)速發(fā)生共振時(shí)對(duì)應(yīng)的頻率.例如,f31表示塔架第一階頻率ω1與3倍風(fēng)輪轉(zhuǎn)速發(fā)生共振時(shí)對(duì)應(yīng)的頻率.圖中的Ωlower和Ωupper分別表示風(fēng)輪工作轉(zhuǎn)速的下限和上限.從圖4中可以看出,fij(i=1,3;j=1,2)的值很容易得到,并且它們之間存在一些特定的關(guān)系,如下:

    f11=ω1

    (21)

    (22)

    f12=ω2

    (23)

    (24)

    圖4 風(fēng)電塔設(shè)計(jì)Campbell圖

    若定義塔架的二階頻率與一階頻率之比為

    (25)

    對(duì)于單管型塔架一般α>3,因此fij(i=1,3;j=1,2)之間存在固定的關(guān)系:

    f31

    (26)

    依據(jù)上面的定義,可以清晰描述不同的塔架設(shè)計(jì)類(lèi)型,如表1所示.為了避免共振,所有的交點(diǎn)fij都須落在風(fēng)輪工作轉(zhuǎn)速范圍的外面.例如,剛-柔塔架設(shè)計(jì)意味著風(fēng)輪的轉(zhuǎn)速頻率的上限小于塔架的基頻ω1,葉片通過(guò)塔架的頻率下限(對(duì)于三葉片風(fēng)機(jī)等于3倍風(fēng)輪轉(zhuǎn)速頻率)大于基頻ω1,即基本的頻率約束為

    f11>Ωupper

    (27)

    同時(shí)次要約束為

    f31<Ωlower

    (28)

    即表示頻率點(diǎn)f31也不能位于風(fēng)機(jī)工作轉(zhuǎn)速范圍之間,這一塔架在啟動(dòng)過(guò)程,風(fēng)輪轉(zhuǎn)速頻率的3倍頻可以和塔的基頻發(fā)生共振,因此對(duì)葉頻來(lái)說(shuō)是一個(gè)柔性設(shè)計(jì).其他的塔架設(shè)計(jì)均在表1中列出.圖5中采用一維數(shù)軸圖直觀地表示了這些不同類(lèi)型塔架設(shè)計(jì).

    但需要注意的是,2a和3b這兩種設(shè)計(jì)只有當(dāng)風(fēng)機(jī)工作轉(zhuǎn)速上下限滿足以下關(guān)系時(shí)才出現(xiàn):

    3Ωlower>Ωupper

    (29)

    因此,在風(fēng)電塔設(shè)計(jì)時(shí),風(fēng)輪轉(zhuǎn)速范圍是否滿足條件(29)也是一個(gè)需要考慮的重要約束.

    表1 不共振塔架設(shè)計(jì)

    圖5 不同塔架設(shè)計(jì)中fij和風(fēng)機(jī)工作轉(zhuǎn)速的關(guān)系

    上述高風(fēng)電塔分類(lèi)方法更清晰地表達(dá)了不同類(lèi)型的塔架設(shè)計(jì)的固有頻率與風(fēng)機(jī)工作轉(zhuǎn)速之間的關(guān)系.結(jié)合本文提出的優(yōu)化方法,能夠說(shuō)明在指定高度和機(jī)型下哪種類(lèi)型塔架更合適,為塔架概念設(shè)計(jì)提供了有價(jià)值的參考.

    3 數(shù)值算例

    3.1 某2.5 MW風(fēng)機(jī)塔架優(yōu)化(算例1)

    這個(gè)算例中采用的風(fēng)電塔模型是參考某公司的鋼-混組合塔架,如圖6所示.塔高108.28 m,塔架下部為混凝土,上部為圓筒鋼塔,兩塔筒之間是剛性連接.鋼和混凝土的材料參數(shù)見(jiàn)表2.整個(gè)塔架的橫截面均為圓環(huán)形,所以設(shè)計(jì)變量為每個(gè)截面的外半徑R和壁厚t,如圖2(a)所示.

    圖6 塔架模型

    表2 鋼和混凝土的材料參數(shù)

    本文將整個(gè)塔架簡(jiǎn)化為非均勻的歐拉-伯努利梁,共分成46個(gè)梁?jiǎn)卧茼敳堪惭b的2.5 MW的風(fēng)機(jī)看成一個(gè)集中質(zhì)量加在梁的末端,塔頂風(fēng)機(jī)的總質(zhì)量為156 209 kg,包括3個(gè)葉片、機(jī)艙、發(fā)電機(jī)等.風(fēng)力發(fā)電機(jī)的工作轉(zhuǎn)速范圍為7.5~13.5 r/min(即0.125~0.225 Hz).

    表3給出了優(yōu)化設(shè)計(jì)時(shí)必要的參數(shù),包括塔架一階和二階頻率約束的上下限以及設(shè)計(jì)變量的邊界.在這個(gè)算例中,目標(biāo)函數(shù)為最小的鋼材體積,并且指定塔架的一階頻率為0.32 Hz,二階頻率大于6倍一階頻率.需注意的是,本算例中只對(duì)塔架的圓筒鋼塔部分進(jìn)行優(yōu)化,即塔架的混凝土部分的半徑和壁厚在優(yōu)化過(guò)程中是不變的.

    表4列出了本算例的優(yōu)化結(jié)果,包括每個(gè)優(yōu)化階段結(jié)束時(shí)塔架前兩階頻率值、塔架鋼段部分的體積以及鋼段體積減少量相對(duì)初始設(shè)計(jì)的百分比.

    表3 優(yōu)化參數(shù)

    從表4中可知,整個(gè)優(yōu)化過(guò)程被分為5個(gè)階段,鋼段部分的初始體積為29.811 0 m3,而最后得到的優(yōu)化設(shè)計(jì)的體積為24.491 6 m3,減少了17.84%.同時(shí),塔架的一階和二階頻率均滿足了約束條件.圖7中給出了鋼塔體積優(yōu)化迭代過(guò)程,并繪出了每一階段結(jié)束時(shí)塔架模型的簡(jiǎn)單輪廓.

    表4 優(yōu)化設(shè)計(jì)結(jié)果

    圖7 鋼塔體積優(yōu)化迭代過(guò)程

    另外,為了了解結(jié)構(gòu)荷載在優(yōu)化過(guò)程中的變化情況,提取出塔底(H=10 m)在每個(gè)優(yōu)化階段中的荷載數(shù)據(jù),如表5所示.從表中可看出,所有的荷載分量在各個(gè)優(yōu)化階段都是不同的,特別是Mxy和Fxy,變化比較大.

    表5 不同優(yōu)化階段中塔底的荷載(H=10 m)

    3.2 不同高度下各類(lèi)型塔架設(shè)計(jì)優(yōu)化(算例2)

    本算例結(jié)合上文中提出的高風(fēng)電塔分類(lèi)方法和塔架優(yōu)化方法,對(duì)不同高度的塔架進(jìn)行分類(lèi)及優(yōu)化,說(shuō)明在指定高度和機(jī)型下哪種類(lèi)型塔架更合適.在本算例中,所有候選塔架均為圓筒鋼塔,所以以塔架每個(gè)截面的半徑R和壁厚t作為設(shè)計(jì)變量.塔頂安裝的風(fēng)機(jī)型號(hào)與算例1中相同.考慮鋼筒運(yùn)輸和制造的限制,設(shè)定半徑允許的變化范圍為1.625~2.150 m,壁厚允許的變化范圍為0.01~0.10 m.

    根據(jù)第2章中對(duì)塔架分類(lèi)的討論,由于此風(fēng)機(jī)的工作轉(zhuǎn)速范圍為7.5~13.5 r/min,可以得到剛-剛設(shè)計(jì)、剛-柔設(shè)計(jì)和柔-柔設(shè)計(jì)的頻率約束條件,見(jiàn)表6.

    表6 不同塔架設(shè)計(jì)的頻率約束條件

    本算例分成兩個(gè)步驟進(jìn)行:第一步,對(duì)指定高度的風(fēng)電塔僅考慮頻率約束和設(shè)計(jì)變量邊界約束進(jìn)行優(yōu)化,即通過(guò)優(yōu)化列式(1a)~(1f)、(1j)先挑選出可能的優(yōu)化設(shè)計(jì);第二步,采用完整的優(yōu)化列式(1)對(duì)可能的優(yōu)化設(shè)計(jì)進(jìn)行優(yōu)化,檢驗(yàn)是否能得到最優(yōu)設(shè)計(jì).

    對(duì)80~170 m的10種不同高度的塔架進(jìn)行第一步的優(yōu)化,統(tǒng)計(jì)出各個(gè)高度塔架的各類(lèi)設(shè)計(jì)是否滿足頻率約束的結(jié)果,見(jiàn)表7.其中,“Yes”表示該優(yōu)化設(shè)計(jì)滿足頻率約束條件,“No”則表示不滿足.從表中可以看出,對(duì)于某一特定的高度,不是所有的塔架類(lèi)型都能滿足頻率約束的,而且隨著高度的增加具有潛力可行設(shè)計(jì)趨向于柔性設(shè)計(jì).

    表7 不同高度塔架第一步優(yōu)化統(tǒng)計(jì)結(jié)果

    通過(guò)第一步挑選之后,排除了每種高度下不滿足頻率約束的設(shè)計(jì)類(lèi)型.然后對(duì)每種高度下潛在的可行設(shè)計(jì)進(jìn)行考慮全部約束的第二步優(yōu)化.本算例中,只選擇120 m高的塔架作為例子進(jìn)行第二步優(yōu)化.從表7中得知,120 m高的塔架只有剛-柔設(shè)計(jì)(2a)和柔-柔設(shè)計(jì)(3a)為潛在的可行設(shè)計(jì),所以對(duì)這兩種設(shè)計(jì)采用完整的優(yōu)化列式(1)進(jìn)行進(jìn)一步的優(yōu)化.計(jì)算結(jié)果如表8所示,經(jīng)過(guò)第二步的優(yōu)化后,可行設(shè)計(jì)只剩剛-柔設(shè)計(jì)(2a),而柔-柔設(shè)計(jì)(3a)不能滿足ω1<0.125 Hz的頻率約束條件.同時(shí),還得到了120 m高度下的剛-柔設(shè)計(jì)塔架的最優(yōu)設(shè)計(jì),其體積為50.793 9 m3,一階和二階固有頻率分別為0.240 0 Hz和1.537 7 Hz.通過(guò)優(yōu)化得到的最優(yōu)設(shè)計(jì)可為風(fēng)電塔設(shè)計(jì)初始階段提供有一定價(jià)值的參考.

    表8 120 m塔架第二步優(yōu)化結(jié)果

    4 結(jié) 語(yǔ)

    本文提出了一種基于頻率控制并要求滿足多種幾何和力學(xué)約束的單管型風(fēng)電塔優(yōu)化方法.在此方法中,單管型風(fēng)電塔塔架簡(jiǎn)化成非均勻的歐拉-伯努利梁,整個(gè)塔架的材料可以是多樣的.該優(yōu)化方法對(duì)不同橫截面形狀的單管型風(fēng)電塔均適用,例如,圓環(huán)形截面、正六邊形截面等.從算例1中可以看到,在多種約束下塔架的材料體積得到了有效的減少,并且最優(yōu)化設(shè)計(jì)的塔架具有指定的固有頻率.

    本文的優(yōu)化方法有幾個(gè)顯著特點(diǎn):首先,在設(shè)計(jì)優(yōu)化時(shí)使用Bladed軟件計(jì)算的結(jié)構(gòu)荷載更符合實(shí)際,使得到的優(yōu)化設(shè)計(jì)更具有實(shí)際意義.其次,該方法把塔架固有頻率控制看成一個(gè)至關(guān)重要的約束條件.最后,本文采用基于梯度的優(yōu)化算法求解該風(fēng)電塔優(yōu)化問(wèn)題,并且除了穩(wěn)定性約束的靈敏度是采用差分法計(jì)算外,其他約束條件及目標(biāo)函數(shù)的靈敏度分析均采用解析的方法得到.

    另外,本文還根據(jù)塔架固有頻率和風(fēng)機(jī)工作轉(zhuǎn)速之間的關(guān)系,對(duì)高風(fēng)電塔的分類(lèi)進(jìn)行了進(jìn)一步的發(fā)展.結(jié)合本文提出的優(yōu)化方法,可以發(fā)現(xiàn)在指定高度和機(jī)型下哪種類(lèi)型塔架更合適,為塔架概念設(shè)計(jì)階段提供有價(jià)值的參考.

    [1] Global Wind Energy Council. Global wind statistics 2015 [EB/OL]. (2016-02-10) [2016-03-23]. http://www.gwec.net/wp-content/uploads/vip/GWEC-PRstats-2015_LR_corrected.pdf.

    [2] JIN Xin, LI Lang, JU Wen-bin,etal. Multibody modeling of varying complexity for dynamic analysis of large-scale wind turbines [J]. Renewable Energy, 2016, 90:336-351.

    [3] Hu Y, Baniotopoulos C, Yang J. Effect of internal stiffening rings and wall thickness on the structural response of steel wind turbine towers [J]. Engineering Structures, 2014, 81(81):148-161.

    [4] Negm H M, Maalawi K Y. Structural design optimization of wind turbine towers [J]. Computers & Structures, 2000, 74(6):649-666.

    [5] Uys P E, Farkas J, Jarmai K,etal. Optimisation of a steel tower for a wind turbine structure [J]. Engineering Structures, 2007, 29(7):1337-1342.

    [6] Standard Policy and Strategy Committee Eurocode 1: Actions on Structures - Part 1-4: General Actions - Wind Actions: EN 1991-1-4-2005 [S]. London: British Standard Institute, 2005.

    [7] Yildirim S, ?zkol I. Wind turbine tower optimization under various requirements by using genetic algorithm [J]. Engineering, 2010, 2(8):641-647.

    [8] American National Standards Institute. Minimum Design Loads for Buildings and Other Structures:ANSI/ASCE 7-88 [S]. Reston:ASCE, 2015:996.

    [9] Zwick D, Muskulus M, Moe G. Iterative optimization approach for the design of full-height lattice towers for offshore wind turbines [J]. Energy Procedia, 2012, 24:297-304.

    [10] Det Norske Veritas. Fatigue Design of Offshore Steel Structures, Recommended Practice: DNV-RP-C203 [S]. Norway:Det Norske Veritas, 2008.

    [11] Nicholson J C, Arora J S, Goyal D,etal. Multi-objective structural optimization of wind turbine tower and foundation systems using Isight:a process automation and design exploration software [C] // 10th World Congress on Structural and Multidisciplinary Optimization. Orlando:WCSMO, 2013.

    [12] Haghi R, Ashuri T, van der Valk P L C,etal. Integrated multidisciplinary constrained optimization of offshore support structures [J]. Journal of Physics:Conference Series, 2014, 555:012046.

    [13] Germanischer Lloyd WindEnergie GmbH. Guideline for the Certification of Wind Turbines Edition 2010 [Z]. Hamburg:Germanischer Lloyd WindEnergie GmbH, 2010.

    [14] DU Jian-bin, Olhoff N. Topological design of freely vibrating continuum structures for maximum values of simple and multiple eigenfrequencies and frequency gaps [J]. Structural and Multidisciplinary Optimization, 2007, 34(2):91-110.

    [15] 徐趙東,馬樂(lè)為. 結(jié)構(gòu)動(dòng)力學(xué)[M]. 北京:科學(xué)出版社, 2007.

    XU Zhao-dong, MA Le-wei. Structural Dynamics [M]. Beijing:Science Press, 2007. (in Chinese)

    [16] European Committee for Standardization. Eurocode 3:Design of Steel Structures, Part 1-6:Strength and Stability of Shell Structures:EN 1993-1-6-2007 [S]. Brussels:CEN, 2007.

    [17] Miner M A. Cumulative damage in fatigue [J]. Journal of Applied Mechanics, 1945, 12(3):159-164.

    [18] Svanberg K. The method of moving asymptotes-a new method for structural optimization [J]. International Journal for Numerical Methods in Engineering, 1987, 24(2):359-373.

    [19] HU Wei-hua, Th?ns S, Rohrmann R G,etal. Vibration-based structural health monitoring of a wind turbine system Part I:Resonance phenomenon [J]. Engineering Structures, 2015, 89:260-272.

    [20] HU Wei-hua, Th?ns S, Rohrmann R G,etal. Vibration-based structural health monitoring of a wind turbine system Part II:Environmental/operational effects on dynamic properties [J]. Engineering Structures, 2015, 89:273-290.

    Frequency-based optimization method for wind turbine tower under multiple constraints

    CHENG Geng-dong*1, XU Xiang-dong1, LIU Xiao-feng2

    ( 1.State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology,Dalian 116024, China;2.Goldwind Science & Technology Company Limited, Beijing 100176, China )

    A frequency-based optimization method for single tubular wind turbine tower under multiple constraints is proposed. The cross section parameters of tower are design variables, and minimizing the material volume of tower, which is considered as a cantilever structure is objective function. The natural frequencies of tower are constrained to ensure the tower being stiff-stiff or stiff-soft or soft-soft design. The wind loads are calculated by commercial software Bladed, and according to specification of wind turbine tower, the strength, stability and fatigue are included in constraints, these factors make the optimum structure conform to the practical design. Since the load evaluation by Bladed is very costly, the whole optimization process is broken into several stages. At the beginning of each stage optimization starts by using the previous optimum design as its initial design and the wind load of the initial tower design is recalculated. In each stage, optimization problem is solved by method of moving asymptotes (MMA) under the fixed wind load, and the sensitivities of the natural frequency, strength and fatigue constraints with respect to design variables are obtained by analytical method. A numerical example for improving the existing tower design demonstrates the optimization method. Furthermore, a classification of high wind turbine tower is developed based on the relationship between tower natural frequency and working frequency range of the wind turbine to be installed. By the above optimization method, it can help to decide that which type of the tower design is the potentially optimum candidate for specific height and turbine, and provide valuable suggestions for optimum design during concept design stage for wind turbine tower.

    wind turbine tower; optimization design; multiple constraints

    2016-04-06;

    2016-09-20.

    國(guó)家自然科學(xué)基金資助項(xiàng)目(11332004).

    程耿東*(1941-),男,教授,博士生導(dǎo)師,中國(guó)科學(xué)院院士,E-mail: chenggd@dlut.edu.cn;徐向東(1989-),男,碩士生,E-mail:xuxiangdong@mail.dlut.edu.cn.

    1000-8608(2016)06-0551-10

    TU359;TK83

    A

    10.7511/dllgxb201606001

    猜你喜歡
    塔架固有頻率風(fēng)電
    長(zhǎng)征六號(hào)甲火箭矗立在塔架旁
    上海航天(2022年5期)2022-12-05 01:55:46
    現(xiàn)場(chǎng)測(cè)定大型水輪發(fā)電機(jī)組軸系的固有頻率
    海上風(fēng)電躍進(jìn)隱憂
    能源(2018年6期)2018-08-01 03:42:00
    分散式風(fēng)電破“局”
    能源(2018年6期)2018-08-01 03:41:56
    風(fēng)電:棄風(fēng)限電明顯改善 海上風(fēng)電如火如荼
    能源(2018年8期)2018-01-15 19:18:24
    重齒風(fēng)電
    風(fēng)能(2016年12期)2016-02-25 08:46:38
    門(mén)式起重機(jī)塔架系統(tǒng)穩(wěn)定性分析
    雙塔式低塔架自平衡液壓提升裝置與吊裝技術(shù)
    風(fēng)力發(fā)電機(jī)設(shè)備塔架設(shè)計(jì)探析
    總溫總壓測(cè)頭模態(tài)振型變化規(guī)律研究
    亚洲国产欧美人成| 秋霞伦理黄片| 插逼视频在线观看| 日日撸夜夜添| 一级毛片我不卡| 热re99久久精品国产66热6| 成人一区二区视频在线观看| 欧美区成人在线视频| 中文欧美无线码| 少妇被粗大猛烈的视频| www.色视频.com| av又黄又爽大尺度在线免费看| 亚洲精品日韩av片在线观看| 精品国产乱码久久久久久小说| 国产黄色免费在线视频| av在线观看视频网站免费| 最黄视频免费看| 我的女老师完整版在线观看| 国产国拍精品亚洲av在线观看| 麻豆国产97在线/欧美| 亚洲国产精品一区三区| 少妇裸体淫交视频免费看高清| 欧美一区二区亚洲| 欧美zozozo另类| 一本一本综合久久| 热re99久久精品国产66热6| 97在线人人人人妻| 亚洲成人av在线免费| 最黄视频免费看| 免费观看性生交大片5| 久久精品久久久久久久性| kizo精华| 国产精品国产三级专区第一集| 欧美精品国产亚洲| 一级毛片我不卡| 久久精品国产a三级三级三级| 亚洲,欧美,日韩| 亚洲精品,欧美精品| 噜噜噜噜噜久久久久久91| 最近2019中文字幕mv第一页| 日韩欧美一区视频在线观看 | 男女下面进入的视频免费午夜| 欧美少妇被猛烈插入视频| 男男h啪啪无遮挡| 日本欧美国产在线视频| 成人美女网站在线观看视频| 在线观看av片永久免费下载| 日韩强制内射视频| 国产精品一二三区在线看| 免费看不卡的av| 22中文网久久字幕| 欧美性感艳星| 午夜精品国产一区二区电影| 在线天堂最新版资源| 日日撸夜夜添| 国产精品熟女久久久久浪| 免费黄频网站在线观看国产| 欧美精品一区二区大全| 搡老乐熟女国产| 麻豆精品久久久久久蜜桃| 韩国高清视频一区二区三区| 亚洲国产日韩一区二区| 少妇人妻 视频| av卡一久久| 内地一区二区视频在线| 日本免费在线观看一区| 色吧在线观看| 亚洲国产色片| 国产毛片在线视频| 26uuu在线亚洲综合色| 啦啦啦中文免费视频观看日本| 国产精品一二三区在线看| 伦精品一区二区三区| 国产精品无大码| 亚洲精品久久午夜乱码| 老师上课跳d突然被开到最大视频| 少妇的逼水好多| 国产黄色免费在线视频| 欧美性感艳星| 亚洲国产精品国产精品| 制服丝袜香蕉在线| 搡老乐熟女国产| 婷婷色综合www| 亚洲人成网站在线播| 九草在线视频观看| 爱豆传媒免费全集在线观看| 午夜激情福利司机影院| 欧美人与善性xxx| 亚洲美女黄色视频免费看| 少妇裸体淫交视频免费看高清| 熟女人妻精品中文字幕| 精品少妇黑人巨大在线播放| 国产一区二区在线观看日韩| 能在线免费看毛片的网站| 丝袜喷水一区| 这个男人来自地球电影免费观看 | 18禁裸乳无遮挡动漫免费视频| 99久久综合免费| 亚洲国产精品国产精品| 精品人妻偷拍中文字幕| 国产精品爽爽va在线观看网站| 国产 精品1| 久久毛片免费看一区二区三区| 久久久久久九九精品二区国产| 美女国产视频在线观看| 国产成人aa在线观看| 18+在线观看网站| 一区二区三区精品91| 色视频在线一区二区三区| 国产精品久久久久久久久免| h日本视频在线播放| 美女国产视频在线观看| 丝袜脚勾引网站| 国产成人午夜福利电影在线观看| 亚洲天堂av无毛| 久久韩国三级中文字幕| 精品一品国产午夜福利视频| 五月伊人婷婷丁香| 91精品国产国语对白视频| 狠狠精品人妻久久久久久综合| 日日啪夜夜爽| 人妻夜夜爽99麻豆av| 日本wwww免费看| 偷拍熟女少妇极品色| 极品教师在线视频| 国产av码专区亚洲av| 精品午夜福利在线看| 最近2019中文字幕mv第一页| 一个人看的www免费观看视频| 99久久精品热视频| 乱系列少妇在线播放| 亚洲精华国产精华液的使用体验| 黄片wwwwww| 日韩免费高清中文字幕av| 亚洲三级黄色毛片| 哪个播放器可以免费观看大片| 免费大片18禁| 高清在线视频一区二区三区| 男人舔奶头视频| 青春草亚洲视频在线观看| 久久99热6这里只有精品| 日韩国内少妇激情av| 最新中文字幕久久久久| 一区二区三区乱码不卡18| 成年av动漫网址| 亚洲欧美成人精品一区二区| 成人毛片a级毛片在线播放| 一级毛片 在线播放| 好男人视频免费观看在线| 成人二区视频| 中文字幕久久专区| 人体艺术视频欧美日本| 亚洲天堂av无毛| 十分钟在线观看高清视频www | 国产精品秋霞免费鲁丝片| 中文字幕精品免费在线观看视频 | 超碰av人人做人人爽久久| 在线亚洲精品国产二区图片欧美 | 少妇精品久久久久久久| av在线观看视频网站免费| 久久韩国三级中文字幕| 国产成人a区在线观看| 欧美精品人与动牲交sv欧美| 91在线精品国自产拍蜜月| av在线app专区| 国内少妇人妻偷人精品xxx网站| 我要看日韩黄色一级片| av黄色大香蕉| 婷婷色麻豆天堂久久| 亚洲欧美日韩卡通动漫| 精品一品国产午夜福利视频| 最新中文字幕久久久久| 国产中年淑女户外野战色| 91狼人影院| 精品久久久久久电影网| 秋霞在线观看毛片| 一级毛片黄色毛片免费观看视频| 欧美另类一区| 中文字幕人妻熟人妻熟丝袜美| av在线蜜桃| 欧美性感艳星| 夜夜爽夜夜爽视频| 3wmmmm亚洲av在线观看| 51国产日韩欧美| 少妇高潮的动态图| 久久久久人妻精品一区果冻| 久久久a久久爽久久v久久| 男人舔奶头视频| 18禁在线无遮挡免费观看视频| 亚洲人成网站高清观看| 日本午夜av视频| 中文字幕精品免费在线观看视频 | 欧美日韩视频高清一区二区三区二| 亚洲精品日本国产第一区| 蜜桃在线观看..| 观看av在线不卡| 国产精品久久久久成人av| 国产精品人妻久久久久久| 九九久久精品国产亚洲av麻豆| 欧美zozozo另类| 国产色婷婷99| 欧美日韩在线观看h| 18禁裸乳无遮挡免费网站照片| 精品视频人人做人人爽| 国产伦精品一区二区三区视频9| 老司机影院成人| 国产精品99久久久久久久久| 麻豆国产97在线/欧美| 欧美极品一区二区三区四区| 国产黄色免费在线视频| 国产老妇伦熟女老妇高清| 日韩国内少妇激情av| 精品午夜福利在线看| 久热久热在线精品观看| 亚洲精品亚洲一区二区| 99久久精品国产国产毛片| 国产伦在线观看视频一区| 国产精品国产三级国产av玫瑰| 精品人妻熟女av久视频| 插逼视频在线观看| 一二三四中文在线观看免费高清| 99久久精品一区二区三区| 日韩中字成人| 久久精品国产亚洲网站| 日韩欧美 国产精品| 亚洲欧美日韩卡通动漫| 国产毛片在线视频| 国产精品欧美亚洲77777| 99久久人妻综合| 亚洲av国产av综合av卡| 久久久色成人| 熟女电影av网| 一二三四中文在线观看免费高清| 免费看不卡的av| 十分钟在线观看高清视频www | 国产91av在线免费观看| 夜夜看夜夜爽夜夜摸| 97超视频在线观看视频| 特大巨黑吊av在线直播| 日本与韩国留学比较| 大香蕉97超碰在线| 波野结衣二区三区在线| 日本爱情动作片www.在线观看| 国产日韩欧美在线精品| 日日撸夜夜添| 欧美成人午夜免费资源| 国产一区二区三区av在线| 日韩亚洲欧美综合| 夫妻性生交免费视频一级片| 天堂8中文在线网| 日本欧美视频一区| 色婷婷av一区二区三区视频| 18禁动态无遮挡网站| 插逼视频在线观看| 国产亚洲午夜精品一区二区久久| 国产 精品1| 一个人看视频在线观看www免费| av在线观看视频网站免费| 一区二区三区四区激情视频| 国产伦理片在线播放av一区| 91久久精品国产一区二区三区| 内地一区二区视频在线| 免费观看a级毛片全部| 狂野欧美激情性bbbbbb| 1000部很黄的大片| 六月丁香七月| 日本黄色日本黄色录像| 丰满迷人的少妇在线观看| 又黄又爽又刺激的免费视频.| 国产一区二区三区av在线| 下体分泌物呈黄色| 大话2 男鬼变身卡| 毛片女人毛片| 欧美最新免费一区二区三区| 嫩草影院新地址| 亚洲av中文字字幕乱码综合| 九色成人免费人妻av| 亚洲av在线观看美女高潮| 自拍偷自拍亚洲精品老妇| 国产片特级美女逼逼视频| 亚洲经典国产精华液单| 久久国产亚洲av麻豆专区| 亚洲精华国产精华液的使用体验| 2018国产大陆天天弄谢| 久久久久久久久久久丰满| 国产午夜精品久久久久久一区二区三区| 日本一二三区视频观看| 99视频精品全部免费 在线| 校园人妻丝袜中文字幕| 人人妻人人爽人人添夜夜欢视频 | 街头女战士在线观看网站| a 毛片基地| 亚洲电影在线观看av| 亚洲aⅴ乱码一区二区在线播放| 久热这里只有精品99| 中国三级夫妇交换| 新久久久久国产一级毛片| 久久精品国产自在天天线| 最近手机中文字幕大全| 小蜜桃在线观看免费完整版高清| 亚洲精品乱码久久久久久按摩| 成人影院久久| 中文在线观看免费www的网站| 夫妻性生交免费视频一级片| 大码成人一级视频| 国产视频内射| 女人久久www免费人成看片| 99久久精品热视频| 婷婷色综合www| 亚洲精华国产精华液的使用体验| 一边亲一边摸免费视频| 久久久a久久爽久久v久久| 国产精品秋霞免费鲁丝片| 嘟嘟电影网在线观看| 亚洲精品久久久久久婷婷小说| 插阴视频在线观看视频| 全区人妻精品视频| 亚洲va在线va天堂va国产| 国产白丝娇喘喷水9色精品| 在线观看美女被高潮喷水网站| av专区在线播放| av又黄又爽大尺度在线免费看| 熟妇人妻不卡中文字幕| 久久精品久久久久久噜噜老黄| 一级黄片播放器| av网站免费在线观看视频| 亚洲欧美清纯卡通| 成年美女黄网站色视频大全免费 | 一级二级三级毛片免费看| 久久这里有精品视频免费| 又粗又硬又长又爽又黄的视频| 成年av动漫网址| 精华霜和精华液先用哪个| 建设人人有责人人尽责人人享有的 | 亚洲av国产av综合av卡| 秋霞在线观看毛片| 自拍偷自拍亚洲精品老妇| 国产中年淑女户外野战色| 亚洲av欧美aⅴ国产| 日本欧美国产在线视频| 成人毛片60女人毛片免费| 热99国产精品久久久久久7| 国产av精品麻豆| 久久午夜福利片| 国产精品.久久久| 免费观看无遮挡的男女| 看十八女毛片水多多多| 亚洲aⅴ乱码一区二区在线播放| 欧美老熟妇乱子伦牲交| 天堂8中文在线网| 高清视频免费观看一区二区| 精品少妇久久久久久888优播| 久久久欧美国产精品| 一级毛片电影观看| 中文资源天堂在线| 久久影院123| 欧美日本视频| 国产综合精华液| 中文资源天堂在线| 欧美精品一区二区免费开放| 日本黄色片子视频| av播播在线观看一区| 麻豆乱淫一区二区| 国产黄频视频在线观看| 少妇熟女欧美另类| 只有这里有精品99| 国产精品伦人一区二区| 国产精品福利在线免费观看| 亚洲精品456在线播放app| av播播在线观看一区| 天堂俺去俺来也www色官网| 永久免费av网站大全| 国产淫片久久久久久久久| 亚洲欧美精品自产自拍| 麻豆成人午夜福利视频| 男女啪啪激烈高潮av片| 在线免费观看不下载黄p国产| 在线看a的网站| 女性被躁到高潮视频| 国产成人精品久久久久久| 久久国内精品自在自线图片| av网站免费在线观看视频| 日韩 亚洲 欧美在线| 国产永久视频网站| 亚洲精品第二区| 中文字幕av成人在线电影| 国产人妻一区二区三区在| 亚洲av欧美aⅴ国产| 亚洲人与动物交配视频| 国产精品成人在线| 蜜桃亚洲精品一区二区三区| 亚洲美女视频黄频| 99久久人妻综合| 国产成人免费观看mmmm| 美女中出高潮动态图| 色哟哟·www| 日韩欧美一区视频在线观看 | av不卡在线播放| 永久免费av网站大全| 最近2019中文字幕mv第一页| 2022亚洲国产成人精品| 成人二区视频| 九九久久精品国产亚洲av麻豆| 国产爽快片一区二区三区| 国内揄拍国产精品人妻在线| 国产大屁股一区二区在线视频| 亚洲美女视频黄频| 国产色爽女视频免费观看| 欧美精品国产亚洲| 人妻 亚洲 视频| 亚洲av电影在线观看一区二区三区| 在线播放无遮挡| 在线观看美女被高潮喷水网站| 国产毛片在线视频| 卡戴珊不雅视频在线播放| 成人亚洲精品一区在线观看 | 亚洲性久久影院| 深爱激情五月婷婷| 女性生殖器流出的白浆| 一本久久精品| 成人毛片a级毛片在线播放| 99热6这里只有精品| 亚洲美女黄色视频免费看| 久久久久久久大尺度免费视频| 3wmmmm亚洲av在线观看| 中文字幕人妻熟人妻熟丝袜美| 久久久久久伊人网av| 极品少妇高潮喷水抽搐| 大片电影免费在线观看免费| 久久人人爽人人爽人人片va| 国产伦精品一区二区三区四那| 国产男女内射视频| 蜜桃在线观看..| 18禁在线播放成人免费| 性高湖久久久久久久久免费观看| 一二三四中文在线观看免费高清| 欧美最新免费一区二区三区| 在线 av 中文字幕| 日产精品乱码卡一卡2卡三| 如何舔出高潮| 色吧在线观看| 午夜视频国产福利| 91在线精品国自产拍蜜月| 欧美最新免费一区二区三区| 在线免费十八禁| 亚洲av国产av综合av卡| 久久久久久久国产电影| 精品午夜福利在线看| 午夜免费男女啪啪视频观看| 日本爱情动作片www.在线观看| 黄色一级大片看看| 插阴视频在线观看视频| 91午夜精品亚洲一区二区三区| 99热6这里只有精品| 国产 一区精品| 久久精品国产鲁丝片午夜精品| 七月丁香在线播放| 欧美变态另类bdsm刘玥| 久久青草综合色| 色网站视频免费| 国产精品久久久久久精品电影小说 | 国产伦精品一区二区三区四那| 国产男人的电影天堂91| 免费看不卡的av| 一区二区三区四区激情视频| 一级毛片黄色毛片免费观看视频| 亚洲国产精品成人久久小说| 国产精品福利在线免费观看| 九九爱精品视频在线观看| 大码成人一级视频| 国产伦理片在线播放av一区| videossex国产| 中文天堂在线官网| 如何舔出高潮| 精品午夜福利在线看| 欧美老熟妇乱子伦牲交| 99久久综合免费| 日韩欧美 国产精品| 久久精品国产亚洲网站| 精品人妻一区二区三区麻豆| 各种免费的搞黄视频| 亚洲人成网站在线观看播放| 国产免费视频播放在线视频| 1000部很黄的大片| 蜜桃亚洲精品一区二区三区| 一本色道久久久久久精品综合| 国产免费又黄又爽又色| 国产精品一区二区在线不卡| 成人亚洲精品一区在线观看 | 亚洲无线观看免费| 美女福利国产在线 | 精品久久久久久久久亚洲| 久久精品久久久久久久性| 亚洲精品日韩在线中文字幕| 在线免费十八禁| 久久久久久九九精品二区国产| 亚洲va在线va天堂va国产| 大陆偷拍与自拍| 日日摸夜夜添夜夜爱| 少妇高潮的动态图| 国产免费一区二区三区四区乱码| 亚洲人成网站在线观看播放| 亚洲av福利一区| 久久精品国产a三级三级三级| 精品人妻视频免费看| 久久精品久久久久久久性| 在线观看免费日韩欧美大片 | 国国产精品蜜臀av免费| 国产欧美亚洲国产| 在现免费观看毛片| 久久国产乱子免费精品| 久久久a久久爽久久v久久| 在线观看国产h片| 亚洲经典国产精华液单| 51国产日韩欧美| 26uuu在线亚洲综合色| 青春草视频在线免费观看| 日本av手机在线免费观看| 亚洲av二区三区四区| 国产成人aa在线观看| 男女下面进入的视频免费午夜| 又黄又爽又刺激的免费视频.| 国产精品欧美亚洲77777| 色综合色国产| 精品一区二区三区视频在线| 一区在线观看完整版| 国产高清三级在线| 偷拍熟女少妇极品色| 全区人妻精品视频| 超碰97精品在线观看| 日韩三级伦理在线观看| 亚洲在久久综合| 一本一本综合久久| 五月伊人婷婷丁香| 国产淫语在线视频| 卡戴珊不雅视频在线播放| 国内少妇人妻偷人精品xxx网站| 美女视频免费永久观看网站| 国产男人的电影天堂91| 久久久久国产网址| 国产精品99久久99久久久不卡 | 男女边摸边吃奶| 99热这里只有精品一区| 2018国产大陆天天弄谢| 性高湖久久久久久久久免费观看| 麻豆成人午夜福利视频| 激情 狠狠 欧美| 国产高清不卡午夜福利| 99热这里只有是精品在线观看| 国内揄拍国产精品人妻在线| 亚洲av中文字字幕乱码综合| 好男人视频免费观看在线| freevideosex欧美| 新久久久久国产一级毛片| 丰满迷人的少妇在线观看| 欧美日韩一区二区视频在线观看视频在线| 亚洲精品成人av观看孕妇| 日韩制服骚丝袜av| 国产伦理片在线播放av一区| av不卡在线播放| 3wmmmm亚洲av在线观看| 日韩一区二区视频免费看| 美女福利国产在线 | 天天躁日日操中文字幕| 国产爱豆传媒在线观看| 九九爱精品视频在线观看| 伦理电影免费视频| 精品国产露脸久久av麻豆| 国产无遮挡羞羞视频在线观看| 最近最新中文字幕免费大全7| 国产午夜精品一二区理论片| 成人无遮挡网站| 汤姆久久久久久久影院中文字幕| 日日啪夜夜爽| 在线观看美女被高潮喷水网站| h日本视频在线播放| 久久人人爽人人片av| 一本—道久久a久久精品蜜桃钙片| 尾随美女入室| 中文乱码字字幕精品一区二区三区| 国产真实伦视频高清在线观看| 欧美成人午夜免费资源| 成人18禁高潮啪啪吃奶动态图 | 亚洲国产毛片av蜜桃av| 日本av手机在线免费观看| 欧美精品人与动牲交sv欧美| 九九爱精品视频在线观看| 青春草亚洲视频在线观看| 寂寞人妻少妇视频99o| 中文字幕亚洲精品专区| 亚洲精品国产色婷婷电影| 亚洲国产高清在线一区二区三| 久久精品国产亚洲网站| xxx大片免费视频| 看十八女毛片水多多多| 国产乱来视频区| 自拍欧美九色日韩亚洲蝌蚪91 | 免费播放大片免费观看视频在线观看| 成人综合一区亚洲| 亚洲高清免费不卡视频| 91aial.com中文字幕在线观看| 高清在线视频一区二区三区| 色视频在线一区二区三区| av网站免费在线观看视频| 中文字幕人妻熟人妻熟丝袜美| 成人漫画全彩无遮挡| 午夜福利网站1000一区二区三区| 晚上一个人看的免费电影| a级毛片免费高清观看在线播放| 国产成人freesex在线| 少妇人妻精品综合一区二区| 国产成人免费观看mmmm| 99国产精品免费福利视频| 国产片特级美女逼逼视频| 日本vs欧美在线观看视频 |