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

平面鋼閘門三維有限元線性分析的若干問題

2023-08-28 01:54:12許中武周建方
中國農(nóng)村水利水電 2023年8期
關(guān)鍵詞:有限元模型

許中武,蔡 偉,周建方

(河海大學(xué)機電工程學(xué)院,江蘇 常州 213022)

0 引 言

閘門是用來關(guān)閉、開啟或局部開啟水工建筑物中過水孔口的活動結(jié)構(gòu),其主要作用是控制水位、調(diào)節(jié)流量。作為水工建筑物的重要組成部分,閘門的安全性和適用性在很大程度上影響著整個水工建筑物的運行效果[1],平面鋼閘門作為閘門中使用頻率最高的類型之一,其設(shè)計計算和安全運行也更值得工程設(shè)計人員關(guān)注。對于平面鋼閘門,其設(shè)計計算方法有三維有限元方法和傳統(tǒng)的平面體系計算方法,目前各大設(shè)計院的關(guān)于平面鋼閘門設(shè)計和復(fù)核計算方法仍以平面體系計算方法為主。平面體系計算方法先將平面鋼閘門簡化分解為若干基本構(gòu)件(如面板、主梁、水平次梁、豎直次梁、頂梁和底梁等),再對每一構(gòu)件進行設(shè)計計算[2]。由于平面體系計算方法忽略了各構(gòu)件的協(xié)同作用,且其計算模型和荷載分配均依據(jù)工程經(jīng)驗,可能會造成某些構(gòu)件的材料浪費或局部偏于危險的現(xiàn)象。

隨著三維有限元方法和相關(guān)軟件的成熟和完善,許多專家學(xué)者開始將三維有限元方法應(yīng)用于平面鋼閘門結(jié)構(gòu)計算。同時,美國的《水工鋼結(jié)構(gòu)設(shè)計》[3]、中國的SL 74-2019《水利水電鋼閘門設(shè)計規(guī)范》[4](以下簡稱規(guī)范)和NB 35055-2015《水電工程鋼閘門設(shè)計規(guī)范》[5]等相關(guān)設(shè)計規(guī)范均指出,除了使用傳統(tǒng)的平面體系計算方法外,三維有限元方法也可供平面鋼閘門設(shè)計計算使用。但工程設(shè)計或科研人員在使用三維有限元方法對平面鋼閘門進行分析計算時并沒有相應(yīng)導(dǎo)則可供參考。不同人員對軟件操作和三維有限元方法掌握程度的差異,可能導(dǎo)致對同一閘門的計算結(jié)果不同,從而降低了計算結(jié)果的可信度。為此,亟需平面鋼閘門三維有限元分析導(dǎo)則指導(dǎo)設(shè)計人員進行相關(guān)三維有限元線性分析。目前,已有國內(nèi)學(xué)者對閘門三維有限元分析原則進行了研究。張雪才[3,6]等分別從閘門面板和支臂結(jié)構(gòu)的合理單元類型、單元尺寸或數(shù)量等方面展開分析研究,并得出了相應(yīng)的分析原則。然而,關(guān)于閘門三維有限元分析原則的研究目前仍處于初步階段,并不完善。

本文以某水電站機組進水口快速平面鋼閘門為主要研究對象,結(jié)合相關(guān)設(shè)計資料,基于ANSYS有限元分析軟件,從有限元網(wǎng)格劃分和漏水孔簡化兩個方面系統(tǒng)研究平面鋼閘門三維有限元線性分析時的若干問題,以期相關(guān)結(jié)論可供工程設(shè)計人員在對此類平面鋼閘門進行設(shè)計計算時提供參考,同時也為我國相關(guān)平面鋼閘門設(shè)計規(guī)范的修訂和三維有限元線性分析原則的制定提供依據(jù)。

1 閘門概況

本文研究的水電站機組進水口快速閘門為平面鋼閘門,孔口尺寸為8.8 m×11.0 m(寬×高),底檻高程736.00 m,設(shè)計水位827.83 m,設(shè)計水頭按92.00 m 計算。該平面鋼閘門的門葉結(jié)構(gòu)布置圖如圖1所示,為方便下文闡述相關(guān)計算數(shù)據(jù),將該平面鋼閘門布置的12 根主梁從上到下分別記為1 號到12 號主梁。其中1 號和2 號主梁為箱形結(jié)構(gòu),記為上節(jié)箱形梁;同理11 號和12 號主梁構(gòu)成下節(jié)箱形梁,其余主梁為工形梁;底梁為矩形截面梁,截面尺寸為40 mm×9 960 mm×200 mm。

