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

基于非線性規劃和有限元法的懸臂式桁架優化設計

2020-06-11 00:34:04馬佳吳世寶張瑞媛高婧
河北工業科技 2020年3期
關鍵詞:有限元法

馬佳 吳世寶 張瑞媛 高婧

摘要:為了對懸臂式桁架結構進行改進與優化,提出了一種更為簡單可行的非線性規劃模型。首先,在充分分析懸臂式桁架使用要求的基礎上,創新性地將懸臂端節點撓度和結構的可靠性分別視為目標函數和限制條件,并通過靜力學分析確定了目標函數和限制條件的具體表達式。然后,利用Python語言編程對優化方案進行求解,發現在給定的限制條件下,當桁架的總長為4 m、總高為26 m且桿件的橫截面積為1 400 mm2時,懸臂端節點的撓度最小。最后,利用ABAQUS軟件并根據有限元法來判斷優化方案的合理性。有限元法的計算結果表明:得出的優化方案是較為準確合理的,且適當減少桁架的總長并增大桁架的總高和桿件的橫截面積,可以降低懸臂端節點的撓度。因此,基于非線性規劃模型提出的目標函數和限制條件為其他類型桁架結構的優化設計提供了新思路。

關鍵詞:非線性規劃;有限元法;懸臂式桁架; Python;ABAQUS

中圖分類號:O342文獻標識碼:A

doi: 10.7535/hbgykj.2020yx03011

Abstract:

In order to improve and optimize the structure of the cantilever truss, a more simple and feasible nonlinear programming model was proposed. First of all, based on fully analyzing the requirements of cantilever truss, the deflection of the cantilever truss end node and the reliability of the structure were innovatively regarded as the objective function and the constraint condition respectively, and the static analysis was carried out to figure out the specific expressions of the objective function and the constraint condition. Subsequently, the optimization scheme was solved by Python language programming, which revealed that with the total length of this truss of 4 m, the total height of

2.6 m and the crosssectional area of 1 400 mm2, the deflection of the cantilever truss end node was minimum under the given constraints. Finally, ABAQUS software was utilized to judge the rationality of the optimization scheme based on finite element method. The calculation results of finite element method demonstrate that the optimization scheme is correct and reasonable to some extent. In addition, decreasing the total length of truss and increasing the total height of truss and the crosssectional area of the bars appropriately can reduce the deflection of this cantilever truss end node. Therefore, the objective function and the constraint condition of the nonlinear programming model can provide a new idea for presenting better optimization design model of other types of truss structures.

Keywords:

nonlinear programming; finite element method; cantilever truss; Python; ABAQUS

作為一種格構化的梁式結構[1],桁架結構以其優越的承載能力以及抗拉(壓)性能,被廣泛應用于建筑、煤礦、機械、航空、船舶等諸多工業領域的大跨度主承力結構中。因而,研究這類大跨度桁架結構的力學特性、結構強度、仿真分析和施工技術對充分發揮此類結構的承載能力和保障結構安全具有十分重要的意義。例如,張國偉等[2]利用靜力試驗和ABAQUS軟件對影響M型鋼桁架組合樓板的保溫與受力性能的參數進行了研究,發現增加型鋼厚度和鋼筋直徑可以有效地提高樓板的承載能力;袁國平等[3]對大型煤場中廣泛使用的大跨度拱形拉索桁架結構的關鍵施工技術進行了總結,指出在施工時應主要對跨中節點進行監測,且應適當增加拉索的實際內力;廖能解等[4]基于SolidWorksWorkbench平臺對某種桁架機器人的2種典型工況進行了靜力學分析和模態分析,指出當伺服電機的激振頻率在40~50 Hz時,不會影響桁架機器人的工作精度以及穩定性;李元章等[5]根據積木式驗證方法對某飛艇的桁架式復合材料龍骨結構進行了強度驗證,認為根據此方法得出的強度結論滿足工程應用需要,且結論偏向保守;崔宏林等[6]采用直接計算法對某艘鋼管桁架式連接橋結構的小型雙體船的連接橋部位進行了局部強度評估,認為在5級海況下,橫向支撐桿件和強、弱框架都無法滿足規定的強度要求。

