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

考慮機械載荷和熱載荷的溫室大棚骨架結構輕量化設計

2024-06-17 13:39:03凡健易振峰姚興智謝錦鵬譚文超王昱
中國農機化學報 2024年6期

凡健 易振峰 姚興智 謝錦鵬 譚文超 王昱

摘要:為探究溫度變化對溫室大棚骨架結構安全性的影響,建立考慮機械載荷和熱載荷的大棚骨架熱彈性結構優(yōu)化設計模型。以機械載荷和熱載荷下結構的最大應力值最小化為目標,以結構的總材料用量為約束,考慮機械載荷與熱載荷聯(lián)合作用下的應力分布,實現(xiàn)對大棚骨架的連續(xù)體結構優(yōu)化設計,使結構在滿足支撐剛度的前提下最大程度緩和結構的應力集中問題。考慮到熱應力優(yōu)化問題中的設計依賴性和中間變量問題,采用密度過濾函數獲得清晰的最優(yōu)拓撲構型。通過兩種典型溫室大棚骨架優(yōu)化算例,對比不同溫度變化幅度和不同材料用量下優(yōu)化結構的最大等效應力,結果表明:該模型在相同體積比下結構最大等效應力優(yōu)化效率最高可達到約15%,同種載荷條件下體積比增加0.1可實現(xiàn)結構最大等效應力優(yōu)化效果提高近1%。該研究為考慮機械載荷和熱載荷的大棚骨架熱彈性結構優(yōu)化設計提供有效的設計方法,對工程應用中的溫室大棚骨架結構設計具有指導意義。

關鍵詞:溫室結構;拓撲優(yōu)化;熱彈性結構;應力

中圖分類號:S625.1

文獻標識碼:A

文章編號:2095-5553 (2024) 06-0077-06

收稿日期:2022年10月28日

修回日期:2022年12月15日

*基金項目:國家自然科學基金資助項目(51705161);廣州市基礎與應用基礎研究項目(202102020870)

第一作者:凡健,男,1996年生,安徽滁州人,碩士研究生;研究方向為結構拓撲優(yōu)化、結構設計。E-mail: fanjian@stu.scau.edu.cn

通訊作者:王昱,女,1987年生,長沙人,博士,副教授;研究方向為結構拓撲優(yōu)化、輕量化。E-mail: yu-wang@scau.edu.cn

Lightweight design of greenhouse frame structure considering mechanical load and thermal load

Fan Jian1, 2, Yi Zhenfeng1, 2, Yao Xingzhi1, 2, Xie Jinpeng1, Tan Wenchao1, 2, Wang Yu1, 2

(1. College of Engineering, South China Agricultural University, Guangzhou, 510642, China;

2. Key Laboratory of Key Technology on Agricultural Machinery and Equipment, Ministry of Education,

South China Agricultural University, Guangzhou, 510642, China)

Abstract: In order to explore the influence of temperature change on the safety of greenhouse skeleton structure, an optimal design model of greenhouse skeleton thermoelastic structure considering mechanical load and thermal load was established. With the goal of minimizing the maximum stress value of the structure under mechanical load and thermal load, and choosing the total material consumption of the structure as the constraint, the stress distribution under the combined action of mechanical load and thermal load was considered to achieve the optimal design of the continuum structure of the greenhouse framework, so that the structure could minimize the stress concentration of the structure under the premise of meeting the support stiffness. Considering the design dependency and intermediate variables in the thermal stress optimization problem, the density filter function was used to obtain a clear optimal topology. Through two typical greenhouse frame optimization examples, the maximum equivalent stress of the optimized structure under different temperature changes and material consumption was compared. The results showed that the maximum equivalent stress optimization efficiency of the model could reach about 15% under the same volume fraction. Under the same load conditions, increasing the volume ratio by 0.1 could achieve an increase of nearly 1% in the maximum equivalent stress optimization effect of the structure. The obtained conceptual design scheme of the greenhouse skeleton structure has guiding significance for the design of the greenhouse skeleton structure in engineering applications.

Keywords: greenhouse structure; topology optimization; thermoelastic structure; stress

0 引言