圖1 門葉結(jié)構(gòu)布置圖Fig.1 Structural drawing of door leaf

2 有限元網(wǎng)格劃分

劃分網(wǎng)格是建立有限元模型的一個重要環(huán)節(jié),所劃分的網(wǎng)格形式對計算精度和計算規(guī)模將產(chǎn)生直接影響[7],因此需要考慮多方面因素。在有限元軟件中網(wǎng)格劃分流程主要包括選定單元屬性、設(shè)置網(wǎng)格尺寸和選擇網(wǎng)格劃分方式等,其中單元屬性包括單元類型、材料類型和單元坐標系等。目前,在平面鋼閘門有限元模型的單元類型選定、網(wǎng)格尺寸設(shè)置和網(wǎng)格劃分方式選擇等方面的相關(guān)依據(jù)和研究較少,科研和工程設(shè)計人員多憑經(jīng)驗劃分網(wǎng)格。本節(jié)將從上述3個方面對平面鋼閘門有限元模型的網(wǎng)格劃分進行討論。

2.1 單元類型

對于平面鋼閘門這樣典型的空間薄壁結(jié)構(gòu),為保證計算精度,在建立其有限元模型時應(yīng)盡量建立完整空間薄壁結(jié)構(gòu)。此時,平面鋼閘門所有構(gòu)件均采用殼單元劃分網(wǎng)格。合理的殼單元類型將決定結(jié)構(gòu)應(yīng)力和撓度的計算精度。文獻[6]以平面鋼閘門面板區(qū)格的幾種計算模型為例,分別驗證了ANSYS軟件中分析薄壁結(jié)構(gòu)常用殼單元的計算精度。然而,文獻[6]僅比較了使用不同單元類型時計算模型撓度的有限元解和理論解,實際上,在平面鋼閘門結(jié)構(gòu)分析時,保證其強度符合規(guī)范要求同樣重要。另一方面,文獻[6]中的計算模型僅討論了平面鋼閘門面板區(qū)格,而實際上平面鋼閘門主梁的簡支梁模型在設(shè)計計算中也很常見。

由于ANSYS 在近幾個版本已將shell43 和shell63 單元從圖形用戶界面操作選項中淘汰,因此本文將選用當前版本中常用的shell181 和shell281 單元。shell181 和shell281 單元均適合對薄的及具有一定厚度的殼體結(jié)構(gòu)進行分析,shell181 單元使用一次插值函數(shù),由4 個結(jié)點組成,每個結(jié)點具有6 個自由度;shell281 單元則使用二次位移插值函數(shù),由8 個節(jié)點組成,各節(jié)點自由度與shell181相同。

平面鋼閘門的平面體系計算方法一般將主梁視為簡支梁模型,將面板區(qū)格視為四邊固定矩形彈性薄板模型,下面將使用這兩種常用計算模型,比較shell181和shell281這兩種單元類型在應(yīng)力和變形上的計算精度。在ANSYS 中分別建立相應(yīng)的計算模型,兩種模型的材料參數(shù)均依據(jù)Q355B 確定。其中,簡支梁截面尺寸為0.3 m×0.05 m(梁高×板厚),支承長度為10 m,均布荷載為10 kN/m;四邊固定矩形彈性薄板尺寸為2 m×1 m×0.05 m(板長×板寬×板厚),均布荷載為0.5 MPa。據(jù)此建立的計算模型如圖2所示。

圖2 有限元計算模型Fig.2 Finite element calculation model

對于受均布荷載的簡支梁,其應(yīng)力和最大撓度可根據(jù)材料力學(xué)公式直接求得。對于受均布荷載的四邊固定矩形彈性薄板,其應(yīng)力計算可依據(jù)規(guī)范和彈性力學(xué)理論[8];在計算其撓度時,根據(jù)彈性力學(xué)理論可知,在長邊與短邊之比為2 時,矩形薄板的最大撓度理論解w為[9]:

式中:q為均布荷載;a為短邊長度;E為彈性模量,取2.1×106MPa;h為薄板厚度。