正如文獻[6]所指出的,一旦桁架結構不再滿足規定的使用要求,則必須對其進行改進與優化,以重新達到目標要求。因而,研究桁架結構的優化模型或優化算法對獲得更為合適的優化結果有著重要的指導價值。例如,喬心州等[7]基于區間變量描述不確定參數提出了一種具有非概率可靠性的桁架結構形狀優化模型,并發現此模型可以快速穩定地得到收斂的最優解,且符合工程設計經驗;郝寶新等[8]改進了傳統的桁架結構拓撲優化模型,并提出了一種綜合考慮多種約束的桁架結構拓撲優化模型——半定規劃模型,此模型有效地解決了在應用傳統拓撲優化模型時出現的優化模型非凸以及多重特征值不存在常規梯度等問題;李彥蒼等[9]通過在傳統的海豚群算法中引入信息熵的方法來減少盲目搜索,并發現改進后的海豚群算法在克服了傳統海豚群算法容易陷入局部最優解問題的同時,表現出更好的收斂速度和尋優精度;李沛豪等[10]通過在粒子群算法中設置更為合理的罰函數機制來對桁架結構截面進行優化設計,并指出與傳統的粒子群算法相比,改進后的粒子群算法穩定性更好,尋優效率更高。

基于上述分析,為了確定某個懸臂式桁架的優化方案,以懸臂端節點撓度最小為目標函數,利用一個非線性規劃數學模型,在控制桿件的尺寸并保證桿件不發生破壞的前提下,結合Python語言編程,求得懸臂端節點的撓度最小時桁架結構中桿件的尺寸。隨后,利用Python計算出的優化尺寸,在ABAQUS軟件中建立此桁架的力學模型并根據有限元法計算出懸臂端節點的撓度、鉸支座支反力和桿件應力,通過對比Python計算出的最小撓度和ABAQUS軟件計算出的最小撓度以及觀察ABAQUS軟件計算出的鉸支座支反力和桿件應力是否滿足相應的約束條件,來判斷優化模型是否合理。

1優化模型的建立

在對懸臂式桁架結構進行優化設計時,一般需要結合實際情況考慮以下幾個方面的問題:1)由于工作空間有限,需要對桁架結構中各個桿件的尺寸(如長度、橫截面積等)進行必要的限制;2)由于懸臂式結構一般都是用鉸支座固連到墻壁上,因而在桁架承受一定的外載荷后,鉸支座處會產生相應的約束力,其大小不能超過規定的許用值;3)由于桁架結構中的桿件在受到外載荷后都會產生拉伸或壓縮變形,其拉應力或壓應力不能超過材料的許用拉應力或壓應力;4)對于懸臂式結構,其最大的變形往往發生在懸臂端,因而在優化懸臂式桁架時,也要求其懸臂端節點的撓度盡可能小。由此可見,這一優化模型的目標函數和限制條件不都是線性函數,其中也包含許多非線性函數。因而,可以基于非線性規劃的原理[11]來建立懸臂式桁架結構的優化模型。

1.1構建目標函數

由上述分析可知,懸臂端節點的撓度最小,可以作為此優化模型的目標函數。在已知懸臂式桁架結構所受的外載荷時,可以使用單位載荷法(莫爾積分)來計算懸臂端節點的撓度[12]。設此桁架結構中共有

n個長度不完全相等的桿件,它們的彈性模量均為E,橫截面都是面積為A的圓,由單位載荷法可得懸臂端節點的撓度y的表達式為

1.2.2限制鉸支座所受的約束力

為了保證結構的安全,需要將鉸支座所受的約束力限制在一定的范圍內。設此桁架結構與墻壁之間共有m個鉸支座,且各個鉸支座所受的總約束力為Fj,其大小可以通過靜力學分析求出。若其許用值為[F],則此約束條件可表示為

2算例分析

2.1工程概況

某懸臂式桁架簡化后的桁架結構及其受到的外載荷如圖1所示。此結構中的桿件均由Q235鋼制成,其彈性模量為211 GPa,屈服極限為235 MPa,A為橫截面積;A—J為鉸接的節點;K和L是鉸支座,其中鉸支座K是固定鉸支座,鉸支座L是可動鉸支座,且水平方向2個相鄰節點之間的距離均為a,豎直方向2個相鄰節點之間的距離為b/2。節點A,C,F,I均受到大小為P=35 kN、方向豎直向下的載荷,節點B,E,H均受到大小為Q=120 kN、方向豎直向上的載荷。