溫室大棚作為一種作物栽培設施,因其建造和運行成本低,在我國北方迅速發(fā)展起來,成為中國設施農業(yè)的主體[1]。按照不同的功能需求,可將常見溫室大棚分為兩大類:連棟溫室和日光溫室。長期以來,溫室結構的有關研究主要集中在溫室結構受到各種機械載荷工況下結構的安全性研究[2-4],忽略了晝夜溫差變化產生的熱應力對大棚結構性能的影響。當熱應力的存在使結構內部應力超過許用應力時,則導致結構破壞,給現(xiàn)實農業(yè)生產帶來安全隱患[5]

除了在典型溫室結構的基礎上進行尺寸優(yōu)化外,拓撲優(yōu)化作為力學工程領域的先進設計方法[6],近年來發(fā)展了一系列的優(yōu)化模型,被應用于解決土木、車輛、儀器設備等[7-9]工程結構的機械載荷和熱載荷作用下的應力問題。桁架結構是溫室大棚骨架結構中最為常見的形式,通過拓撲優(yōu)化方法可合理布置桁架單元的空間分布,可以提升其承載能力[10]。然而,鮮有文獻在針對溫室大棚骨架結構研究的同時,綜合考慮機械載荷和熱載荷的結構應力問題。

本文將基于連續(xù)體拓撲優(yōu)化方法,針對設施農業(yè)對大棚骨架結構的剛度和熱變形以及輕量化的工程需求,建立考慮機械載荷和熱載荷作用下結構應力分布的設計模型,在輕量化的同時使得結構具有較高的剛度和較小的熱變形。為驗證該優(yōu)化模型的可行性,對兩種典型溫室大棚骨架優(yōu)化設計算例,針對常見溫室主要承重結構和工作環(huán)境特點給出優(yōu)化設計方案,獲得滿足實際農業(yè)生產需要的結構。

1 考慮機械載荷和熱載荷的結構優(yōu)化模型

1.1 連續(xù)體結構拓撲優(yōu)化方法

當前連續(xù)體結構優(yōu)化領域的主要方法有:尺寸優(yōu)化、形狀優(yōu)化和拓撲優(yōu)化。拓撲優(yōu)化作為一種啟發(fā)式的結構設計方法,指在滿足已知的約束、載荷和邊界條件的情況下,通過將所求結構的某種力學性能轉化成目標函數,建立其與結構參數的數學模型,進而采用數值優(yōu)化算法優(yōu)計算出最優(yōu)的結構參數,最終獲得最優(yōu)結構的過程。主流的拓撲優(yōu)化方法可分為以均勻化為主的材料插值類方法[11]和描述結構形狀函數類方法,如水平集法[12]。結構的拓撲優(yōu)化旨在設計區(qū)域內尋找一個給定體積的子區(qū)域,使得該區(qū)域對應的目標函數(如結構柔順度、結構位移等)取得極值。引入離散變量的材料密度函數,如式(1)所示。

x=1x∈Ωsolid

0x∈Ωvoid(1)

式中: Ωsolid——實體材料域;

Ωvoid——孔洞材料域。

由于整數模型的計算求解非常困難,通常采用變量連續(xù)化方法將0~1整數變量問題變?yōu)?、1間的連續(xù)變量優(yōu)化模型。上述變量連續(xù)化后的模型是一個病態(tài)問題,優(yōu)化解表現(xiàn)為中間密度值、棋盤格和網格依賴性,可采用中間密度懲罰模型和棋盤格控制措施解決上述問題。常用的中間材料懲罰模型有SIMP(Solid Isotropic Material with Penalization Model)[13]和RAMP(Rational Approximation of Material Properties)[14]模型。由于熱載荷問題的復雜性,本文的拓撲優(yōu)化模型基于RAMP模型建立,將在材料插值模型部分詳細闡述。

1.2 熱彈性結構優(yōu)化模型

考慮機械載荷和熱載荷的的熱彈性結構輕量化設計問題,結構變形為機械載荷與熱載荷聯(lián)合作用下的總體變形,以結構的最大應力最小化作為目標函數,以結構的體積比(即總材料用量)為約束,建立相應的熱彈性結構優(yōu)化模型如式(2)所示,示意圖見圖1。

本文模型采用優(yōu)化領域廣泛應用的漸進優(yōu)化法MMA(Method of moving asymptotes)[15]進行求解。

minσmaxpn

s.t.KU=Fm+Fth

V=∑Ni=1xiv0≤V0

0min≤xi≤1 i=1,2,3,…,N(2)

式中: σmax——結構最大應力;