對受均布荷載的簡支梁和四邊固定矩形彈性薄板模型,在劃分網(wǎng)格建立有限元模型時,分別選用shell181 和shell281 單元,且單元尺寸設(shè)置等均相同,計算后提取相應(yīng)結(jié)果,并與理論解比較,如表1和表2所示。

表1 簡支梁計算結(jié)果Tab.1 Calculation results of simply supported beam

表2 四邊固定矩形彈性薄板計算結(jié)果Tab.2 Calculation results of rectangular elastic thin plate with four edges fixed

由表1 和表2 可知,在荷載、約束和網(wǎng)格尺寸等其他設(shè)置相同的情況下,對于簡支梁和四邊固定矩形彈性薄板模型,shell281 單元無論在應(yīng)力還是變形計算上精度都更高。以簡支梁最大正應(yīng)力為例,使用shell181 單元的有限元解相對于理論解誤差為-20.23%,而shell281 單元有限元解誤差僅為-0.06%。特別地,從四邊固定矩形彈性薄板的計算結(jié)果可知,選用兩種不同的單元類型對撓度計算結(jié)果影響不大,與文獻[6]的結(jié)論一致,說明了上述計算結(jié)果的正確性。但是比較應(yīng)力計算結(jié)果的理論值和有限元解,shell281 單元的精度顯著提高,相差一個量級,計算時間相差不大。因此對平面鋼閘門這樣的空間薄壁結(jié)構(gòu)進行有限元計算時應(yīng)優(yōu)先選取shell281單元。

2.2 網(wǎng)格尺寸

有限元模型網(wǎng)格尺寸對計算精度影響較大,模型網(wǎng)格尺寸的減小在一定程度上會提高計算精度,但當網(wǎng)格尺寸減小到一定程度時,精度提高會很有限,而計算時耗卻會大大增加。此時,可以認為網(wǎng)格尺寸對計算精度的影響可以忽略不計,即有限元計算結(jié)果與網(wǎng)格尺寸之間具有無關(guān)性。在平面鋼閘門有限元計算時,獲取網(wǎng)格無關(guān)性的解是工程設(shè)計或科研人員需要關(guān)注的關(guān)鍵問題。本文將選取shell281 單元,計算分析不同網(wǎng)格尺寸下的平面鋼閘門有限元模型,探討適用于本文所研究的平面鋼閘門的合理網(wǎng)格尺寸范圍。

分別建立7 種不同網(wǎng)格尺寸的平面鋼閘門有限元模型,各模型的網(wǎng)格尺寸如表3 所示。本文在建立有限元模型時,使用ANSYS 默認的全局直角坐標系(笛卡爾全局坐標系),其中x軸沿平面鋼閘門跨度方向,y軸沿平面鋼閘門高度方向,z軸沿水流方向,坐標原點位于平面鋼閘門豎板底部中點處。

表3 模型尺寸分組Tab.3 Model size grouping

該平面鋼閘門采用Q355B 材料,計算時材料彈性模量E取2.06×105MPa,密度ρ取7 850 kg/m3,泊松比v取0.3。平面鋼閘門作用水頭按設(shè)計水頭設(shè)置。此外,1 號主梁腹板承受額外的垂直水柱壓力。根據(jù)平面鋼閘門實際工作狀況,滑塊上受水流方向約束(Uz=0),平面鋼閘門底部受鉛直方向約束(Uy=0)。為保證計算模型的幾何不變性,假定平面鋼閘門底部中點沿跨度方向無位移(Ux=0)。下文所有平面鋼閘門有限元計算的荷載、相關(guān)參數(shù)和約束設(shè)置均相同。

以網(wǎng)格尺寸為0.40 m×0.40 m 的模型1 為例,計算后得到的整體等效應(yīng)力和變形云圖如圖3 所示,提取所有平面鋼閘門模型中主要構(gòu)件的最大等效應(yīng)力和撓度列于表4和表5。

表5 平面鋼閘門主要構(gòu)件最大撓度mTab.5 Maximum deflection of main components of plane steel gate

圖3 模型1有限元計算結(jié)果Fig.3 Finite element calculation results of Model 1