由于工作條件發生了變化,圖1所示的懸臂式桁架原有的結構不再滿足新的使用要求,需要對其進行改進與優化。由于空間有限,在改造該懸臂式桁架時,對其結構尺寸提出了一些新的要求,包括其總長4a不少于4 m但不超過8 m,其總高b不少于1 m但不超過3 m,每個桿件的橫截面積A不少于700 mm2但不超過1 400 mm2,且K和L兩處鉸支座的約束力均不得超過300 kN,并取安全系數S為15。在滿足上述要求的同時,還應滿足在給定載荷的作用下,使節點A的撓度盡可能小。

2.2將目標函數具體化

要想確定具體的目標函數和限制條件,首先需要對圖1所示的桁架結構進行靜力學分析,以求出K和L兩處的約束力以及各個桿件的軸力。取此桁架結構整體為研究對象,解除鉸支座K和L處的約束,由于K是固定鉸支座,解除該處的約束后可以用水平和豎直2個相互正交的約束力代替;L是可動鉸支座,因而只有1個水平方向的約束力[13]。桁架所受的各個力形成了平面任意力系,因而可以列出3個獨立的平衡方程,并由此解得鉸支座K和L兩處的約束力,其中,K處的水平約束力為FKX=(6Qa-10Pa)/b,豎直約束力為FKY=(3Q-4P);L處的約束力為FL=(6Qa-10Pa)/b。

然后利用節點法,分別求出此桁架結構中各個桿件在原始載荷作用下的軸力,計算結果如表1所示。

由表1可知,在原始載荷的作用下,桁架結構中的各個桿件所受軸力的方向并不相同。由工程概況中給出的載荷大小,可以判斷出此桁架結構中的各個節點都有向上運動的趨勢,因而位于該結構上方的桿件一般都承受壓縮載荷,結構下方的桿件一般都承受拉伸載荷,且越靠近鉸支座的桿件,所受的軸力越大。

去掉該結構所受的所有載荷,并在節點A施加1個方向豎直向上,大小為1的單位載荷R,利用類似的方法,可以得到在此單位載荷R的作用下,此結構中各個桿件的軸力,計算結果如表2所示。

由表2可知,由于施加的單位載荷R的方向是豎直向上,因而該桁架結構總體的變形趨勢與原始載荷作用下該桁架結構的變形趨勢大致相同。除此之外,在單位載荷R的作用下,該桁架結構中各個桿件軸力的方向也與原始載荷作用下該桁架結構中各個桿件軸力的方向不完全相同。事實上,這是一種正常現象,在計算時只需要按照單位載荷法規定的符號要求代入莫爾積分的表達式即可。

由表3可知,當水平方向相鄰2節點之間的距離a=1 m,鉸支座K和L之間的距離b=2.6 m且桿件的橫截面積A=1 400 mm2時,懸臂端節點A的撓度取得最小值為3.473 mm。由一般的工程經驗可知,對于一個懸臂式結構,外伸的長度越短,其懸臂端的撓度越小。除此之外,對于等直桿一類的構件,橫截面積越大,其抗拉(壓)剛度越大,越有利于抵抗外力產生的變形。因而,可以初步判斷表3所示的優化結果是合理的。

3優化模型的合理性驗證

為了判斷上述用于優化此懸臂式桁架的非線性規劃模型的合理性,本文利用ABAQUS軟件對該模型的優化結果進行驗證。具體的驗證思路:1)利用表3所示的已經優化好的該懸臂式桁架中各個桿件的尺寸(水平方向相鄰2節點之間的距離a、鉸支座K和L之間的距離b、桿件的橫截面積A),在ABAQUS軟件中建立此懸臂式桁架結構的有限元模型;2)給建立好的有限元模型施加工程概況中所給出的各個載荷;3)通過計算在給定載荷下ABAQUS軟件利用有限元法[1415]得出的節點A的撓度與表3中所示的節點A的撓度之間的相對誤差,同時觀察ABAQUS軟件根據有限元法計算出的鉸支座K和L處的支反力以及各個桿件中的最大應力σmax是否滿足2.3節中的各項約束條件,來判斷此非線性規劃模型是否合理。

3.1有限元模型的建立過程

3.1.1創建幾何模型并劃分網格

