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

基于渦格法的螺旋槳剖面拱線設(shè)計(jì)

2016-09-02 03:08:16饒志強(qiáng)楊晨俊
中國(guó)艦船研究 2016年4期
關(guān)鍵詞:方法設(shè)計(jì)

饒志強(qiáng),楊晨俊

1上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海2002402高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,上海200240

基于渦格法的螺旋槳剖面拱線設(shè)計(jì)

饒志強(qiáng)1,2,楊晨俊1,2

1上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海200240
2高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,上海200240

船舶螺旋槳的葉剖面形狀與其負(fù)荷形式、效率、空泡等水動(dòng)力性能密切相關(guān)。剖面拱線與剖面的負(fù)荷分布有直接關(guān)系。為了設(shè)計(jì)給定負(fù)荷分布的剖面拱線,提出用升力分布定義螺旋槳剖面拱線的負(fù)荷分布,在此基礎(chǔ)上對(duì)升力分布進(jìn)行參數(shù)化表達(dá),并進(jìn)行拱線設(shè)計(jì)。拱線的設(shè)計(jì)采用Newton-Raphson迭代方法使拱線的升力分布滿足給定值,拱線的性能計(jì)算采用基于薄翼理論的二維渦格法。計(jì)算表明:渦格法與Newton-Raphson迭代方法相結(jié)合能夠精確并快速設(shè)計(jì)滿足給定升力分布的拱線。最后,修改升力分布的部分參數(shù),并由此設(shè)計(jì)了SJ系列拱線供剖面設(shè)計(jì)參考。

螺旋槳葉剖面;拱線設(shè)計(jì);升力分布;渦格法;SJ系列拱線

網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/42.1755.TJ.20160729.0945.012.html期刊網(wǎng)址:www.ship-research.com

引用格式:饒志強(qiáng),楊晨俊.基于渦格法的螺旋槳剖面拱線設(shè)計(jì)[J].中國(guó)艦船研究,2016,11(4):38-43.

RAO Zhiqiang,YANG Chenjun.Camber lines design approach of propeller blade sections based on the vortex lattice method[J].Chinese Journal of Ship Research,2016,11(4):38-43.

0 引 言

