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

磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型復合辨識方法

2017-06-05 14:20:15張華強王英廣趙學濤
中國慣性技術(shù)學報 2017年2期
關(guān)鍵詞:模態(tài)有限元分析

張華強,王英廣,趙學濤

(1. 山東理工大學 機械工程學院,淄博 255049;2. 北京控制工程研究所,北京 100094)

磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型復合辨識方法

張華強1,王英廣2,趙學濤1

(1. 山東理工大學 機械工程學院,淄博 255049;2. 北京控制工程研究所,北京 100094)

針對磁懸浮撓性轉(zhuǎn)子的純有限元分析模型精度低的問題,提出一種磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型復合辨識方法。該方法首先采用有限元方法把磁懸浮撓性轉(zhuǎn)子劃分為多個Timoshenko梁單元,建立磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型,并對其撓性臨界頻率、阻尼、剛度和振型等關(guān)鍵參數(shù)進行分析,進而得到磁懸浮撓性轉(zhuǎn)子的等效降階數(shù)學模型;然后采用魯棒自適應方法分析磁懸浮撓性轉(zhuǎn)子系統(tǒng)的動態(tài)特性;最后,采用變LEVY方法從動態(tài)特性分析數(shù)據(jù)中對磁懸浮撓性轉(zhuǎn)子進行系統(tǒng)辨識,校正有限元分析得到的降階數(shù)學模型。實驗結(jié)果表明,本方法可以得到磁懸浮撓性轉(zhuǎn)子較為準確的系統(tǒng)模型。

撓性轉(zhuǎn)子;有限元;魯棒自適應;變LEVY方法;復合辨識

在磁懸浮離心機、渦輪機等超高速應用領(lǐng)域,轉(zhuǎn)子一般設(shè)計成撓性的[1]。磁懸浮撓性轉(zhuǎn)子系統(tǒng)是一個復雜的機電一體化系統(tǒng),包括傳感器、控制器、功率放大器、執(zhí)行器(磁軸承)和撓性轉(zhuǎn)子等組成部件[2-3]。獲知磁懸浮轉(zhuǎn)子系統(tǒng)模型是進行后續(xù)過臨界轉(zhuǎn)速控制的前提條件[4-5]。

磁懸浮轉(zhuǎn)子剛性模型只是轉(zhuǎn)子在低轉(zhuǎn)速時的一種近似情況,是建立在轉(zhuǎn)子自身無形變這一前提條件之下的,并只有六個自由度[6-7]。在進入撓性轉(zhuǎn)速后,轉(zhuǎn)子上任意兩點的相對空間位置都是隨時間變化的,其自身的形變不可忽略。撓性轉(zhuǎn)子模型還受加工精度和裝配條件影響較大,即便微小加工裝配的誤差,也會形成較大的模型差異[8]。

本文首先采用有限元方法把磁懸浮撓性轉(zhuǎn)子劃分為多個Timoshenko梁單元,建立磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型,并對其撓性臨界頻率、阻尼、剛度和振型等關(guān)鍵參數(shù)進行分析,得到轉(zhuǎn)子系統(tǒng)的等效降階數(shù)學模型;然后,使用魯棒自適應方法對磁懸浮撓性轉(zhuǎn)子系統(tǒng)進行動態(tài)特性分析;最后,采用變LEVY方法從動態(tài)特性分析數(shù)據(jù)中對磁懸浮撓性轉(zhuǎn)子進行系統(tǒng)辨識,校正有限元分析得到的降階數(shù)學模型,從而獲得較為精確的磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型。

1 磁懸浮撓性轉(zhuǎn)子系統(tǒng)建模

1.1 磁懸浮撓性轉(zhuǎn)子有限元模型的建立及特性分析