在“部件”模塊的“創建部件”對話框中,設置模型空間為“二維平面”、類型為“可變形”、基本特征為“線”,即可開始繪制此懸臂式桁架結構的草圖。根據表3中給出的桿件尺寸,可以確定圖1所示的懸臂式桁架結構中各個節點在二維平面內的坐標,分別為A(4 000 mm, 0),

然后給圖2所示的幾何模型劃分網格。由于桁架結構中的各個桿件在承受外載荷時都只受軸向的拉力或壓力,因而在劃分網格時,可以將每個桿件視為1個網格單元,而桿件內部不必再細化網格。ABAQUS軟件中默認桁架結構的網格類型為T2D2,這是一種專門為二維二節點桁架劃分網格的網格單元類型。由于圖2中的結構共有21個桿件,因而共需要劃分21個網格單元。

3.1.2定義材料和截面屬性

首先給各個桿件賦予相應的材料。由于此懸臂式桁架結構中的各個桿件在承受外載荷時發生的變形都是彈性變形,且沿各個方向的力學行為相同,因而需要在“屬性”模塊中創建一種各向同性材料,即用于制造各個桿件的Q235鋼,其楊氏模量為211 000 MPa,泊松比為0.3。給各個桿件創建相應的截面,主要的創建步驟如圖3所示。

由圖3 a)可知,在設置桿件的截面類型時,可以選擇“梁”類別中的“桁架”類型,且在ABAQUS軟件中,“桁架”類型的截面形狀默認為圓形。由圖3 b)可知,桁架中桿件的橫截面積被設置為1 400 mm2,與表3中給出的橫截面積相同。

3.1.3設置邊界條件和載荷

首先需要在初始分析中為鉸支座K和L施加邊界條件。由于鉸支座K是固定鉸支座,需要限制其在X和Y方向上的位移;鉸支座L是可動鉸支座,只需要限制其在X方向的位移。隨后,在節點A,B,C,E,F,H,I處按工程概況施加相應的載荷,即可完成載荷和邊界條件的設置。所有定義完成的載荷和邊界條件如圖4所示。

由圖4可知,鉸支座K處出現了相互正交的2個小三角形,表示限制了X和Y方向上的位移;鉸支座L處出現了1個水平方向的小三角形,說明對X方向的位移進行了限制。此外,節點A,C,F,I處出現了豎直向下的箭頭,節點B,E,H處出現了豎直向上的箭頭,表示相應的載荷已經加載到了對應的節點上。

3.2有限元計算

3.2.1計算結果

將此懸臂式桁架結構的有限元模型提交,ABAQUS軟件根據有限元法計算出了此桁架結構中各個節點在Y方向的位移、鉸支座K和L的支反力以及各個桿件所受的應力,如圖5所示。

由圖5 a)—c)左上角的圖例可知,桿件在某處的顏色越深,代表在該處桿件的相關物理量的數值越大。從圖5 a)中看出,桁架的整體結構發生了向上的位移,且懸臂端的幾個節點發生的位移較大;從圖5 b)中發現,鉸支座K處的支反力大于鉸支座L處的支反力;由圖5 c)可知,靠近鉸支座K和L的桿件應力要大于遠離鉸支座K和L的桿件應力,且桿件HJ所受的應力最大。

3.2.2結果對比

利用ABAQUS軟件的查詢功能,可以獲得節點A在豎直方向的位移、鉸支座K的支反力以及桿件HJ所受的應力。經查詢,節點A在豎直方向的位移大小為3.478 mm,與表3所示的計算結果之間的相對誤差為0.14%;與此同時,鉸支座K的支反力大小為262 kN,在規定的小于300 kN的范圍內;另外,桿件HJ所受的應力的大小為133.9 MPa,也小于其許用值157 MPa。由此可見,本文所提出的非線性規劃模型的優化結果是合理的。

4結語

通過建立一個以懸臂端節點撓度最小為目標函數,以桿件的長度、橫截面積等結構尺寸、鉸支座所受的約束力和桿件的最大應力為限制條件的非線性規劃模型,初步確定了某懸臂式桁架結構的優化方案,并利用ABAQUS軟件,根據有限元法判斷該優化方案是否合理。主要的結論如下。

1)當水平方向相鄰2節點之間的距離為1 m,兩鉸支座之間的距離為2.6 m,且桿件的橫截面積均為1 400 mm2時,此懸臂式桁架的懸臂端節點的撓度最小。