剖面形狀對(duì)螺旋槳的水動(dòng)力性能有直接影響。在剖面設(shè)計(jì)中,NACA系列的剖面[1]因其優(yōu)良的抗空化性能而被設(shè)計(jì)者采用,比如NACA-66(mod)厚度分布和a=0.8拱線組合。但NACA剖面提出的時(shí)間很早,之后螺旋槳剖面(或稱翼剖面)設(shè)計(jì)也已經(jīng)發(fā)展了幾種方法,這些方法主要有單點(diǎn)設(shè)計(jì)法、多點(diǎn)設(shè)計(jì)法、面元法和優(yōu)化設(shè)計(jì)方法等[2]。單點(diǎn)法以Lighthill[3]方法為代表,做法是給定一個(gè)攻角和壓力分布,通過保角變換設(shè)計(jì)剖面。多點(diǎn)法以Eppler[4]方法為代表,它是在剖面表面分段給出不同攻角的速度分布、設(shè)計(jì)剖面。由于分段處理,它為改善剖面在設(shè)計(jì)和非設(shè)計(jì)工況下的空泡性能提供了潛力。一些學(xué)者對(duì)Eppler方法在剖面設(shè)計(jì)中的應(yīng)用進(jìn)行了研究。Eppler等[5]和Shen等[6]將該方法推廣到水翼設(shè)計(jì)中。王大政等[7-8]基于Eppler方法提出了可以控制最大厚度和拱度的剖面設(shè)計(jì)方法,并研究了剖面參數(shù)對(duì)空泡特性的影響。Zhou等[9]將Eppler方法設(shè)計(jì)的剖面直接應(yīng)用到螺旋槳上。面元設(shè)計(jì)方法是對(duì)初始剖面進(jìn)行擾動(dòng)并做迭代計(jì)算,使其滿足給定壓力分布。處理方法主要有2種:一種是直接在剖面上布置控制點(diǎn),并對(duì)其進(jìn)行攝動(dòng)和迭代計(jì)算,這些工作以蘇玉民等[10]和蔡昊鵬等[11]的工作為代表;另一種是用B樣條擬合并擾動(dòng)擬合曲線控制點(diǎn),譚廷壽等[12]和李俊華等[13]都采用該方法設(shè)計(jì)剖面,所不同的是前者擬合拱線而厚度不變,后者對(duì)表面進(jìn)行擬合。前述設(shè)計(jì)方法的共同點(diǎn)都是需要按給定壓力或者速度分布進(jìn)行設(shè)計(jì),但在設(shè)計(jì)初始階段合理的壓力和速度分布可能并不容易給定。給定值有缺陷可能對(duì)幾何或計(jì)算壓力分布的光順性造成影響,文獻(xiàn)[11]討論了這個(gè)問題。另外,給定壓力分布直接設(shè)計(jì)剖面可能會(huì)產(chǎn)生非常厚的剖面,文獻(xiàn)[13]中提出該問題并通過減小空泡裕度解決剖面過厚的問題。面對(duì)設(shè)計(jì)中的問題,優(yōu)化方法成為剖面設(shè)計(jì)的另一種選擇。黃勝等[14]和常欣等[15]在剖面上下表面疊加基函數(shù)并采用多目標(biāo)優(yōu)化算法設(shè)計(jì)剖面。曾志波[16]和龍文等[17]基于Eppler方法對(duì)螺旋槳剖面進(jìn)行了優(yōu)化設(shè)計(jì)。在優(yōu)化過程中,文獻(xiàn)[14-15]和文獻(xiàn)[16]分別對(duì)基函數(shù)的疊加量或者剖面最大厚度做了約束。

分析部分針對(duì)空泡性能的研究成果,特別是Eppler優(yōu)化設(shè)計(jì)結(jié)果[16-17],可以看出滿足約束的最優(yōu)設(shè)計(jì)剖面有如下特征:最大厚度位置向?qū)н吳耙疲黾尤~背和葉面片空泡裕度,最大拱度位置向隨邊后移將負(fù)荷向隨邊移動(dòng)[16]。本文作者曾針對(duì)固定厚度分布不變,采用面元法設(shè)計(jì)了具有不同拱線的剖面,并考察載荷分布對(duì)空泡性能的影響[18]。計(jì)算表明:載荷的移動(dòng)對(duì)剖面空泡斗的上下位置有明顯的影響。如此可以理解:載荷的移動(dòng)在一定程度上能夠改變空泡性能。這樣問題的關(guān)鍵變成了如何設(shè)計(jì)具有給定載荷的拱線。如能實(shí)現(xiàn)給定載荷設(shè)計(jì),拱線的幾何設(shè)計(jì)問題實(shí)際上也就變成了載荷設(shè)計(jì)。

本文將針對(duì)上述問題對(duì)載荷分布用升力分布來定義,對(duì)其進(jìn)行參數(shù)化表達(dá),并采用Newton-Raphson迭代方法設(shè)計(jì)滿足給定升力分布要求的拱線,拱線性能的計(jì)算采用二維渦格法。最后,修改給定升力分布的部分參數(shù),設(shè)計(jì)SJ系列拱線供參考。

1 拱線設(shè)計(jì)方法

1.1升力分布的定義

升力分布(CLx)定義為升力在拱線弦向的密度,具體定義為微段弦長(zhǎng)上的升力系數(shù)與微段弦長(zhǎng)的比,如式(1)所示。

式中:ΔCL為對(duì)應(yīng)的升力系數(shù);Δc為微段弦長(zhǎng)。由于拱線沒有厚度,因此可用薄翼理論來求解。

在零厚度薄翼理論中,薄翼的升力系數(shù)(CL)為[19]