采用有限元方法把磁懸浮撓性轉(zhuǎn)子分割成一個個小單元,通過經(jīng)典轉(zhuǎn)子模型建模,再加上有限元之間的連接關(guān)系就組成單個有限元模型。通過一個個有限元模型的組成可以構(gòu)建描述整個轉(zhuǎn)子模型的矩陣,通過分析該矩陣即可獲得磁懸浮撓性轉(zhuǎn)子的特性參數(shù)。磁懸浮電機轉(zhuǎn)子結(jié)構(gòu)圖如圖1所示。

圖1 磁懸浮電機轉(zhuǎn)子結(jié)構(gòu)圖Fig.1 Rotor structure of magnetic levitated motor

單個有限元模型通常使用工程力學中的梁模型,目前的梁理論主要有Euler-Bernoulli梁和Timoshenko梁[9]。由于Timoshenko梁模型還涉及回轉(zhuǎn)慣性和梁截面剪切變形,其表達式包含梁橫截面回轉(zhuǎn)所產(chǎn)生的動能和剪切變形所引起的彈性勢能,因此采用Timoshenko梁模型可以得到很好的分析結(jié)果。

從磁懸浮撓性轉(zhuǎn)子中劃分 Timoshenko梁單元及其受力分析示意圖如圖2所示。

從圖2中得到轉(zhuǎn)軸撓曲q(z,t)動力學模型如下:

式中,E為彎曲彈性模量,I為有限元截面對中心軸的慣性矩,ρ為有限元材料密度,A為截面面積,Gshear為剪切彈性模量。式(1)中,第一項為彎曲變形勢能項,第二項為徑向運動動能項,第三項為主轉(zhuǎn)動慣量動能項,第四項為主剪切變形勢能項,第五項為合并的剪切變形和轉(zhuǎn)動慣量耦合項。

圖2 單個有限元劃分及其受力分析Fig.2 Single finite element division and its stress analysis

圖3 AR4模型描述單個有限元Fig.3 Single finite element described by AR4 model

式(2)后兩項分別為磁軸承項和不平衡擾動項。根據(jù)Timoshenko梁理論的AR4模型,可把式(2)變換成下列表示形式:

傳感器輸出為

式中,S為觀測矩陣。設(shè)式(3)特解的形式為

代入式(3)得其特征方程為

因陀螺耦合陣G不能對角化,從式(6)解得的特征值為復數(shù)形式,虛部對應轉(zhuǎn)子模態(tài)頻率。

對應的特征向量Φn1、Φn2(n=1,2,…,N)亦為復數(shù)形式,表現(xiàn)為空間三維曲線。

使用RotFE軟件對4 kW磁懸浮電機轉(zhuǎn)子進行有限元分析。將轉(zhuǎn)子劃分為55個有限元,前兩階撓性模態(tài)分析結(jié)果如圖4、圖5所示,一階彎曲模態(tài)頻率為691 Hz,二階彎曲模態(tài)頻率為1 442 Hz。如式(6)所示,隨轉(zhuǎn)速變化,轉(zhuǎn)子的模態(tài)頻率因陀螺效應出現(xiàn)分叉現(xiàn)象,升速時與前向模態(tài)頻率(高頻)相交,降速時與后向模態(tài)頻率(低頻)相交。使用RotFE軟件繪制Compbell圖如圖6所示,升速臨界轉(zhuǎn)速43 000 r/min,降速臨界轉(zhuǎn)速在41 000 r/min。

圖4 撓性轉(zhuǎn)子一階彈性模態(tài)空間運動軌跡Fig.4 Space motion trajectory of the first-order elastic modal

圖5 撓性轉(zhuǎn)子二階彈性模態(tài)空間運動軌跡Fig.5 Space motion trajectory of the second-order elastic modal

圖6 4 kW電機轉(zhuǎn)子Campbell圖Fig.6 Campbell diagram of 4 kW motor rotor

1.2 磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型的等效降階處理

首先,將特征向量組成矩陣Φ:

高階彈性模態(tài)頻率一般遠高于系統(tǒng)額定轉(zhuǎn)速,而且因眾多低通延時環(huán)節(jié)的存在,系統(tǒng)會對高階模態(tài)失去控制能力。另外,真實系統(tǒng)中往往只有前兩到三階彈性模態(tài)對系統(tǒng)穩(wěn)定性有影響,取Φ前r列組成rΦ對式(8)所描述系統(tǒng)進行降階處理:

得到系統(tǒng)動力學方程:

觀測方程為:

上述轉(zhuǎn)換矩陣除Gr外均為對角形式,對式(10)進行拉氏變換,得到系統(tǒng)頻域方程:

The Development of County Rural Tourism under the Background of the All-For-One Tourism——A Case Study of Fengning Manchu Autonomous County,Hebei Province____________________________XU Xinguo,WANG Wenxuan 7

從式(14)可見,轉(zhuǎn)子特征頻率處運動方程為:

圖7所示為全階系統(tǒng)和降階系統(tǒng)不平衡相應比較,激勵點旋轉(zhuǎn)在最左端,響應點選擇在轉(zhuǎn)子中心,系統(tǒng)模型階數(shù)從55降為8,降階后只包含兩個剛性模態(tài)和兩個撓性模態(tài)。在額定轉(zhuǎn)速的兩倍轉(zhuǎn)速120000r/min 以內(nèi),不平衡響應函數(shù)幾乎完全重合,而且功放帶寬有限,系統(tǒng)對高階模態(tài)幾乎不響應。因此,使用降階模型已經(jīng)足夠精確,但計算量大大降低。

圖7 全階轉(zhuǎn)子和降階轉(zhuǎn)子模型的不平衡響應Fig.7 Unbalance response of full-order rotor and reduced-order rotor model

系統(tǒng)模型寫成狀態(tài)方程的形式為:

2 磁懸浮撓性轉(zhuǎn)子系統(tǒng)辨識方法

上文對磁懸浮撓性轉(zhuǎn)子的模型進行了理論分析。撓性轉(zhuǎn)子模型受加工精度和裝配條件影響較大,最終導致基于線性—彈性模型的有限元分析精度較低,造成分析結(jié)果與實際結(jié)果有較大誤差。因此,需要使用系統(tǒng)辨識方法對有限元模型進行校正。

系統(tǒng)辨識采用頻域方法,并使用頻率分析手段區(qū)分各個模態(tài),對包含多模態(tài)的復雜系統(tǒng)有較好的效果。頻率響應方法(FRF)使用掃頻激勵方式,可選擇性地在各個頻率給轉(zhuǎn)子提供激勵,給予各個模態(tài)以更多細節(jié)上的反應[10]。FRF模態(tài)辨識方法,由兩部分組成,一是測得系統(tǒng)的頻率響應,二是從頻率響應中提出模態(tài)參數(shù)。下面將對其分別進行介紹。

2.1 磁懸浮撓性轉(zhuǎn)子系統(tǒng)動態(tài)特性分析

如圖8所示,磁懸浮轉(zhuǎn)子頻率響應主要包括激勵和辨識兩個環(huán)節(jié)。激勵采用正弦掃頻信號形式,頻率變化可選擇遞乘形式。

在激勵接入點之后的被測轉(zhuǎn)子系統(tǒng)兩端引出兩個受激勵的信號,送入自適應辨識環(huán)節(jié)。自適應辨識環(huán)節(jié)采用基于魯棒自適應的方法對磁懸浮撓性轉(zhuǎn)子系統(tǒng)的動態(tài)特性進行分析,其結(jié)構(gòu)圖如圖9所示。

自適應辨識環(huán)節(jié)通過兩個引出信號與激勵信號間的幅值之比和相位之差,得到被測系統(tǒng)的動態(tài)特性。采用LMS自適應算法提取被辨識信號中的同頻成分,擬合公式如下:

圖8 磁懸浮轉(zhuǎn)子頻率響應分析結(jié)構(gòu)圖Fig.8 Analysis structure of magnetic flexible rotor frequency response

圖9 魯棒自適應辨識結(jié)構(gòu)圖Fig.9 Robust adaptive identification structure

LMS是一種遞歸收斂方法,式(17)的自適應步長采用固定步長形式。μ取值大,收斂速度快,但收斂誤差大;μ取值小,收斂誤差小,但收斂速度慢。而且不同頻率端對收斂誤差的要求不同。撓性模態(tài)附近頻率響應函數(shù)變化快,要求收斂速度盡量快,但頻率響應幅值大,對收斂精度敏感度低,要求μ取大。而在其他頻率段,頻率響應幅值小,辨識結(jié)果對收斂精度很敏感,要求μ取小。為平衡收斂速度和辨識誤差間的關(guān)系,在此使用魯棒變步長 LMS算法,以辨識誤差均方差最小為原則,實時調(diào)整步長。調(diào)整公式為

式中,εmin為擬合誤差的均方差。

在原始信號中疊加有大量噪聲,影響辨識精度,需要對其進行預濾波處理。濾波器在二階Butterworth陷波器進行改進,激勵信號頻率f和通帶寬度控制參數(shù)sharp計算得到二階帶通濾波器系數(shù)(A,B),計算公式為

其中,A=(A[0],A[1],A[2])為二階帶通濾波器分母多項式的系數(shù),B=(B[0],B[1],B[2])為二階帶通濾波器分子多項式的系數(shù)。

得到被辨識模塊輸入輸出信號在激勵正余弦信號構(gòu)造的直角坐標系中的投影坐標后,計算相頻幅頻變化公式如下:

得到控制電流到轉(zhuǎn)子位移的幅頻相頻曲線如圖10所示。為更好地對比魯棒自適應和所設(shè)計單通濾波器的效果,圖10將普通自適應、魯棒自適應、帶通濾波+魯棒自適應進行對比。普通自適應方法噪聲大,所產(chǎn)生的累積誤差導致高頻幅值偏大。魯棒自適應在跟蹤速度上雖然比普通自適應方法略差,但噪聲明顯減小,高頻累積誤差大幅降低。從圖 10可見,加入帶通濾波器后,噪聲進一步降低。

圖10 磁懸浮轉(zhuǎn)子頻率響應結(jié)果Fig.10 Frequency response of magnetic flexible rotor

2.2 基于變LEVY方法的系統(tǒng)辨識方法

得到轉(zhuǎn)子系統(tǒng)頻率特性后,獲知轉(zhuǎn)子模型還需從頻率響應函數(shù)中提取轉(zhuǎn)子的模態(tài)參數(shù),得到轉(zhuǎn)子系統(tǒng)模型,本文使用變LEVY方法實現(xiàn)這一過程[11]。轉(zhuǎn)子系統(tǒng)傳遞函數(shù)如下:

上式簡寫作:

因[Z(s)]非奇異可逆,由式(24)可得:

式中:

式(28)略去了下標i j,式中,

為使誤差方程線性化,將式(30)乘以Dk,稱為加權(quán)誤差:

Levy法矩陣方程變?yōu)?/p>

求解如上方程得到{a}、{b}之后即可解得轉(zhuǎn)子模態(tài)參數(shù)。

3 系統(tǒng)模型復合辨識方法的效果及其分析

使用變LEVY方法對上述測試結(jié)果進行辨識。為提高辨識精度,采用分段擬合方法。首先對一階模態(tài)頻率附近區(qū)域進行辨識,辨識結(jié)果如圖11所示。測試結(jié)果與辨識結(jié)果相差較大,主要是傳感器與磁軸承軸向位置不同所致。對初始相位角修正-213°后,得到如圖12所示的辨識結(jié)果,修正后的辨識結(jié)果與測試結(jié)果吻合度很好。一階模態(tài)辨識結(jié)果如下:

同理,對二階模態(tài)頻率附近幅相曲線進行辨識得到如下結(jié)果,擬合曲線見圖13。

圖11 未修正初始相角時轉(zhuǎn)子一階撓性模態(tài)辨識效果圖Fig.11 The first-order flexible modal identification effect without correcting the initial phase angle

圖12 修正初始相角后轉(zhuǎn)子一階撓性模態(tài)辨識效果圖Fig.12 The first-order flexible modal identification effect after correcting the initial phase angle

圖13 磁懸浮轉(zhuǎn)子二階撓性模態(tài)辨識效果圖Fig.13 The second-order flexible modal identification effect of magnetic flexible rotor

對低頻段附近幅相曲線進行辨識得到如下結(jié)果,擬合曲線見圖14。

將低頻段,一階模態(tài)、二階模態(tài)進行綜合得到如下傳遞函數(shù):

圖14 磁懸浮轉(zhuǎn)子低頻段辨識效果圖Fig.14 Identification effect in low-frequency range of the magnetic flexible rotor

測試幅相特性曲線與辨識幅相特性曲線比較見圖 15。本文目的是穿越磁懸浮撓性轉(zhuǎn)子的臨界轉(zhuǎn)速,關(guān)心撓性模態(tài)頻率附近的模型精度,而對其它頻率段的辨識精度要求不高。從圖15辨識結(jié)果看,撓性模態(tài)頻率附近辨識模型與測試結(jié)果擬合程度很高。

圖15 磁懸浮轉(zhuǎn)子辨識結(jié)果圖Fig.15 Identification effect of the magnetic flexible rotor

4 結(jié) 論

獲知磁懸浮轉(zhuǎn)子系統(tǒng)模型是進行后續(xù)過臨界轉(zhuǎn)速控制的前提條件。本文針對磁懸浮撓性轉(zhuǎn)子的純有限元分析模型精度低的問題,提出了一種磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型復合辨識方法。該方法首先采用有限元方法對磁懸浮撓性轉(zhuǎn)子進行了多個 Timoshenko梁單元的劃分,并建立磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型,對其撓性臨界頻率、阻尼、剛度和振型等關(guān)鍵參數(shù)進行分析,通過深入分析得到了磁懸浮撓性轉(zhuǎn)子的等效降階數(shù)學模型;然后采用魯棒自適應方法分析了磁懸浮撓性轉(zhuǎn)子系統(tǒng)的動態(tài)特性;最后,采用變LEVY方法從動態(tài)特性分析數(shù)據(jù)中對磁懸浮撓性轉(zhuǎn)子進行了系統(tǒng)辨識,校正有限元分析得到的降階數(shù)學模型。實驗結(jié)果表明,該方法可以得到較為準確的磁懸浮撓性轉(zhuǎn)子系統(tǒng)模型,具有很強的理論研究價值和工程意義。

(References):

[1] Sun Zhen, He Ying, Zhao Jing-jing, et al. Identification of active magnetic bearing system with a flexible rotor[J]. Mechanical Systems and Signal Processing, 2014, 49(1-2): 302-316.

[2] 房建成, 張會娟, 劉虎. 磁懸浮剛性轉(zhuǎn)子系統(tǒng)振動機理分析與動力學建模[J]. 控制理論與應用, 2014, 31(12): 1707-1713. Fang Jian-cheng, Zhang Hui-juany, Liu Hu. Vibration mechanism analysis and dynamic model development of magnetically suspended rigid rotor system[J]. Control Theory & Applications, 2014, 31(12): 1707-1713.

[3] Zhao Jie, Zhang Hai-tao, Fan Ming-can, et al. Control of a constrained flexible rotor on active magnetic bearings [J]. IFAC-Papers On Line, 2015, 48(28): 156-161.