σpn——結構等效最大應力;

K——總剛度矩陣;

U——總節(jié)點位移矩陣;

Fm——作用在結構上的機械載荷;

Fth——變化溫度場產生的熱載荷;

V——結構優(yōu)化后的體積;

xi——第i個單元的單元密度;

v0——單元的體積;

V0——給定材料體積上限;

xmin——單元的最小密度值,取0.001以避免總剛度矩陣的奇異;

N——設計變量數目。

2 優(yōu)化模型的實施

2.1 應力計算

針對基于密度的最大應力最小化結構優(yōu)化問題,拓撲優(yōu)化迭代過程中出現(xiàn)的奇異現(xiàn)象存在于基于密度的拓撲優(yōu)化中,即當單元密度為很低的值時依然存在應變導致局部出現(xiàn)人造大應力的現(xiàn)象。應力松弛方法[16]以減輕低密度單元的應力,可使優(yōu)化穩(wěn)定進行。同時,在衡量結構的應力時,只能得到每個單元的局部應力大小,無法建立結構整體應力水平和局部單元應力大小的聯(lián)系,但考察每一個單元則大大加重了計算成本。因此,使用P范數近似應力最大值[17],將各個單元應力分量整合成一個衡量值,提高計算效率。

對于熱彈性結構而言,第i個單元的應力矢量σi如式(3)所示。

σi=D0(BUi0ΔT?T)(3)

式中: D0——固體材料的彈性矩陣;

B——應變—位移矩陣;

α0——固體材料的熱膨脹系數;

ΔT——結構受到的均勻溫升場;

?T——列向量,?=[1 1 0]。

對于二維平面問題,單元應力σi是一個3×1的列向量,不便于優(yōu)化的開展。采用von·mises應力準則,將σi整合成一個等效應力σvmi