式中:L為剖面總升力;ρ為流體密度;V為來流速度;c為弦長(zhǎng);xˉ為無因次的弦向坐標(biāo);γ()xˉ為弦向環(huán)量密度。

根據(jù)式(1)升力分布的定義,容易得到零厚度拱線的升力分布為

即在基于薄翼理論的拱線設(shè)計(jì)中,給定升力分布可以理解為2倍的給定環(huán)量密度在弦向的分布。采用二維渦格法可數(shù)值求解拱線的環(huán)量密度和升力分布。另外,通過升力分布的定義可知,剖面升力系數(shù)可由升力分布在弦向上的積分得到。

1.2升力分布的參數(shù)化表達(dá)

獲得升力分布的定義后,需要給定升力分布的形式。這里采用5段線來給定升力分布的形式。如圖1所示,升力分布由3條直線(L1~L3)和2條相切圓弧(C1和C2)構(gòu)成,整條升力分布曲線各段光順連接。

圖1 升力分布的參數(shù)化表達(dá)Fig.1 Parametric representation of lift distribution

用5個(gè)參數(shù)可確定整個(gè)升力分布曲線,包括:

1)第1個(gè)角點(diǎn)A的位置xa;

2)第2個(gè)角點(diǎn)B的位置xb;

3)圓弧C1的相對(duì)半徑,r1=Ar1r1 max;

4)圓弧C2的相對(duì)半徑,r2=Ar2r2 max;

5)第3段直線L3的斜率角θ3。

式中:r1 max和r2 max分別為2段圓弧的最大半徑,此時(shí),L1和L2這2條直線長(zhǎng)度為0。參數(shù) Ar1和 Ar2為最大半徑的比率,取值在0~1之間,用于控制C1和C2的半徑。5個(gè)參數(shù)外加給定升力系數(shù)CL,即可通過迭代的方法獲得整條升力分布曲線。在這種參數(shù)化定義中,拱線導(dǎo)邊和隨邊的升力分布都為0,即載荷為0。第3段線的角度θ3可以控制升力分布的前移或者后移。如果調(diào)整A角和B角點(diǎn)位置,以及圓弧C1和C2的半徑大小還可以對(duì)最大升力分布的位置進(jìn)一步進(jìn)行移動(dòng)。這是一種直接對(duì)載荷分布進(jìn)行控制的方法。盡管調(diào)整環(huán)量分布可以達(dá)到載荷調(diào)整的效果,但升力分布的定義更直觀。

1.3拱線設(shè)計(jì)過程

拱線的設(shè)計(jì)以升力分布為目標(biāo)。首先,選擇一個(gè)初始拱線,在初始拱線上設(shè)置M+1個(gè)控制點(diǎn),如圖2所示。

圖2 拱線上控制點(diǎn)布置Fig.2 Control nodes arrangement on camber line

圖中,x/c和f/c分別為無因次弦長(zhǎng)位置和拱度弦長(zhǎng)比,導(dǎo)邊和隨邊為固定控制點(diǎn),其他為可動(dòng)控制點(diǎn)。對(duì)攻角和可動(dòng)控制點(diǎn)進(jìn)行攝動(dòng),以獲得攝動(dòng)量對(duì)弦向某些位置(稱為評(píng)價(jià)位置)升力分布的影響。本文中評(píng)價(jià)位置與可動(dòng)控制點(diǎn)的位置重合。由于增加了攻角的攝動(dòng),因此在導(dǎo)邊和第一個(gè)可動(dòng)控制點(diǎn)之間增加一個(gè)評(píng)估位置來構(gòu)建M階雅可比矩陣,如式(4)所示。

對(duì)新的拱線和攻角再進(jìn)行攝動(dòng),直到評(píng)價(jià)位置上的升力分布計(jì)算值和給定值達(dá)到容差為止,本文容差為 1.0×10-4,這個(gè)設(shè)計(jì)過程即為Newton-Raphson迭代過程。在設(shè)計(jì)過程中,攝動(dòng)初始拱線拱度為0,初始攻角為0.5°。