[4] 黃梓嫄, 韓邦成, 周銀鋒. 非線性接觸下磁懸浮電機柔性轉(zhuǎn)子系統(tǒng)模態(tài)分析[J]. 中國電機工程學報, 2014, 34(15): 2438-2444. Huang Zi-yuan, Han Bang-cheng, Zhou Yin-feng. Modal analysis of the flexible rotor system of magnetic levitation motors under nonlinear contact[J]. Proceedings of the CSEE, 2014, 34(15): 2438-2444.

[5] Wang Y G, Fang J C, Zheng S Q. A field balancing technique based on virtual trial-weights method for a magnetically levitated flexible rotor[J]. Journal of Engineering for Gas Turbines and Power, 2014, 136(9): 1-7.

[6] 劉超, 劉剛, 蓋玉歡. 基于磁力等效原理的剛性磁懸浮轉(zhuǎn)子系統(tǒng)高精度在線動平衡[J]. 振動與沖擊, 2016, 35(4): 67-71. Liu Chao, Liu Gang, Ge Yu-huan. Field balancing for a magnetically suspended rigid rotor based on magnetic forces equivalence principle[J]. Journal of vibration and shock, 2016, 35(4): 67-71.

[7] Tiwari R, Chougale A. Identification of bearing dynamic parameters and unbalance states in a flexible rotor system fully levitated on active magnetic bearings[J]. Mechatronics, 2014, 24(3): 274-286.

[8] Fang jian-cheng, Zheng Shi-qiang. AMB vibration control for structural resonance of double-gimbal control moment gyro with high-speed magnetically suspended rotor [J]. IEEE-ASME Transactions on Mechatronics, 2013, 18(1): 32-43.

[9] Shafiei N, Kazemi M, Ghadiri M. Comparison of modeling of the rotating tapered axially functionally graded Timoshenko and Euler-Bernoulli microbeams[J]. Physica E: Low-dimensional Systems and Nanostructures, 2016, (83): 74-87.

[10] 萬文軍, 李軍. 一種基于LCR發(fā)散振蕩響應的控制系統(tǒng)頻率特性辨識方法[J]. 電機與控制學報, 2014, 18(11): 113-119. Wan Wen-jun, Li Jun. A method of frequency domain identification for control system based on process response in LCR diverging oscillation[J]. Electric Machines and Control, 2014, 18(11): 113-119.