由表4和表5可知:①隨著網(wǎng)格尺寸的減小,平面鋼閘門面板和各主梁的最大等效應(yīng)力呈現(xiàn)增大趨勢,最大等效應(yīng)力出現(xiàn)位置均相同,但網(wǎng)格尺寸對面板最大等效應(yīng)力影響最大。對于面板,從模型1 到模型4,即網(wǎng)格尺寸從0.40 m×0.40 m 到0.25 m×0.25 m 時,面板最大等效應(yīng)力增長平緩;但從模型4 到模型5,即網(wǎng)格尺寸從0.25 m×0.25 m 到0.20 m×0.20 m 時,其最大等效應(yīng)力增長率最大,達到7.56%;從模型5 到模型7,即網(wǎng)格尺寸小于0.20 m×0.20 m 時,面板最大等效應(yīng)力繼續(xù)增長平緩,趨于收斂。此時,可以認為網(wǎng)格尺寸和平面鋼閘門面板最大等效應(yīng)力計算結(jié)果具有無關(guān)性。②網(wǎng)格尺寸對平面鋼閘門主要構(gòu)件撓度影響很小。各模型面板和主梁最大撓度偏差較小。以面板為例,隨著網(wǎng)格尺寸的減小,各模型最大撓度增長率最大僅為0.39%,即在網(wǎng)格尺寸小于0.40 m×0.40 m 時,網(wǎng)格尺寸和平面鋼閘門主要構(gòu)件最大撓度計算結(jié)果即具有網(wǎng)格無關(guān)性。綜上,在對本文所研究的平面鋼閘門進行有限元計算時,其合理網(wǎng)格尺寸應(yīng)小于或等于0.20 m×0.20 m。

2.3 網(wǎng)格劃分方式

對于面,ANSYS 為用戶提供了兩種網(wǎng)格劃分方式:自由網(wǎng)格劃分和映射網(wǎng)格劃分。自由網(wǎng)格劃分通常由有限元計算軟件自動劃分,易出現(xiàn)畸變網(wǎng)格和退化網(wǎng)格[10]。而映射網(wǎng)格劃分則需滿足特定規(guī)則,具有規(guī)則形狀且明顯成排。一般映射劃分得到的網(wǎng)格質(zhì)量優(yōu)于自由劃分得到的網(wǎng)格質(zhì)量,映射劃分有限元模型的計算精度理論上也優(yōu)于自由劃分有限元模型的計算精度。將映射劃分應(yīng)用于平面鋼閘門有限元模型的建立上,需要在前期投入大量工作,對某些復(fù)雜的面需要進行切割,對各個面也只能逐一劃分網(wǎng)格。如果映射劃分方法對平面鋼閘門在應(yīng)力和撓度的計算精度上提高甚微,則會得不償失。

本文使用shell281單元,網(wǎng)格尺寸設(shè)置為0.10 m×0.10 m,分別對該平面鋼閘門進行自由網(wǎng)格劃分和映射網(wǎng)格劃分,生成的局部有限元模型如圖4 所示,其中自由網(wǎng)格劃分得到的有限元模型出現(xiàn)了大量退化網(wǎng)格。兩種網(wǎng)格劃分方式所得有限元模型的單元數(shù)和節(jié)點數(shù)如表6所示。

表6 平面鋼閘門有限元模型節(jié)點和單元數(shù)Tab.6 Node and element number of plane steel gate finite element model

圖4 平面鋼閘門局部有限元模型Fig.4 Local finite element model of plane steel gate

經(jīng)有限元計算后提取平面鋼閘門主要構(gòu)件的最大等效應(yīng)力和撓度列于表7。

表7 平面鋼閘門主要構(gòu)件應(yīng)力和撓度計算結(jié)果Tab.7 The stress and deflection calculation results of the main components of plane steel gate

兩種網(wǎng)格劃分方式下,計算得到的平面鋼閘門整體的應(yīng)力和變形分布均類似,最大等效應(yīng)力和最大變形值偏差較小。由表7可知:與映射劃分模型相比,平面鋼閘門各主梁和面板在自由劃分下的最大等效應(yīng)力值誤差僅為-3.38%,最大撓度誤差不超過1%,即無論是應(yīng)力還是變形,兩種劃分方式下的計算結(jié)果相對誤差都很小。

但實際上,在使用映射劃分方式時,平面鋼閘門有限元模型的單元數(shù)和節(jié)點數(shù)都更多,計算規(guī)模更大,并且映射劃分在前期建立模型階段還需要花費大量時間和精力。此外,雖然通過映射劃分方式得到的網(wǎng)格質(zhì)量更佳,但從計算結(jié)果來看精度提高并不明顯。因此使用有限元方法對此類平面鋼閘門進行計算分析時采用自由劃分方式即可。