σvmi=(σi12i22+3σi32-σi1σi212(4)

式中: σi1、σi2、σi3——σi的三個分量,表示為水平方向應力、豎直方向應力和切應力。

為解決應力奇異問題,這里引入q松弛法對單元等效應力進行懲罰

σvmi=xqiσvmi(5)

式中: q——懲罰因子,這里取0.5。

為了實現(xiàn)目標函數連續(xù)可導,采用P范數近似全局最大應力

σpn=∑Ni=1σvmiP1P(6)

式中: P——范數系數,本文P取8。

2.2 材料插值

熱彈性結構相對傳統(tǒng)只受機械載荷的優(yōu)化問題更加復雜,因其敏度與設計變量有關,在優(yōu)化過程中存在符號改變,給優(yōu)化過程帶來不穩(wěn)定性。傳統(tǒng)的SIMP插值方法在表示單元偽密度的設計變量x趨近于0時,其導數為0,這阻礙了熱應力敏度符號的轉變,影響優(yōu)化過程的順利進行。因此,本文引入RAMP插值方法分別懲罰剛度和熱應力。

對于第i個單元,其單元剛度矩陣和節(jié)點熱應力矢量表達如式(7)所示。

Ki=∫ΩBiTDiBidΩ(7)

Fthi=∫ΩBiTDiεthidΩ(8)

式中: Di——密度為xi的單元等效彈性矩陣;

εthi——單元i的熱應變矩陣。

Di=E(xi)D0(9)

εthi=α(xi)ΔT?T(10)

式中: E(xi)——材料的插值函數;

α(xi)——單元密度為xi時的等效熱膨脹系數。

因此,式(8)可被改寫為

Fthi=E(xi)α(xi)ΔT∫ΩBiTDi?TdΩ(11)

其中只有αi和Di與變量xi有關,引入表征材料固有特性的熱應力系數(TSC)[18]概念來表達三者關系。

β(xi)=E(xi)α(xi)(12)

根據RAMP插值方法,等效彈性矩陣Di和熱應力系數β(xi)可以表達成

E(xi)=xi1+SE(1-xi)E0(13)

β(xi)=xi1+Sβ(1-xi)E0α0(14)

式中: SE——楊氏模量插值系數;

Sβ——TSC插值系數;

E0——固體材料的楊氏模量;

α0——固體材料的熱膨脹系數。

SE和Sβ值的選取影響迭代的穩(wěn)定性,并且當SE<Sβ時迭代產生震蕩[19]。根據數值試驗經驗,這里選取SE=8,Sβ=0。

2.3 密度過濾

基于變密度法的拓撲優(yōu)化方法是通過將0-1離散變量轉化成連續(xù)變量,這將導致在優(yōu)化過程中產生無物理意義的具有中間密度的單元。選用合理的材料插值模型參數,可懲罰中間密度往0-1方向逼近。在產生數學層面上局部最優(yōu)解的同時,結構的局部區(qū)域伴隨產生材料密度為0和1周期分布的狀態(tài),即棋盤格現(xiàn)象,增加了結構加工成型的難度。針對上述問題以及2.1節(jié)提到的應力非線性問題,單元的應力極易受其相鄰元素的影響。通過施加密度過濾,可有效防止棋盤格現(xiàn)象并且降低相鄰元素對于單元應力的影響,使優(yōu)化收斂于全局解。

本研究中采用的密度過濾可定義為

x~(k,l)=∑m,nH(m,n)x(k-m,l-n)∑m,nH(m,n)(15)

式中: x(i,j)——第i行第j列的單元密度;

x~(i,j)——第i行第j列單元過濾后的偽密度;

H(m,n)——離散變量m,n的過濾核函數。

H(m,n)=max[0,rmin-δ(m,n)](16)

式中: δ(m,n)——相距m行n列兩單元的中心距離;

rmin——最小過濾半徑。

rmin的大小影響著過濾后結構細節(jié)枝干的數量,同時決定結構的最小尺寸。

3 典型大棚骨架結構輕量化設計

針對實際農業(yè)應用領域兩種典型溫室結構:連棟溫室和日光溫室,如圖2所示,分別選取其主要承重結構作為設計域進行優(yōu)化。建立兩種溫室主體框架部分二維模型,以均勻密度場作為初始設計,設定外圍輪廓作為非設計域,代入本文的模型進行優(yōu)化。選擇現(xiàn)實生活中常用的Q235鋼材作為材料,其楊氏模量E為200GPa,泊松比μ為0.3,熱膨脹系數α為1.20×10-5,許用應力為235MPa。

3.1 連棟溫室主體承重結構優(yōu)化

連棟溫室主體結構部分如圖3所示。該部分為軸對稱結構,為減少優(yōu)化過程的計算量,截取框架右半邊作為設計區(qū)域。結構長l=800cm,高h=400cm,厚度t取1cm,非設計域(圖3中黑色邊框)寬度為20cm,整個區(qū)域離散成80×40的平面四節(jié)點有限單元。結構同時受到均勻溫度場ΔT和分布載荷t-作用,分布載荷隨著外輪廓傾角的變化呈1000N/m到10N/m漸變。右下方采用固定約束,選擇0.3作為體積比。為探究不同溫度變化對結構應力的影響,ΔT分別選擇0℃、10℃、20℃和50℃,優(yōu)化后拓撲結構和應力分布見圖4,優(yōu)化結果見表1(表中表達式上標init和opt分別表示初始構型和優(yōu)化構型),η表示優(yōu)化效率,計算如式(17)所示。

η=σinitpn-σoptpnσinitpn(17)

由表1可知,當連棟溫室結構只受到機械載荷時,優(yōu)化后的結構與實際農業(yè)生產應用的溫室結構相似,其應力水平遠低于材料的許用應力。加入熱載荷后,結構的最大等效應力根據溫升大小的變化同幅度的增加,伴隨著優(yōu)化效率的小幅提升。從圖4可知,等效最大應力值大于結構實際應力最大值,對于P范數等效最大應力的方法造成的與真實應力的差距,可通過增大范數系數P的值縮小。

此外,升溫幅度的不同,結構拓撲也有相應的改變,但始終保持三條棱的樣式。溫升為10℃時,優(yōu)化后的結構應力水平在許用應力之下;20℃及以上的溫升,結構的最大應力則嚴重超過了許用應力。研究發(fā)現(xiàn),優(yōu)化后的結構最大應力出現(xiàn)在右下角支撐區(qū)域處,通過對初始框架模型右下角進行局部優(yōu)化后再進行拓撲優(yōu)化,可有效降低結構應力集中。

3.2 日光溫室主體承重結構優(yōu)化

日光溫室主體框架部分如圖5所示。結構長l=800cm,高h=400cm,厚度t取1cm,非設計域寬度為20cm。整個區(qū)域同樣離散成80×40的平面四節(jié)點有限單元。左側斜坡受到均布載荷800N/m,右側斜坡受到均布載荷t-為100N/m。為探究不同體積比對優(yōu)化結果的影響,設定溫升ΔT為10℃不變,分別選取體積比等于0.2、0.3、0.4,對應優(yōu)化后的拓撲結構和應力分布見圖6,結果如表2所示。

通過表2可知,在同一溫升下,體積比增加0.1,優(yōu)化后結構的最大等效應力減小了約1%且幅度逐漸變小。因此,在實際農業(yè)生產應用中進行優(yōu)化設計時,應當權衡好應力大小和材料用量的關系,選擇最合適的參數組合。

4 結論

本文建立考慮機械載荷和熱載荷的大棚骨架熱彈性結構優(yōu)化設計模型,將其應用于溫室大棚骨架的輕量化設計,同時探究熱載荷對溫室大棚骨架的結構應力分布的影響。通過探究不同溫升和不同材料用量對優(yōu)化后結構應力狀況的影響。

1) 熱載荷的加入大幅提高結構的最大應力:在10℃溫升時,結構最大應力小于材料許用應力,設計滿足使用要求;溫升在20℃及以上時,最大應力超出許用應力,最大應力位置不變指導設計需要對結構初始進行改良。

