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

船用起重機多繩減搖系統的非線性穩定控制

2025-08-29 00:00:00趙庭祺程宏宇孫茂凱王生海王浩然陳海泉
中國機械工程 2025年7期

中圖分類號:U664

DOI:10.3969/j.issn.1004-132X.2025.07.011

開放科學(資源服務)標識碼(OSID):

Nonlinear Stability Control of Multi-cable Anti-sway Systems for Marine Cranes

ZHAO Tingqi1 CHENG Hongyu2SUN MaokailWANG Shenghai1* WANG Haoran1CHEN Haiquan1 1.College of Marine Engineering,Dalian Maritime University,Dalian,Liaoning,116026 2.AECC Shenyang Liming Aero-Engine Co.,Ltd.,Shenyang,110043

Abstract: In order to improve the robustness of mechanical anti-sway control algorithms for marine cranes,and avoid the overload of anti-sway cables caused by the ineffective synchronization motions between anti-sway cables and main hoisting cable in multi-cable anti-sway systems,a fast terminal sliding mode nonlinear anti-sway controler with improved reaching law(FTSMASC-IRL) and a fuzzy adaptive sliding mode synchronization controllr(FASMSC) were developed. The results of simulation experiments show that, FTSMASC-IRL reduces the payload swing angles by 80% ,while the constant tension controller reduces the payload swing angles by 70% . Compared to the traditional augmented PD controller(APDC),FSMASC improves the performance in error control over 75% :

Key words: nonlinear anti-sway controller; synchronization controller; multi-cable anti-sway system;marine crane

0 引言