拱線的計(jì)算方法采用二維渦格法,即將拱線分成若干單元(本文取120個(gè)單元)。在每個(gè)單元的1/4長(zhǎng)度處布置一個(gè)點(diǎn)渦,并將3/4位置作為單元控制點(diǎn)的位置,物面不可穿透邊界條件在單元控制點(diǎn)上滿足。渦格法的具體實(shí)施可見文獻(xiàn)[19-21]的相關(guān)章節(jié)。

這里舉一個(gè)例子詳細(xì)說明迭代結(jié)果。在這個(gè)例子中,θ3=0°,xa=0.1,xb=0.85,Ar1=0.2,Ar2=0.2,設(shè)計(jì)升力系數(shù)為CL=1.0。總控制點(diǎn)個(gè)數(shù)為31個(gè)(對(duì)應(yīng)雅可比矩陣階數(shù)M=30),其弦向位置由升力分布的5段線決定。圖3為設(shè)計(jì)過程中拱線隨迭代次數(shù)的變化。可以看出,迭代3次后計(jì)算結(jié)果收斂。圖4為給定升力分布、設(shè)計(jì)拱線的升力分布以及評(píng)價(jià)位置處升力分布的比較。可以看出,設(shè)計(jì)結(jié)果和給定的目標(biāo)十分吻合。這說明本方法迭代計(jì)算量少,設(shè)計(jì)結(jié)果精確。

圖3 設(shè)計(jì)過程中拱線的變化Fig.3 Camber lines during iterations

圖4 給定和設(shè)計(jì)拱線的升力分布比較Fig.4 Comparison of given and designed lift distributions

在文獻(xiàn)[18]中也討論過剖面設(shè)計(jì),但剖面的性能計(jì)算基于面元法。如果將剖面的厚度限定得很小,攝動(dòng)剖面的拱線,也能通過Newton-Raphson方法設(shè)計(jì)出來。但由于厚度的影響不能完全忽略,實(shí)際計(jì)算時(shí)發(fā)現(xiàn)采用面元法的設(shè)計(jì)過程不太穩(wěn)定,設(shè)計(jì)拱線也可能會(huì)出現(xiàn)不光順的情況。這可能是由于文獻(xiàn)[18]中升力分布的定義和計(jì)算是在剖面的上下表面對(duì)應(yīng)位置插值造成的,而本文渦格法的計(jì)算中升力分布和環(huán)量密度直接相關(guān),無需插值計(jì)算。由于渦格法能夠獲得穩(wěn)定、光順的拱線,是一種比較理想的方法,這里僅討論渦格法的計(jì)算結(jié)果。

2 SJ系列拱線

針對(duì)不同的第3段直線斜率角θ3,給定不同的升力分布曲線,設(shè)計(jì)了系列拱線,其中,升力分布曲線的其他參數(shù)仍為:xa=0.1,xb=0.85,Ar1=0.2,Ar2=0.2,CL=1.0。這一系列的拱線命名為SJ系列拱線,其中SJ后第1位數(shù)字代表斜率角θ3的正負(fù)(0為負(fù),1為正),后2個(gè)數(shù)字代表角度。如SJ015,代表θ3=-15°時(shí)的設(shè)計(jì)拱線。圖5為部分SJ系列拱線的計(jì)算升力分布曲線。由于設(shè)計(jì)結(jié)果與給定升力分布十分吻合,故沒有提供給定升力分布的結(jié)果。

圖5 SJ系列拱線的計(jì)算升力分布Fig.5 Calculated lift distributions of SJ series camber lines

圖6為對(duì)應(yīng)的部分SJ系列的無因次拱線圖,其中f/fmax為弦向拱度與最大拱度比值。從圖中可以看出,拱線的最大拱度位置隨升力分布的前移或后移而分別前后移動(dòng)。