2)ABAQUS軟件利用有限元法計算出的結果與Python的計算結果基本吻合,且符合限制條件,因此可以證明本文提出的優化模型與優化方案都是合理的。

3)本文提出的優化模型不僅對優化懸臂式桁架結構具有一定的參考價值,也為解決其他類型的桁架結構(如簡支式、外伸式或空間桁架)的優化問題提供了新思路。

4)在實際工程中,懸臂式桁架結構可能還會受到其他條件的制約,因而本文提出的優化模型中的限制條件還需要在未來的研究中進一步補充,以獲得更為合適的優化結果。

參考文獻/References:

[1]耿雪宵,黃旭就. 基于ABAQUS軟件的三維桁架有限元分析[J]. 煤礦機械,2012,33(4):113115.

GENG Xuexiao, HUANG Xujiu. Finite element analysis of threedimensional truss based on ABAQUS software[J]. Coal Mine Machinery, 2012, 33(4): 113115.

[2]張國偉,喬東需,丁夢婷,等. 基于ABAQUS的M型型鋼桁架樓板有限元分析[J]. 鋼結構,2019,34(11):6064.

ZHANG Guowei, QIAO Dongxu, DING Mengting, et al. Finite element analysis of Msection steel truss floor based on ABAQUS[J]. Steel Construction, 2019, 34(11): 6064.

[3]袁國平,蔣永揚,吳曉宇. 大跨度拱形拉索桁架結構封閉煤場關鍵施工技術[J]. 鋼結構,2019,34(12):108113.

YUAN Guoping, JIANG Yongyang, WU Xiaoyu. Key construction techniques of a closed coal storage with longspan arched cable truss structure[J]. Steel Construction, 2019, 34(12): 108113.

[4]廖能解,馬平,李健洪,等. 上下料桁架機器人仿真研究[J]. 機床與液壓,2019,47(21):914.

LIAO Nengjie, MA Ping, LI Jianhong, et al. Simulation study on loading and unloading material truss robot[J]. Machine Tool & Hydraulics, 2019, 47(21): 914.

[5]李元章,魯國富,任三元,等. 基于積木式的飛艇桁架式復合材料龍骨結構驗證方法[J]. 復合材料科學與工程,2020(1):6771.

LI Yuanzhang, LU Guofu, REN Sanyuan, et al. Verification method for composite keel structure of airship truss based on building blocks[J]. Composites Science and Engineering, 2020 (1): 6771.

[6]崔宏林,李輝. 小型雙體船桁架式連接橋結構強度評估及優化[J]. 艦船科學技術,2019,41(12):1619.

CUI Honglin, LI Hui. Strength assessment and optimization of trussed corss structure for a small catamaran[J]. Ship Science and Technology, 2019, 41(12): 1619.

[7]喬心州,王兵,彭先龍. 桁架結構非概率可靠性形狀優化設計[J]. 應用力學學報,2020,37(1):176182.

QIAO Xinzhou, WANG Bing, PENG Xianlong. Nonprobabilistic reliabilitybased shape optimization design of truss structures[J]. Chinese Journal of Applied Mechanics, 2020, 37(1): 176182.

[8]郝寶新,周志成,曲廣吉,等. 桁架結構拓撲優化的半定規劃建模與求解[J]. 哈爾濱工業大學學報,2019,51(10):1121.

HAO Baoxin, ZHOU Zhicheng, QU Guangji, et al. Modeling and solving of truss topology optimization problems based on semidefinite programming[J]. Journal of Harbin Institute of Technology, 2019, 51(10): 1121.

[9]李彥蒼,王旭. 基于信息熵的改進海豚群算法及其桁架優化[J]. 重慶大學學報,2019,42(5):7685.

LI Yancang, WANG Xu. Improved dolphin swarm algorithm based on information entropy and its truss optimization[J]. Journal of Chongqing University, 2019, 42(5): 7685.

[10]李沛豪,李東. 基于改進粒子群算法的桁架結構優化設計[J]. 空間結構,2018,24(4):1622.

LI Peihao, LI Dong. Optimal design of truss structure based on the improved particle swarm optimization algorithm[J]. Spatial Structures, 2018, 24(4): 1622.

[11]王俊嶺,鐘敬康,楊明霞,等. 基于非線性規劃的綜合管廊斷面尺寸設計優化[J]. 隧道建設,2019,39(12):20302035.