2) 同等溫升條件下,總材料用量增加0.1,結構的最大等效應力減小約1%。因此,在設計時應結合實際條件需要,選擇出最佳的結構應力大小和總的材料用量組合。

本文不僅為考慮機械載荷和熱載荷的大棚骨架熱彈性結構優(yōu)化設計提供有效的設計方法,還給出具有指導意義的符合工程需求的大棚骨架結構概念設計方案。

參 考 文 獻

[1]劉志杰, 鄭文剛, 胡清華, 等. 中國日光溫室結構優(yōu)化研究現(xiàn)狀及發(fā)展趨勢[J]. 中國農學通報, 2007(2): 449-453.

Liu Zhijie, Zheng Wengang, Hu Qinghua, et al. Current situation and development on structure optimization of solar greenhouse in China [J]. Chinese Agricultural Science Bulletin, 2007(2): 449-453.

[2]肖林剛, 鄒志榮, 吳樂天, 等. 寒冷干旱地區(qū)日光溫室結構的優(yōu)化設計[J]. 農學學報, 2014, 4(5): 52-55.

Xiao Lingang, Zou Zhirong, Wu Letian, et al. Structural optimization and design of the sunlight greenhouse in cold and arid Regions [J]. Journal of Agriculture, 2014, 4(5): 52-55.

[3]劉麗霞, 雷法究. 兩種典型日光溫室結構安全性能分析[J]. 林業(yè)機械與木工設備, 2016, 44(11): 18-23.

Liu Lixia, Lei Fajiu. Safety performance analysis of two typical greenhouse structure [J]. Forestry Machinery & Woodworking Equipment, 2016, 44(11): 18-23.

[4]Briassoulis D, Dougka G, Dimakogianni D, et al. Analysis of the collapse of a greenhouse with vaulted roof [J]. Biosystems Engineering, 2016, 151: 495-509.

[5]李世萍, 秦小楠, 馬倩, 等. 北京力學會第18屆學術年會論文集: 工程應用: 輕鋼溫室結構溫度應力及其對策[C]. 北京: 北京力學會, 2012.

Li Shiping, Qin Xiaonan, Ma Qian, et al. Proceedings of the 18th annual academic conference of beijing mechanics society: Engineering application: Temperature stress of light steel greenhouse structure and its countermeasures [C]. Beijing: Beijing Society of Theoretical and Applied Mechanics, 2012.

[6]Bends?e M, Sigmund O. Topology optimization: Theory, method and applications [M]. New York: Springer-Verlag Berlin Heideberg, 2003.

[7]He Y, Cai K, Zhao Z L, et al. Stochastic approaches to generating diverse and competitive structural designs in topology optimization [J]. Finite Elements in Analysis and Design, 2020, 173: 103399.

[8]Szepessy Z, Zoltan I. Thermal dynamic model of precision wire-wound resistors [J]. IEEE Transactions on Instrumentation and Measurement, 2002, 51(5): 930-934.

[9]Mankame N D, Ananthasuresh G K. Topology synthesis of electrothermal compliant mechanisms using line elements [J]. Structural and Multidisciplinary Optimization, 2004, 26(3): 209-218.