[11] 朱萌, 曹國武, 張志偉, 等. 基于LEVY法的氣動舵機系統(tǒng)辨識[J]. 彈箭與制導學報, 2011, 31(6): 69-72. Zhu Meng, Cao Guo-wu, Zhang Zhi-wei, et al. The system identification of pneumatic actuator based on levy method[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2011, 31(6): 69-72.

[12] 范達, 范春石, 賀楊, 等. 磁化懸浮-感應驅(qū)動型動量球系統(tǒng)[J]. 中國慣性技術(shù)學報, 2016, 24(4): 524-530. Fan Da, Fan Chun-shi, He Yang, et al. Momentum sphere system of magnetization suspension and induction driving [J]. Journal of Chinese Inertial Technology, 2016, 24(4): 524-530.

Compound identification method of magnetic flexible rotor model

ZHANG Hua-qiang1, WANG Ying-guang2, ZHAO Xue-tao1
(1. School of Mechanical Engineering, Shandong University of Technology, Zibo 255049, China;
2. Beijing Institute of Control Engineering, Beijing 100094, China)

Aiming at the problem of low accuracy of the pure finite element analysis model for the magnetic flexible rotor, a new compound identification method of magnetic flexible rotor model is proposed. First, a finite element method is used to divide the magnetic flexible rotor into a plurality of Timoshenko beam elements. Then, a magnetic flexible rotor system model is built, the key parameters of the flexible rotor such as critical frequency, damping, stiffness and vibration mode are analyzed, and then the equivalent reducedorder mathematical model is obtained. Based on this, the dynamic characteristics of the flexible rotor system are analyzed using a robust adaptive method. At last, a correction LEVY method is used to identify the magnetic flexible rotor system by analyzing the dynamic characteristics data of the rotor, and the equivalent reduced-order mathematical model obtained by the finite element method is calibrated. Experimental results show that this method can obtain a more accurate system model of the magnetic flexible rotor.

flexible rotor; finite element; robust adaptive; correction LEVY method; compound identification

V414

A

1005-6734(2017)02-0249-07

10.13695/j.cnki.12-1222/o3.2017.02.021

2017-01-011;

2017-03-24

山東省自然科學基金(ZR2015FL012);國家自然科學基金(51605031);山東省自然科學基金教育廳聯(lián)合專項(ZR2014JL027)

張華強(1982—),男,博士,講師,從事檢測與導航技術(shù)研究。E-mail: huaqiang.zhang@163.com

猜你喜歡
模態(tài)有限元分析
隱蔽失效適航要求符合性驗證分析
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動化發(fā)展趨勢分析
國內(nèi)多模態(tài)教學研究回顧與展望
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
磨削淬硬殘余應力的有限元分析
由單個模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
巨型總段吊裝中的有限元方法應用
船海工程(2013年6期)2013-03-11 18:57:27
主站蜘蛛池模板: 夜夜高潮夜夜爽国产伦精品| 欧美日在线观看| 亚洲动漫h| 国产国拍精品视频免费看| 亚洲人成在线免费观看| 不卡网亚洲无码| AV在线麻免费观看网站| 国产成人乱码一区二区三区在线| 亚洲无码在线午夜电影| 成人福利视频网| 欧美中文字幕在线播放| 无码人中文字幕| 国产丰满大乳无码免费播放| 久草视频一区| 男人天堂亚洲天堂| 午夜激情婷婷| 国产又大又粗又猛又爽的视频| 午夜成人在线视频| 日本成人一区| 成人精品亚洲| 亚洲最黄视频| 国产亚洲欧美日韩在线观看一区二区 | 综合色区亚洲熟妇在线| 婷婷综合亚洲| 亚洲av成人无码网站在线观看| 精品视频一区二区三区在线播| 88av在线| 国产精品分类视频分类一区| 免费观看无遮挡www的小视频| 欧美福利在线| 国产午夜精品一区二区三| 国产精品天干天干在线观看| 久久亚洲国产最新网站| 91视频精品| 99热这里只有精品国产99| 亚洲精品午夜无码电影网| 国产成人AV综合久久| 色综合天天视频在线观看| 欧美人在线一区二区三区| 亚洲AV无码一区二区三区牲色| 国产在线精彩视频论坛| 国产另类视频| 成年女人a毛片免费视频| 国产精品久久久久无码网站| 91精品国产一区自在线拍| 国产一级毛片网站| a毛片基地免费大全| 免费 国产 无码久久久| 午夜视频在线观看免费网站| 欧美成人免费午夜全| 精品国产欧美精品v| 噜噜噜综合亚洲| 美女扒开下面流白浆在线试听| 免费一级无码在线网站 | 成人亚洲天堂| 欧美精品高清| 91娇喘视频| 日本三区视频| 激情国产精品一区| 国产成人精品第一区二区| 亚洲黄网在线| 国产三级a| 欧美 亚洲 日韩 国产| 欧美精品啪啪| 免费观看男人免费桶女人视频| 麻豆AV网站免费进入| 久久久久亚洲av成人网人人软件| 国内精品九九久久久精品 | 99九九成人免费视频精品| 热99re99首页精品亚洲五月天| 日韩区欧美区| 成人福利在线看| 欧美乱妇高清无乱码免费| 亚洲人成影视在线观看| 在线播放国产一区| 久久视精品| 久久久噜噜噜| 伊人久久大香线蕉成人综合网| 五月天福利视频| 日韩无码视频专区| 国产又色又刺激高潮免费看| 亚洲一区二区三区香蕉|