圖6 SJ系列無因次拱線圖Fig.6 Non-dimensional SJ series camber lines

表1為SJ系列拱線的攻角α和拱度弦長(zhǎng)比fmax/c。表2為SJ系列拱線的型值表,可供后續(xù)的設(shè)計(jì)參考。

表1 SJ系列拱線攻角和拱度弦長(zhǎng)比(CL=1.0)Tab.1 Angle of attack and camber ratio of SJ series camber lines(CL=1.0)

表2 SJ系列拱線型值表Tab.2 Offset table of SJ series camber lines

3 結(jié) 論

本文對(duì)螺旋槳剖面拱線的升力分布進(jìn)行了定義,并提出了一種升力分布的參數(shù)化方法。采用Newton-Raphson方法設(shè)計(jì)了給定升力分布的系列拱線。主要結(jié)論和工作如下:

1)對(duì)螺旋槳剖面拱線的載荷用升力分布來表達(dá),并提出了一種升力分布的參數(shù)化形式。

2)結(jié)合Newton-Raphson迭代方法和渦格法可以快速、精確地設(shè)計(jì)給定升力分布的拱線。

3)設(shè)計(jì)了SJ系列拱線并給出了幾何型值供設(shè)計(jì)參考。

[1]ABBOTT I H,VON DOENHOFF A E.Theory of wingsections:including a summary of airfoil data[M].Dover,DE:Dover Publications Inc.,1959.

[2]王獻(xiàn)孚.船用翼理論[M].北京:國(guó)防工業(yè)出版社,1998.

[3]LIGHTHILL M J.A new method of two-dimensional aerodynamic design:M2112[R].London,UK:Aeronautical Research Council,1945.

[4]EPPLER R,SOMERS D M.A computer program for thedesignandanalysisoflowspeedairfoils:NASA-TM-80210[R].Hampton,VA:NASA Langley Research Center,1980.

[5]EPPLER R,SHEN Y T.Wing sections for hydrofoils-part 1:symmetrical profiles[J].Journal of Ship Research,1979,23(3):209-217.

[6]SHEN Y T,EPPLER R.Wing sections for hydrofoils-part 2:nonsymmetrical profiles[J].Journal of Ship Research,1981,25(3):191-200.

[7] 王大政,王言英,柴樹紅.新型葉剖面設(shè)計(jì)及葉剖面參數(shù)對(duì)空泡特性影響的研究——(1)新型葉剖面設(shè)計(jì)[J].水動(dòng)力學(xué)研究與進(jìn)展,1999,14(2):189-200.

WANG Dazheng,WANG Yanying,CHAI Shuhong. Study on a new approach to blade section design and effects of its parameters on cavitation characteristics partⅠ:a new approach to blade section design[J].Journal of Hydrodynamics,1999,14(2):189-200.

[8]王大政,王言英.新型葉剖面設(shè)計(jì)及葉剖面參數(shù)對(duì)空泡特性影響的研究—(2)剖面參數(shù)對(duì)空泡特性的影響[J].水動(dòng)力學(xué)研究與進(jìn)展,2000,15(3):320-328.

WANG Dazheng,WANG Yanying.Study on a new approach to blade section design and effects of its parameters on cavitation characteristics partⅡ:effects of blade section parameters on cavitation characteristics [J].Journal of Hydrodynamics,2000,15(3):320-328.

[9] ZHOU W X,WU Y H,DONG S T.A propeller design method with new blade section for improving cavitation inception under unsteady condition[J].Journal of Ship Mechanics,2001,5(6):18-26.

[10]SU Y M,IKEHATA M,KAI H.A numerical method for designing three-dimensional wing based on surface panel method[J].Journal of the Society of Naval Architects of Japan,1997,182:39-47.

[11]蔡昊鵬,蘇玉民,沈海龍.基于面元法的三維機(jī)翼數(shù)值設(shè)計(jì)[J].海軍工程大學(xué)學(xué)報(bào),2009,21(6):30-35.