WANG Junling, ZHONG Jingkang, YANG Mingxia, et al. Design optimization of standard crosssection size of utility tunnel based on nonlinear programming[J]. Tunnel Construction, 2019, 39(12): 20302035.

[12]劉鴻文. 材料力學Ⅱ[M]. 6版. 北京:高等教育出版社,2017.

[13]哈爾濱工業大學理論力學教研室. 理論力學Ⅰ[M]. 8版. 北京:高等教育出版社,2016.

[14]劉慶剛,魏青,韓偉信,等. 基于有限元法的V 型缺口平板應力集中系數研究[J]. 河北工業科技,2019,36(4):240245.

LIU Qinggang, WEI Qing, HAN Weixin, et al. Study of the stress concentration factors of a Vnotched plate by using finite element method[J]. Hebei Journal of Industrial Science and Technology, 2019, 36(4): 240245.

[15]QI Zhenchao, LIU Yong, CHEN Wenliang. An approach to predict the mechanical properties of CFRP based on crossscale simulation[J]. Composite Structures, 2019, 210(15): 339347.

猜你喜歡
有限元法
正交各向異性材料裂紋疲勞擴展的擴展有限元法研究
三維有限元法在口腔正畸生物力學研究中發揮的作用
集成對稱模糊數及有限元法的切削力預測
有限元法在機械設計方向中的教學實踐
基于HCSR和CSR-OT的油船疲勞有限元法對比分析
船海工程(2013年6期)2013-03-11 18:57:25
有限元法在基坑邊坡穩定分析中應用
主站蜘蛛池模板: 强乱中文字幕在线播放不卡| 亚洲精品波多野结衣| 人人爽人人爽人人片| 国产精品.com| 久青草国产高清在线视频| 国产91麻豆视频| 欧美亚洲另类在线观看| 全色黄大色大片免费久久老太| 极品国产在线| 97人妻精品专区久久久久| 国产制服丝袜91在线| 国模极品一区二区三区| 国产Av无码精品色午夜| 国产在线视频欧美亚综合| 美女一级毛片无遮挡内谢| 欧美国产日韩在线观看| 国模沟沟一区二区三区 | 99伊人精品| 中文字幕av一区二区三区欲色| 88av在线看| 亚洲av无码成人专区| 美女内射视频WWW网站午夜 | 国产福利一区视频| 丁香婷婷激情网| 国产精品无码影视久久久久久久 | 国产主播喷水| 青青草91视频| 欧美午夜久久| 亚洲综合18p| 欧美中出一区二区| 欧美区国产区| 亚洲成aⅴ人在线观看| 在线看片国产| 国内毛片视频| 国产清纯在线一区二区WWW| 免费jizz在线播放| 国产69囗曝护士吞精在线视频| 女人爽到高潮免费视频大全| 98超碰在线观看| 精品国产www| 中文字幕伦视频| 欧美亚洲国产视频| 亚洲视频免费在线看| 国产主播在线一区| 亚洲国产天堂久久九九九| 黄色免费在线网址| 538国产视频| 午夜精品国产自在| 婷婷亚洲综合五月天在线| 国产又粗又猛又爽视频| 国产精品一老牛影视频| 久久天天躁狠狠躁夜夜躁| 久久久久久久久久国产精品| 亚洲 欧美 日韩综合一区| 日韩免费成人| 国产成人1024精品下载| 亚洲av成人无码网站在线观看| 国产精品亚洲综合久久小说| 国产亚洲日韩av在线| 欧美日韩一区二区在线免费观看| 亚洲综合久久一本伊一区| 四虎精品黑人视频| 久久人人妻人人爽人人卡片av| 久久精品无码专区免费| 国产亚洲精品无码专| 精品视频在线一区| 亚洲欧美国产高清va在线播放| 成人蜜桃网| 九月婷婷亚洲综合在线| 亚洲妓女综合网995久久 | 天堂av高清一区二区三区| 亚洲无码精品在线播放| 一级爱做片免费观看久久| 亚洲国产一区在线观看| 欧美亚洲国产视频| 在线观看无码a∨| 国产精品成人免费视频99| 又粗又大又爽又紧免费视频| 欧洲精品视频在线观看| 久久99精品久久久大学生| 国产好痛疼轻点好爽的视频| 2020久久国产综合精品swag|