[10]Xu L, Min H. Optimum design of cold formed steel residential roof trusses [J]. International Specialty Conference on Cold-Formed Steel Structures, 2000.

[11]Suzuki K, Kikuchi N. A homogenization method for shape and topology optimization [J]. Computer Methods in Applied Mechanics and Engineering, 1991, 93(3): 291-318.

[12]Wang M Y, Wang X, Guo D. A level set method for structural topology optimization [J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1): 227-246.

[13]Rietz A. Sufficiency of a finite exponent in SIMP (power law) methods [J]. Structural and Multidisciplinary Optimization, 2001, 21(2): 159-163.

[14]Stolpe M, Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization [J]. Structural and Multidisciplinary Optimization, 2001, 22(2): 116-124.

[15]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.

[16]Cheng G D, Guo X. ε-relaxed approach in structural topology optimization [J]. Structural optimization, 1997, 13(4): 258-266.

[17]Duysinx P, Bends?e M P. Topology optimization of continuum structures with local stress constraints [J]. International Journal for Numerical Methods in Engineering, 1998, 43(8): 1453-1478.

[18]Deaton J D, Grandhi R V. Stiffening of restrained thermal structures via topology optimization [J]. Structural and Multidisciplinary Optimization, 2013, 48(4): 731-745.

[19]Gao T, Zhang W. Topology optimization involving thermo-elastic stress loads [J]. Structural and Multidisciplinary Optimization, 2010, 42(5): 725-738.

主站蜘蛛池模板: 在线观看免费国产| 成年网址网站在线观看| 日韩高清一区 | 青青草原国产免费av观看| 午夜国产不卡在线观看视频| 72种姿势欧美久久久久大黄蕉| 亚洲不卡网| 国产成人调教在线视频| 欧美激情福利| 亚洲欧美日韩中文字幕一区二区三区 | 国产手机在线ΑⅤ片无码观看| 国产精品自拍露脸视频| 一级爆乳无码av| 成人午夜网址| 91国内视频在线观看| 波多野结衣在线se| 中文成人无码国产亚洲| 第一区免费在线观看| 国产99在线观看| 欧美综合中文字幕久久| 国产成人综合日韩精品无码不卡| 青青网在线国产| 国产成人福利在线视老湿机| 亚洲国产成人久久77| 亚洲国产中文欧美在线人成大黄瓜| 久久婷婷人人澡人人爱91| 久久黄色一级视频| 国产91全国探花系列在线播放| 97国产精品视频自在拍| 国产精品亚洲αv天堂无码| 国产精品香蕉在线| 97国产在线播放| 老司机久久99久久精品播放 | 欧美第九页| 中文字幕永久视频| 亚洲国产天堂久久综合| 91精品免费高清在线| 日本一区二区三区精品国产| 国产aⅴ无码专区亚洲av综合网| 国产午夜福利片在线观看| 亚洲有码在线播放| 亚洲码一区二区三区| 国产成人精品一区二区| 亚洲国产天堂久久九九九| 日本a∨在线观看| 波多野结衣中文字幕一区二区 | 99re精彩视频| 97视频精品全国在线观看| 日韩无码精品人妻| 久久午夜夜伦鲁鲁片不卡| 亚洲国产综合精品一区| 国产精品亚洲欧美日韩久久| 日韩欧美视频第一区在线观看| 91福利一区二区三区| 91青青草视频| 在线精品亚洲国产| 亚洲91精品视频| 国产无遮挡猛进猛出免费软件| 日本黄色a视频| 国产午夜精品鲁丝片| 精品国产福利在线| 国产自在线播放| 国产a网站| 国产精品私拍在线爆乳| 日韩在线网址| 国产自在线播放| 亚洲福利一区二区三区| 亚洲精品无码久久毛片波多野吉| 69综合网| 99热这里只有成人精品国产| 老司机精品99在线播放| 国产后式a一视频| 国产欧美日韩资源在线观看| 国产成人综合日韩精品无码首页 | 伊人欧美在线| 国产亚洲欧美日韩在线一区二区三区| 五月婷婷亚洲综合| 免费全部高H视频无码无遮掩| 日韩一区精品视频一区二区| 精品精品国产高清A毛片| 高清久久精品亚洲日韩Av| 奇米影视狠狠精品7777|