船用起重機是船舶高效裝卸和轉運貨物的重要工具。由于工作環境的特殊性,相對于陸地起重機,船用起重機承受的外部激勵更加復雜,因此吊重擺動的抑制更加困難。目前,減搖方式分為電子式減搖和機械式減搖[。電子式減搖通常采用多種傳感器、執行器、外加電源來達到減搖的效果,而不改變船用起重機的機械結構[2-4]。電子式減搖的功率配置要求較高,不適用于現有船用起重機。機械式減搖通過改變船用起重機的機械結構來達到減搖的目的。SUN等[5引入多纜減搖動系統(MCAS)來改善船用起重機的動態性能,但使用的是傳統的恒張力控制器。機械式減搖的現有研究集中于結構的安裝與布置[6-9],對控制算法的研究較少。

船用起重機采用多繩減搖系統時,如果不能有效控制主吊索與減搖索的同步運動,容易使減搖索拉偏、主吊索松動,造成減搖索過載甚至斷裂,嚴重影響船用起重機的正常操作。

繩驅并聯機器人(cable-drivenparallelrobot,CDPR)具有多個運動自由度并能在大工作空間內實現高速運動[10-12]。ETIENNE等[13]研究的繩索驅動并聯起重機不僅降低了起重機的重量和慣性,而且沒有寄生傾斜。LYU等[14]設計出一種具有抗擺動特性的六自由度纜驅并聯機器人。TEMPEL等[15]建立了帶有細長桿結構的多繩索系統動力學模型。AMARE等[16]在多繩索系統的動力學模型中引入極小值優化算法,有效求解了索張力和索長度的優化問題。HAMANN等[17]提出一種用于三自由度可重構繩索并聯機器人的幾何參數辨識模型驅動方法。上述研究集中于陸地式多繩索機械設備,對船舶多繩索機械設備上的CDPR研究較少。

機械式減搖的控制通常為線性控制。采用多繩減搖系統的船用起重機要保證主吊索與減搖索在減搖過程中同步運動,因此開發非線性減搖控制器來實現主吊索與減搖索的同步運動十分必要。本文首先建立海洋環境下船用起重機的剛柔耦合模型,然后設計出帶有改進趨近律的快速終端滑模非線性減搖控制器(fastterminalslidingmode nonlinear anti-sway controller with im-provedreachinglaw,FTSMASC-IRL)和模糊自適應滑模同步控制器(fuzzyadaptiveslidingmodesynchronizationcontroller,FASMSC)。仿真試驗證明,FTSMASC-IRL可顯著減小吊重的擺動并提高機械式減搖控制算法的魯棒性,FASMSC可實現主吊索與減搖索的同步運動,避免主吊索松動導致的減搖索過載甚至斷裂。

1配備多繩減搖系統的船用起重機模型

1.1 運動學模型

圖1中, O0X0Y0Z0 為慣性坐標系,O1X1Y1Z1 為受到外部激勵后的船舶坐標系,O2X2Y2Z2 為船用起重機坐標系, O2'X2'Y2'Z2' 為執行回轉作業后的船用起重機坐標系。設定船用起重機先橫搖再縱搖;起重機的吊重先在O2Y2Z2 面內擺動(產生的擺角為面內角 θ1 ),再在 O2X2Z2 面內擺動(產生的擺角為面外角 θ2 )。H 點為變幅轉動軸心, K 為減搖臂與主吊臂的交點, M,N,S 為減搖臂與減搖索的交點,減搖臂與主吊臂相互垂直。 D 點為吊臂頭點, P 點為吊重點,本文忽略繩索的彈性和質量,并將吊重設為質量為 Ψm 的質點。3根減搖索和主吊索共同組成船用起重機的多繩減搖系統。

圖1配備多繩減搖系統的船用起重機模型 Fig.1 Marine crane model with multi-cable anti-swaysystem

在船用起重機坐標系中, P 點的位置矢量為

ZD-lcosθ2cosθ1]T

式中:為主吊索的長度。

同理,得到點 D,M,N,S 的位置矢量:

[LKMsinβLHKcosα+LKMcosβLOH+LHKsinα]T

其中, L 表示兩點的間距, L 的下標為間距對應的兩點。減搖索的長度矢量為

式中: q0 為減搖索的初始長度; q1,0…q2,0…q3,0 分別為減搖索 1~3 的初始長度; q1,q2,q3 分別為減搖索 1~3 在絞車轉動后的長度; ψ 為絞車的轉動矢量; ψ1、ψ2、ψ3 分別為與減搖索 1~3 的轉動量; N 為絞車的傳動比矩陣; n1 、n2Ωνn3 分別為與減搖索 1~3 的傳動比。

減搖索跟隨主吊索同步運動時,減搖索的期望長度為

式中: (Xi,Yi,Zi?(XP,YP,ZP) 分別為點 i 和點 P 的坐標, Ψi=M,N,S 。對式(7)進行一次微分,得到減搖索的期望速度:

式中: J 為雅可比矩陣; 為點 P 的速度。

1.2 動力學建模

多繩系統動力學分析如圖2所示,其中, F1 !F2、F3 分別為減搖索 1~3 的張力, FR 為主吊索的張力, G 為吊重的重力。為防止3根減搖索過載,設定主吊索承受吊重的 90% 力,減搖索系統承擔吊重的 10% ,因此引入比例系數 Kp (本文設為0.1)。

圖2多繩系統動力學分析

Fig.2Dynamicanalysiofmulti-cablesystem

吊重的動能為

吊重的勢能為

Ep=mgZP

結合式(9)、式(10)及拉格朗日函數可得

拉格朗日方程為

式中: Qi 為對應于位置矢量 P 的廣義力。

將式(11)代人式(12),得到未包含絞車的減搖索動力學模型:

式中: 為吊重的加速度。

帶有絞車的減搖索動力學為

式中: 分別為絞車的慣性矩陣、黏滯摩擦系數矩陣和驅動力矩。

結合式(6)、式(15),重新定義帶有絞車的減搖索動力學模型:

式中: 分別為減搖索的實際速度和實際加速度。

結合式(13)、式(16)易得帶有絞車的減搖索動力學模型:

主吊索受力 FR 在 X2、Y2、Z2 方向的分量為

由牛頓第二定律可得

式中: 分別為吊重在 X2、Y2、Z2 方向的加速度。

對式(1)進行兩次微分可得

將式(14)、式(18)、式 (22)~ 式(24)代人式(19) ~ 式(21),可得角加速度:

2 控制策略

船用起重機只受外部激勵或執行回轉作業時,主吊索長度一般不變化,減搖索不需要跟隨主吊索動作,此時不需要考慮減搖索與主吊索的同步運動,只需要對吊重的擺動進行抑制,因此FASMSC不工作,只需FTSMASC-IRL工作。

船用起重機執行起升和變幅作業時,主吊索長度變化,減搖索需跟隨主吊索同步運動,否則易造成主吊索松繩,導致減搖索拉偏甚至被拉斷。這個過程中,通常先用FASMSC控制減搖索與主吊索的同步運動,再用FTSAMASC-IRL抑制吊重的擺動。

2.1 減搖控制器設計

滑模控制具有較強的魯棒性和克服系統不確定性的能力,可有效抑制吊重的擺動??焖俳K端滑模能保證系統在有限時間內達到穩定狀態。由于傳統的指數趨近律易在原點產生抖振,因此需對趨近律進行改進,提出了采用改進趨近律的快速終端滑模非線性減搖控制器,如圖3所示,圖中, A~D 為包含 θ1 和 θ2 的矩陣, τ1…τ2…τ3 為減搖控制器輸出的力矩 τ 的分量。

圖3減搖控制器框圖

Fig.3Diagram of the anti-sway controller

定義終端滑模面函數:

式中: e 為擺角誤差; 為擺角誤差對時間的導數; c,d,φ 為增益系數, cgt;0,dgt;0,1/3lt;φlt;1;θ,θd 分別為吊重實際的擺角和期望的擺角。

滑??刂破鞯脑O計包括滑模面的選取和趨近律的選擇。傳統的滑模趨近律通常指數趨近律:

式中: 為滑模面 s 對時間的導數: 為增益系數, εgt;0 0kgt;0 。

由于指數趨近律易使系統在原點產生抖振,因此需要對其進行優化,優化后的趨近律為

對式(27)求一次微分可得

系統狀態接近平衡點時,收斂時間主要由 決定;系統狀態遠離平衡點時,收斂時間主要由 決定。

將式(28)代入式(31),得到

結合式(25)式(26)式(31)、式(32),并將 α 設為常數,得到減搖控制器表達式

ψD=[XDψDZDZDψ]T

為證明控制器的穩定性,構建李雅普諾夫函數

對式(38)進行一次微分得

根據李雅普諾夫穩定性定理,減搖控制器是穩

定的。

2.2 同步控制器的設計

多繩系統的減搖索與主吊索的同步控制依賴動力學模型的精度,但受機械加工誤差等因素影響,很難得到一個精確的數學模型。因此,為適應不確定的系統模型并提高對系統參數變化和外部干擾的魯棒性,設計出一種模糊自適應滑模同步控制器,如圖4所示。

為實時調整增益系數 cs ,建立規則庫,通過模糊推理、去模糊化實時調整增益系數 cs 。

設計的模糊規則如表1所示,根據模糊規則,以滑模面 ss 為輸入,以增益系數 cs 為輸出, Z PS、PM、PB分別表示零、正小、正中、正大。

為證明同步控制器的穩定性,構建李雅普諾夫函數

表1 模糊規則表

Tab.1 Rulebaseoftheparameter

對式(44)進行一次微分

根據李雅普諾夫穩定性定理可知同步控制器穩定。

3仿真分析

采用MATLAB/Simulink驗證FASMSC和的 FTSMASC-IRL有效性,設定 LOH、LHD、LHS 、LKM,LKN 分別為 15m,26m,31m,5.5m,5.5 m ,吊重的質量為 10t ,變幅角 α 為 30° 。在FTSMASC-IRL中,設置 c=1,ε=0.5,k=3 。在FASMSC中,設置 分別為0.05、0.6、0.1。

船用起重機為 45t 起重機,如圖5所示。起重機的慣量矩陣 Imt=diag(1.94,1.94,1.94) (20kg?m2 ,摩擦參數矩陣 Fvt=diag(0.2,0.2,0.2) (2號N?s ,傳動比矩陣 N=diag(0.17,0.17,0.17)m 吊重的初始擺角設為 0°

圖545噸船用起重機在船上的工程應用Fig.5The engineering application of 45-ton crane inthe ship

設定海況為4級海況時,海浪的橫搖、縱搖幅度分別為 5° 和 2° ,橫搖、縱搖周期分別為14.8s和6s[18] 。選擇PM波浪譜為計算模型,本文設計的不規則波公式為

其中,有義波高 h 為 4m ,頻率 ω 的變化范圍為0.4~4Hz ,初始相位角 εw=0,g=9.8m/s2 。不規則波形如圖6所示。

圖6 不規則波Fig.6The irregular waves

3.1 減搖控制效果分析

基于空氣阻尼原理的恒張力控制器(constanttension controller,CTC)是機械式減搖中常見的一種控制器。吊重向減搖臂端點 K 擺動時, F3=Fδ3Y ,減搖索1和2為松弛狀態,減搖索3為張緊狀態,抑制吊重擺動;吊重向減搖臂端點 S 擺動時, F22Y,F1=Fδ1Y ,減搖索3為松弛狀態,減搖索1和2為張緊狀態,抑制吊重擺動;吊重向減搖臂端點 N 擺動時, Fδ3X ,減搖索2為松弛狀態,減搖索1和3為張緊狀態,抑制吊重擺動;吊重向減搖臂端點 M 擺動時, 州 F3=FδδδX ,減搖索1為松弛狀態,減搖索2和3為張緊狀態,抑制吊重擺動。其中,Fδ1X、Fδ2X、Fδ3X 和 Fδ1Y?Fδ2Y?Fδ3Y 分別為減搖索1、減搖索2和減搖索3為滿足減搖索保持拉力狀態下的預緊力在 X2 和 Y2 方向上的分量。為驗證FTSMASC-IRL的減搖控制效果,對比控制器CTC和FTSMASC-IRL的擺角抑制效果。

船用起重機在設定不規則波下的回轉角度如圖7所示。

圖7 船用起重機的回轉角度

Fig.7The slewing angle of marine crane

實際工程應用場景中,主吊索太短會影響貨物的轉運。因此,本文探究船用起重機執行回轉作業時的吊重擺動,如圖8所示,此時的主吊索長度(初始長度)分別為 15m.20m.25m 。

主吊索初始長度為 15m 時,吊重的面內角、面外角、吊重在 O2X2Y2 面的運動軌跡如圖9所示。相較于SMC,FTSMASC-IRL的軌跡更加光滑并且收斂更快,證實FTSMASC-IRL可有效減小系統抖振。

圖8主吊索的長度

Fig.8The length range of the main cable

?

圖9吊重的擺動以及吊重的運動軌跡(主吊索長 15m )

Fig.9The swing and the motion trajectory of the payload(length of main hosting cable is as 15m )

主吊索初始長度為 20m.25m 時,吊重面內角、面外角、吊重在 O2X2Y2 面的運動軌跡如圖10、圖11所示。

如圖 9~ 圖11所示,無控制時,主吊索初始長度 15m,20m,25m 的面內角最大值分別為54°,34°,24° ,面外角最大值分別為 28°,30°,26° 。CTC作用下,主吊索初始長度 15m.20m.25m 的面內角最大值分別為 5°,5°,8° ,面外角最大值分別為 12°,12°,8° 。FTSMASC-IRL作用下,主吊索初始長度 15m20m25m 的面內角最大值分別為 3°,2°,1.5° ,面外角最大值約分別 5°,4°,4° 。因此,相對于無控制,FTSMASC-IRL的吊重擺動幅度減小 80% ,CTC的吊重擺動幅度減小70% 。外部激勵相同的情況下,主吊索的長度增加時,吊重的擺幅逐漸減小。

圖10吊重的擺動以及吊重的運動軌跡(主吊索長 20m )

Fig.10The swing and the motion trajectory of the payload(length of main hosting cable isas 20m )

3.2 同步運動的控制效果分析

減搖索的實際位置越接近期望位置就越能實現主吊索與減搖索的同步運動。船用起重機受外部激勵時,通過控制減搖索的長度和速度使減搖索的實際位置接近期望位置。

3.2.1 起升作業的同步運動控制效果分析

本文設定船用起重機執行起升作業時受不規則波的影響,通過控制主吊索的釋放與收縮來完成貨物的起吊作業。起升作業的主吊索長度變化如圖12所示。起升作業下,不控制的減搖索長度誤差、速度誤差如圖13所示,抑制擺角的減搖索長度誤差、繩速誤差如圖14所示。

由圖13、圖14可以看出,在FASMSC的作用下,減搖索的實際位置迅速接近其期望位置,實現了減搖索與主吊索的同步運動。圖14中的減搖索長度誤差和速度誤差產生抖動的原因是FASMSC和FTSAMASC-IRL同時作用于絞車會使絞車頻繁正反轉。

3.2.2 變幅作業的同步運動控制效果分析

本文設定船用起重機執行變幅作業時會受不規則波的影響。變幅作業主要通過改變主吊臂與水平面的夾角來完成變幅動作,設定變幅角的變化速度為 0.5°/s 。變幅作業下,無控制的減搖索長度誤差、速度誤差如圖15所示,抑制擺角的減搖索長度誤差、速度誤差如圖16所示。

由圖15、圖16可以看出,在FASMSC的作用下,減搖索的實際位置迅速接近其期望位置,實現了減搖索與主吊索的同步運動。圖16中,減搖索的長度誤差和速度誤差不能趨于穩定的原因是FASMSC和FTSAMASC-IRL同時作用于絞車會使絞車頻繁的正反轉。

3.2.3 同步運動控制效果對比分析

為驗證FASMSC對主吊索和減搖索的同步運

"

動性能改善的效果,對比了增廣PD控制器(aug-mentedPDcontroller,APDC)與FASMSC的同步運動控制效果。本文設計的增廣PD控制器為

式中: 分別為減搖索長度誤差和減搖索速度誤差的增益矩陣。

圖17為APDC的控制框圖。

由圖18、圖19可知,相對于APDC,FASMSC具有更好的誤差控制性能、更高的收斂速度。起升和變幅作業下的誤差控制性能及收斂速度如表

2~ 表5所示。相對于APDC,FSMASC對減搖索長度誤差的控制性能提高 75% ,在收斂速度至少提高了 60% 。

"
"

4實驗結果

在船用起重機縮比樣機實驗平臺上進行減搖控制實驗。實驗平臺包括船用起重機縮比樣機、液壓單元、六自由度穩定平臺、電子控制柜等,如圖20所示。

"

設定2組實驗的船舶激勵為規則波,其中,起重機受到的橫搖激勵 θ1x 和縱搖激勵 θ1y 均為4sin(πt/3) 和 6sin(πt/3) ,設定吊重的初始擺角為 0° ,吊重的質量為 50kg 。

CTC對吊重擺動的抑制效果如圖21、圖22所示。相對于沒有控制,CTC作用下的吊重擺動幅度減小 70% 。

"

5結論

1)FTSMASC-IRL、CTC的吊重擺幅分別比無控制的吊重擺幅減小 80% 和 70% 。FTSMASC-IRL、CTC都能有效抑制掛載的擺動,但FTSMASC-IRL的魯棒性更高。

2)船用起重機進行回轉作業時,主吊索越長,吊載擺動幅度越小。船用起重機起升和變幅作業時,FASMSC的減搖索與主吊索同步誤差趨于零,有效實現了減搖索與主吊索的同步運動。

3)將FASMSC與APDC的控制效果進行對比可得,FASMSC的最大誤差波動減小 75% ,誤差收斂速度提高 60% ,驗證了FASMSC對主吊索與減搖索的同步運動具有良好的控制效果。

參考文獻:

[1]孫茂凱,王生海,韓廣冬,等.船用起重機多柔索減搖系統的動力學分析與工程應用[J.中國機械工 程,2024,35(7):1308-1317. SUN Maokai,WANG Shenghai, HAN Guangdong,et al. Dynamics Analysis and Engineering Applications of Multi-tagline Anti-swing System for Marine Cranes[J]. China Mechanical Engineering, 2024,35(7):1308-1317.

[2] RAJA ISMAIL R M T,THAT N D,HA Q P. Modelling and Robust Trajectory Following for Offshore Container Crane Systems[J]. Automation in Construction,2015,59:179-187.

[3]LIU Zhuoqing,FU Yu, SUN Ning,et al. Collaborative Antiswing Hoisting Control for Dual Rotary Cranes with Motion Constraints[J]. IEEE Transactions on Industrial Informatics,2022,18(9):6120-6130.

[4]ZHAO Bingqing,OUYANG Huimin, IWASAKI M. Motion Trajectory Tracking and Sway Reduction for Double-pendulum Overhead Cranes Using Improved Adaptive Control without Velocity Feedback [J]. IEEE/ASME Transactions on Mechatronics, 2022,27(5):3648-3659.

[5]SUN Maokai, WANG Shenghai, HAN Guangdong,et al. Multi-cable Anti-swing System for Cranes Subject to Ship Excitation and Wind Disturbance:Dynamic Analysis and Application in Engineering[J]. Ocean Engineering,2023,281:114518.

[6]MARTIN I A, IRANI R A. Evaluation of both Linear and Non-linear Control Strategies for a Shipboard Marine Gantry Crane[C]// OCEANS 2019 MTS/IEEE SEATTLE. Seattle,2019:1-10.

[7]WANG Jianli, WANG Shenghai, CHEN Haiquan, et al.Dynamic Modeling and Analysis of the Telescopic Sleeve Antiswing Device for Shipboard Cranes[J]. Mathematical Problems in Engineering, 2021,2021(1):6685816.

[8]YUAN Guo hui,HUNT B R,GREBOGI C,et al. Design and Control of Shipboard Cranes[C]// ASME 1997 Design Engineering Technical Conferences. Sacramento,1997:10.1115/detc97/vib-4095.

[9]WANG Jianli, WANG Shenghai, CHEN Haiquan, et al.Dynamic Modeling and Analysis of the Telescopic Sleeve Antiswing Device for Shipboard Cranes[J].Mathematical Problems in Engineering, 2021,2021(1):6685816.

[10]IDA E,MATTIONI V. Cable-driven Parallel Robot Actuators:State of the Art and Novel Servowinch Concept[J].Actuators,2022,11(10):290.

[11] CARRICATO M,MERLET J P. Stability Analysis of Underconstrained Cable-driven Parallel Robots[J]. IEEE Transactions on Robotics,2013,29 (1):288-296.

[12]GOUTTEFARDE M,COLLARD JF,RIEHL N,et al.Geometry Selection of a Redundantly Actuated Cable-suspended Parallel Robot[J]. IEEE Transactions on Robotics,2015,31(2):501-510.

[13]ETIENNE L,CARDOU P,METILLON M,et al.Design of a Planar Cable-driven Parallel Crane without Parasitic Tilt[C]// ASME 2021 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference,2021.

[14].LYU Wei, TAO Limin, JI Zhengnan. Research on Anti-swing Characteristic of Redundancy Cabledriven Parallel Robot[C]// 2017 IEEE 2nd Advanced Information Technology, Electronic and AutomationControl Conference ( IAEAC). Chongqing,2017:1504-1508.

[15]TEMPEL P,LEE D,TRAUTWEIN F,et al. Modeling of Elastic-flexible Cables with Time-varying Length for Cable-driven Parallel Robots[C]// Fourth International Conference on Cable-Driven Parallel Robots.Krakow,2019:295-306.

[16]AMARE Z,ZI Bin,QIAN Sen,et al. Three-dimensional Static and Dynamic Stiffness Analyses of the Cable Driven Parallel Robot with Non-negligible Cable Mass and Elasticity[J]. Mechanics Based Design of Structures and Machines,2018, 46(4):455-482.

[17]HAMANN M,NUSSE P M,WINTER D,et al. Towards a Precise Cable-driven Parallel Robot - a Model-driven Parameter Identification Enhanced by Data-driven Position Correction[C]// Fourth International Conference on Cable-Driven Parallel Robots.Krakow,2019:367-376.

[18]LOVE L. Compensation of Wave-induced Motion and Force Phenomena for Ship-based High Performance Robotic and Human Amplifying Systems [R]. Oak Ridge:Oak Ridge National Lab,2003.

(編輯張洋)

作者簡介:趙庭祺,男,1997年生,碩士研究生。研究方向為船用起重機吊裝減搖,發表論文6篇。E-mail:15847880562@163.com。王生海*(通信作者),男,1988年生,副教授。研究方向為起重機吊裝補償技術、繩索驅動機器人。發表論文50余篇。E-mail:shenghai_wang@dlmu.edu.cn。

本文引用格式:

趙庭祺,程宏宇,孫茂凱,等.船用起重機多繩減搖系統的非線性穩定控制[J].中國機械工程,2025,36(7):1487-1496.ZHAO Tingqi,CHENGHongyu,SUNMaokai,etal.NonlinearStability Control ofMulti-cable Anti-sway Systems forMarineCranes[J].ChinaMechanical Engineering,2025,36(7):1487-1496.

主站蜘蛛池模板: 亚洲黄网在线| 一区二区自拍| 视频二区中文无码| 免费播放毛片| 最新国产高清在线| 国产亚洲欧美另类一区二区| 久久中文电影| 黄片一区二区三区| 性欧美精品xxxx| 99久久人妻精品免费二区| 亚洲成AV人手机在线观看网站| 色综合天天综合中文网| 国产成人做受免费视频| 国产乱论视频| 亚州AV秘 一区二区三区| 99精品免费欧美成人小视频| 九月婷婷亚洲综合在线| 91成人试看福利体验区| 999福利激情视频| 谁有在线观看日韩亚洲最新视频 | 国产成人亚洲综合A∨在线播放| 国产精品大尺度尺度视频| 久久午夜夜伦鲁鲁片无码免费 | 日本人妻丰满熟妇区| 爽爽影院十八禁在线观看| 国产乱人免费视频| 美女毛片在线| 欧美精品黑人粗大| 成人噜噜噜视频在线观看| 高潮毛片无遮挡高清视频播放| 婷婷开心中文字幕| 国产午夜看片| 欧美亚洲一二三区| 欧美成一级| 午夜综合网| 国产乱子伦一区二区=| 91探花在线观看国产最新| 热这里只有精品国产热门精品| 婷婷99视频精品全部在线观看| 久久精品国产精品国产一区| 丁香六月综合网| 国产一区二区在线视频观看| 女人18毛片水真多国产| 欧美视频在线观看第一页| 欧美三级视频网站| 亚洲天堂视频网站| 97精品久久久大香线焦| 91 九色视频丝袜| 99热这里只有精品免费| 秋霞午夜国产精品成人片| 欧美成人午夜视频免看| 无码一区二区波多野结衣播放搜索| 免费看的一级毛片| 国内精品自在自线视频香蕉| JIZZ亚洲国产| 亚洲欧美色中文字幕| 亚洲综合一区国产精品| 99国产在线视频| 欧美日韩一区二区三| 亚洲第一色网站| 在线亚洲天堂| 国产啪在线91| 色综合天天娱乐综合网| 国产夜色视频| 亚洲欧美成人在线视频| 狠狠色香婷婷久久亚洲精品| 亚洲精品波多野结衣| 中文国产成人久久精品小说| 国产在线精品人成导航| 日韩大乳视频中文字幕 | 亚洲va视频| 国产青榴视频| 中文无码日韩精品| 久久久久久高潮白浆| 国产污视频在线观看| 最新精品国偷自产在线| 国产精品999在线| 中文字幕一区二区人妻电影| 国产97视频在线| 爆乳熟妇一区二区三区| 3D动漫精品啪啪一区二区下载| 久久99热这里只有精品免费看 |