3 漏水孔結(jié)構(gòu)簡化

由于平面鋼閘門結(jié)構(gòu)復(fù)雜,構(gòu)件眾多,細節(jié)特征也較多,因此在有限元軟件ANSYS中建立其模型時通常需要適當簡化,使模型與計算機性能和分析目的相匹配。簡化平面鋼閘門模型一方面可以減小有限元求解的規(guī)模,提高計算效率,另一方面也可以降低工程設(shè)計或科研人員在前期的建模工作量。

為了防止平面鋼閘門各主梁上的積水影響啟閉甚至導(dǎo)致其構(gòu)件銹蝕,進而影響平面鋼閘門結(jié)構(gòu)性能,主梁腹板上一般會開設(shè)有漏水孔。如果漏水孔對平面鋼閘門應(yīng)力或變形計算影響不大,則可在建模階段忽略這些漏水孔,減小有限元分析工作量,提高效率。

本文研究的平面鋼閘門在主梁腹板上開設(shè)有若干直徑為200 mm的漏水孔,其中2號至11號各節(jié)主梁腹板分別均勻布置6 個,12 號主梁腹板均勻布置8 個。根據(jù)規(guī)范,一般主梁腹板應(yīng)盡量避免開大孔,漏水孔直徑Ф與主梁截面高度(下文簡稱“梁高”)H之比(下文簡稱“徑高比”)應(yīng)小于0.5,必要時還需對漏水孔采取補強措施,保證結(jié)構(gòu)強度。因此漏水孔的孔徑尺寸是判斷其是否可以簡化的重要參數(shù)。

本文將針對該平面鋼閘門,在保持漏水孔位置不變的情況下,分別建立徑高比為0.1、0.15、0.2、0.25和0.3的平面鋼閘門模型,相應(yīng)的孔徑和梁高參數(shù)設(shè)置如表8所示。同時,還建立了無漏水孔的平面鋼閘門模型,用以與相應(yīng)的計算結(jié)果對比。開展有限元計算并提取平面鋼閘門主梁應(yīng)力和撓度的計算結(jié)果,統(tǒng)計繪制出圖5和圖6。

表8 孔徑和梁高設(shè)置Tab.8 Setting of aperture and beam height

圖5 各模型主梁最大等效應(yīng)力Fig.5 Maximum equivalent stress of main beam for different models

圖6 各模型主梁最大撓度Fig.6 Maximum deflection of main beam for different models

由圖5 和圖6 可知:①在徑高比為0.1 到0.15 時,各模型主梁最大等效應(yīng)力值之間相差很小,出現(xiàn)的位置相同,且各組計算結(jié)果與簡化模型計算結(jié)果基本相同;在徑高比為0.2 到0.3時,各計算組主梁在兩側(cè)漏水孔附近均出現(xiàn)了較為明顯的應(yīng)力集中現(xiàn)象,且奇點應(yīng)力值大于簡化模型主梁最大等效應(yīng)力,并隨徑高比的增大而增大,部分主梁兩側(cè)漏水孔附近奇點應(yīng)力超出相應(yīng)的許用應(yīng)力;②不同徑高比模型的主梁最大撓度均出現(xiàn)在跨中位置,與簡化模型計算結(jié)果發(fā)生的位置相同。隨著徑高比的增大,主梁最大撓度相對偏差逐漸增大。最大偏差出現(xiàn)在徑高比為0.3模型的下節(jié)箱形梁,為10.12%,但各模型主梁撓度都滿足規(guī)范設(shè)計要求。

綜上,當徑高比小于或等于0.15 時,在建立平面鋼閘門有限元模型時可忽略主梁腹板漏水孔,簡化模型;而當徑高比大于0.15 時不可忽略漏水孔,并且需要根據(jù)有限元計算結(jié)果采取相應(yīng)的補強措施。

4 結(jié) 論

為了使工程設(shè)計和科研人員在使用三維有限元方法對平面鋼閘門進行有限元線性分析時有據(jù)可依,且更加科學(xué)高效,本文從網(wǎng)格劃分和模型簡化等方面展開對平面鋼閘門有限元計算過程的若干問題進行了討論研究,主要得到了以下結(jié)論。