CAI Haopeng,SU Yumin,SHEN Hailong.Numerical design of three-dimensional wing based on surface panel method[J].Journal of Naval University of Engineering,2009,21(6):30-35.

[12]譚廷壽,劉飛.非均勻流場(chǎng)中槳葉剖面設(shè)計(jì)方法[J].武漢理工大學(xué)學(xué)報(bào)(交通科學(xué)與工程版),2006,30(3):437-440.

TAN Tingshou,LIU Fei.Design method of blade section in non-uniform flow field[J].Journal of Wuhan University of Technology(Transportation Science& Engineering Edition),2006,30(3):437-440.

[13]李俊華,唐登海,董世湯.指定壓力分布的螺旋槳設(shè)計(jì)[J].船舶力學(xué),2010,14(1/2):10-19.

LI Junhua,TANG Denghai,DONG Shitang.Propeller design by prescribed pressure distribution[J]. Journal of Ship Mechanics,2010,14(1/2):10-19.

[14]黃勝,任萬龍,王超.水翼剖面多目標(biāo)粒子群算法優(yōu)化[J].哈爾濱工程大學(xué)學(xué)報(bào),2014,35(12):1451-1457.

HUANGSheng,RENWanlong,WANGChao. Multi-objective particle swarm optimization of hydrofoil sections[J].Journal of Harbin Engineering University,2014,35(12):1451-1457.

[15]常欣,郭春雨,孟祥印,等.基于粒子群算法的水翼剖面優(yōu)化設(shè)計(jì)[J].船舶工程,2010,32(5):1-3,26.

CHANG Xin,GUO Chunyu,Meng Xiangyin,et al. Design of hydrofoil section based on particle swarm optimization[J].Ship Engineering,2010,32(5):1-3,26.

[16] 曾志波.最大空泡起始航速及最小空泡脈動(dòng)壓力螺旋槳設(shè)計(jì)[D].北京:中國(guó)艦船研究院,2014.

[17]龍文,謝偉,楊向暉,等.基于新型剖面設(shè)計(jì)的螺旋槳多目標(biāo)性能權(quán)衡優(yōu)化[J].艦船科學(xué)技術(shù),2011,33(11):56-60.

LONG Wen,XIE Wei,YANG Xianghui,et al. Multi-objective optimization of open water performance of propeller based on new anti-cavitation airfoil profile design[J].Ship Science and Technology,2011,33(11):56-60.

[18] 饒志強(qiáng),楊晨俊,李巍.基于升力分布的二維翼型設(shè)計(jì)方法[C]//2015年船舶水動(dòng)力學(xué)學(xué)術(shù)會(huì)議.哈爾濱:[s.n.],2015.

[19] 王國(guó)強(qiáng),董世湯.船舶螺旋槳理論與應(yīng)用[M].哈爾濱:哈爾濱工程大學(xué)出版社,2007.

[20] 徐華舫,張炳暄,朱自強(qiáng),等.亞、超音速定常位流的面元法[M].北京:國(guó)防工業(yè)出版社,1981.

[21]KATZ J,PLOTKIN A.Low-speed aerodynamics [M].2nd ed.New York:Cambridge University Press,2001.

Camber lines design approach of propeller blade sections based on the vortex lattice method

RAO Zhiqiang1,2,YANG Chenjun1,2
1 State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China
2 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration,Shanghai 200240,China

Blade sections correlates closely with marine propeller hydrodynamic performance,including loading status,efficiency,cavitation,and other factors,and there is a direct relationship between camber lines and the loading status of blade sections.For the sake of designing camber lines of blade sections with specified loading,a lift distribution is first defined as loading description.A parametric method is then proposed for lift distribution expression,after which the camber line is designed for the specified lift distribution through Newton-Raphson iteration.Finally,the flow characteristics of camber lines are calculated with 2D vortex lattice method.The results indicate that the lift distribution of camber lines designed by combining Newton-Raphson iteration and vortex lattice method matches the given one.Furthermore,certain parameters of the lift distribution are modified to generate camber lines called SJ series for the future design.