(1)本文以平面鋼閘門平面體系計算方法中常用的簡支梁和四邊固定矩形彈性薄板模型為算例,討論了常用殼單元的計算精度。結(jié)果表明,在同樣網(wǎng)格尺寸下,采用shell281單元的應(yīng)力計算精度顯著提高,相差約一個量級,而計算時間并無明顯增加,對平面鋼閘門進行有限元計算時宜選用該單元。

(2)針對本文所研究的平面鋼閘門,在劃分網(wǎng)格建立有限元模型時,合理的網(wǎng)格尺寸設(shè)置應(yīng)小于或等于0.20 m×0.20 m;映射劃分相對于自由劃分得到的網(wǎng)格質(zhì)量更高,但兩種劃分方式下閘門主要構(gòu)件的應(yīng)力和變形計算結(jié)果相對偏差很小,鑒于映射劃分方式前期工作量較大,為提高效率,采用自由劃分即可。

(3)在對平面鋼閘門使用三維有限元方法進行計算分析時,當徑高比小于或等于0.15 時,在建立平面鋼閘門模型時可忽略漏水孔;在徑高比大于0.15 時,由于部分主梁兩側(cè)漏水孔邊緣出現(xiàn)應(yīng)力集中現(xiàn)象,且奇點應(yīng)力大于許用應(yīng)力,此時應(yīng)建立平面鋼閘門完整模型,并需要根據(jù)有限元計算結(jié)果對漏水孔采取相應(yīng)的補強措施。

猜你喜歡
有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
新型有機玻璃在站臺門的應(yīng)用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 在线观看欧美精品二区| 天堂网亚洲系列亚洲系列| 男人天堂伊人网| 国产精品极品美女自在线| 色成人亚洲| 欧美午夜一区| 91网在线| 精品福利国产| 精品一区二区三区四区五区| 成人午夜免费观看| 91国内在线视频| 婷婷亚洲最大| 国产一区免费在线观看| 亚洲狠狠婷婷综合久久久久| 欧美亚洲激情| 亚洲无码91视频| 1024你懂的国产精品| 国产成人1024精品| 欧美精品成人一区二区在线观看| 3D动漫精品啪啪一区二区下载| 国产美女91视频| 毛片网站在线看| 呦女亚洲一区精品| 午夜毛片免费看| 欧美一级专区免费大片| 91精品视频播放| 国产精品综合久久久| 97人人模人人爽人人喊小说| 日韩在线播放欧美字幕| 中文纯内无码H| 成年人国产视频| 亚洲国产精品一区二区第一页免| 日本成人福利视频| a色毛片免费视频| 亚洲精品第一在线观看视频| 狠狠色丁香婷婷综合| 中文字幕有乳无码| 婷婷色中文网| 9久久伊人精品综合| 欧美精品综合视频一区二区| 日本国产精品一区久久久| 国产激情无码一区二区APP| 国产噜噜噜视频在线观看| 白丝美女办公室高潮喷水视频| 成人av专区精品无码国产| 露脸真实国语乱在线观看| 天堂岛国av无码免费无禁网站| 亚洲区视频在线观看| 亚洲AV一二三区无码AV蜜桃| 香蕉eeww99国产在线观看| 人与鲁专区| 色亚洲激情综合精品无码视频| 国产又色又刺激高潮免费看| 精品一區二區久久久久久久網站| 国产日本欧美亚洲精品视| 国产精品香蕉在线| 伊人天堂网| 国产精品刺激对白在线| 成人福利在线观看| 毛片最新网址| 国内精品久久久久久久久久影视 | 伊人久久青草青青综合| 毛片免费视频| 国产成熟女人性满足视频| 天堂av高清一区二区三区| 欧美精品二区| 人妻丰满熟妇αv无码| 欧美成人怡春院在线激情| 亚洲一级毛片在线观播放| 中文字幕免费播放| 免费精品一区二区h| 四虎亚洲国产成人久久精品| 亚洲成人免费在线| 99在线观看国产| 拍国产真实乱人偷精品| 日韩精品毛片人妻AV不卡| 色色中文字幕| AV老司机AV天堂| 国产欧美日韩视频怡春院| 精品一区二区三区视频免费观看| 欧美视频二区| 操美女免费网站|