propeller blade section;camber line design;lift distribution;vortex lattice method;SJ series camber lines

U664.33

A

10.3969/j.issn.1673-3185.2016.04.006

2015-09-30網(wǎng)絡(luò)出版時(shí)間:2016-7-29 9:45

饒志強(qiáng)(通信作者),男,1986年生,博士生。研究方向:船舶推進(jìn)器水動(dòng)力學(xué)。

E-mail:hellostar@126.com

楊晨俊,男,1964年生,博士,教授,博士生導(dǎo)師。研究方向:船舶推進(jìn)器水動(dòng)力學(xué)。

E-mail:cjyang@sjtu.edu.cn

猜你喜歡
方法設(shè)計(jì)
何為設(shè)計(jì)的守護(hù)之道?
《豐收的喜悅展示設(shè)計(jì)》
流行色(2020年1期)2020-04-28 11:16:38
學(xué)習(xí)方法
瞞天過海——仿生設(shè)計(jì)萌到家
設(shè)計(jì)秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設(shè)計(jì)叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 试看120秒男女啪啪免费| 欧美笫一页| 青青草原国产精品啪啪视频| 亚洲日本精品一区二区| 国产精品99r8在线观看| 国产亚洲精品91| 国产精品一区二区不卡的视频| 91娇喘视频| 国产日韩精品一区在线不卡| 九九热在线视频| 色偷偷av男人的天堂不卡| 国产在线拍偷自揄拍精品| 免费毛片视频| 狠狠做深爱婷婷久久一区| 麻豆精品国产自产在线| 亚洲无码久久久久| 日韩精品免费在线视频| 性网站在线观看| 免费国产黄线在线观看| 亚洲av无码久久无遮挡| 中国丰满人妻无码束缚啪啪| 美女裸体18禁网站| 国产成人免费手机在线观看视频 | 亚洲中文字幕国产av| 中文字幕精品一区二区三区视频 | 黄色网址免费在线| 黄色网在线| 亚洲成a人片77777在线播放| 黄色一级视频欧美| 天天躁夜夜躁狠狠躁图片| 亚洲成人播放| 欧美一区二区三区不卡免费| 亚洲成人高清在线观看| 亚洲美女久久| 99热国产这里只有精品9九 | 国产特级毛片| 久久永久精品免费视频| 九色视频线上播放| 日韩精品一区二区三区swag| 免费观看男人免费桶女人视频| a级毛片网| 999福利激情视频| 亚洲视频色图| 国产91透明丝袜美腿在线| 国产黄网永久免费| 日韩A级毛片一区二区三区| 国产成人AV综合久久| 99视频在线看| 欧美一级大片在线观看| 国产激情第一页| 人妻无码中文字幕第一区| 亚洲国产精品一区二区高清无码久久| 另类专区亚洲| 久久综合亚洲鲁鲁九月天| 国产午夜无码片在线观看网站 | 一区二区欧美日韩高清免费| 国产麻豆91网在线看| 国产一级毛片网站| 91精品日韩人妻无码久久| 久久网综合| 无码专区第一页| 丝袜久久剧情精品国产| a级毛片在线免费| 国产精品丝袜视频| 国产xx在线观看| 影音先锋丝袜制服| 国产在线一二三区| 久久网欧美| 国产青青操| 99草精品视频| 日韩一区精品视频一区二区| 欧美97欧美综合色伦图| 动漫精品啪啪一区二区三区| 国产精品视频第一专区| 中文字幕永久视频| 久久不卡国产精品无码| 国产精欧美一区二区三区| 国产呦精品一区二区三区网站| 国产91高跟丝袜| 久久青草精品一区二区三区| 成年人国产视频| 国产精品美女网站|