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

附加自由阻尼梁高頻響應的能量有限元方法模型

2016-01-15 03:13:00孔祥杰,陳花玲,祝丹暉
振動與沖擊 2015年17期

第一作者孔祥杰男,博士生,1984年7月生

通信作者陳花玲女,教授,博士生導師,1954年8月生

附加自由阻尼梁高頻響應的能量有限元方法模型

孔祥杰1,2, 陳花玲1,2,祝丹暉1,2,張文博1,2

(1.西安交通大學機械工程學院,西安710049; 2.西安交通大學機械結構強度與振動國家重點實驗室,西安710049)

摘要:能量有限元方法(EFEA)是一種預示結構高頻響應的新方法。為了利用能量有限元方法準確的預示附加阻尼結構的高頻響應,將復剛度法與能量有限元的相關理論相結合,同時對現有的能量有限元推導思路進行了大結構阻尼條件下的修正,推導了大阻尼工況下,附加自由阻尼梁結構高頻振動響應的能量密度控制方程。同時,通過對阻尼處理交界面處能量轉移關系的分析,建立了經過局部附加阻尼處理的梁結構高頻響應的能量有限元模型。通過與各自工況下模態解析解的對比,所建立的能量有限元方法模型可以準確的預示大結構阻尼工況下附加阻尼梁的高頻響應。

關鍵詞:能量有限元;附加阻尼結構;高頻響應預示

基金項目:國家自然科學基金項目(11172222)

收稿日期:2014-05-23修改稿收到日期:2014-08-19

中圖分類號:TB123

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2015.17.016

Abstract:Energy Finite Element Analysis (EFEA) is a method developed for high-frequency structural response prediction in recent years. To study high-frequency vibration responses of beam structures with free-layer damping (FLD) treatment, the governing equation of energy density for bending vibration of beam structures with FLD treatment under large damping condition was derived based on the model of equivalent complex flexural stiffness and the theory of EFEA with special treatments and modifications for high structural damping. Meanwhile, the EFEA model of beam structures with partial FLD treatment was also built by studying the energy transfer relationship at the interface of damping treatment. Numerical simulations verified the proposed model through comparing with the analytical modal solutions. The results showed that the proposed EFEA model can predict the high-frequency vibration responses of beam structures with FLD treatment under high damping condition accurately.

Energy finite element ananlysis for high-frequency vibration of beams with free-layer damping treatment

KONGXiang-jie1,2,CHENHua-ling1,2,ZHUDan-hui1,2,ZHANGWen-bo1,2(1. School of Mechanical Engineering, Xi’an Jiaotong University, Xi’an 710049, China;2. State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China)

Key words:energy finite element analysis; free layer damping (FLD) treatment; high-frequency response

阻尼減振降噪技術是目前常用的振動與噪聲控制方法之一,作為其中最常用的一種處理手段,附加阻尼處理對薄壁類大尺寸構件的減振降噪效果特別顯著,因此廣泛的應用在汽車的外殼、車身立柱,飛機機身艙壁大尺寸薄壁結構的振動與噪聲控制中。

在結構共性上,此類大尺寸輕薄的附加阻尼結構的振動響應更多表現為高頻響應。目前求解結構高頻響應問題的主要預示方法為有限元法、統計能量法和能量有限元方法。

有限元法是目前結構響應分析中使用最廣泛的方法。但是,在高頻時,為保證運算精度,有限元法所劃分單元的個數需要成倍增加,求解時所需的計算時間和計算成本急劇增加,甚至使得求解不切實際而無法實現。同時,高頻下結構響應對結構的材料特征微小變化的敏感度遠高于低頻,因此,有限元法一般用于低頻分析,無法完全勝任此類附加阻尼結構的高頻響應預示。統計能量法(SEA)是一種較為成熟的高頻聲振響應的預示方法。它以模態理論為基礎,以子系統平均振動響應的能量密度為研究對象,只需用一組線性代數方程便可描述子系統間的響應能量關系,計算效率高。但是,由于統計特征的限制使子系統劃分粗大,僅能給出結構粗略的平均響應,無法預測子系統中局部位置的精確響應及響應隨位置的空間變化,且不能考慮能量密度在子系統內的變化和阻尼的不均勻性,因此,統計能量法無法滿足精確描述附加阻尼結構高頻響應的要求。

能量有限元方法(EFEA)是近年來提出的一種高頻響應預示的新方法,它基于波動的方法以能量形式描述系統的結構響應,可以描述結構響應隨空間變化的情況,結構的幾何及阻尼特征可以得到充分表達。能量有限元方法采用有限元數值求解,只需要很少的網格便可以較為準確的預示結果,克服了有限元法在預示高頻問題時計算成本高計算時間長的缺點及統計能量法無法預測子系統中局部位置的精確響應及響應隨位置的空間變化的不足,是一種非常具有研究價值和發展前景的結構高頻響應預示方法。

能量有限元方法由Nefske等[2]提出,他們以梁結構為研究對象,采用能量密度及能量流描述結構的高頻響應,最終得到了以能量密度為變量的、近似于熱傳導方程的二階微分方程。同時采用有限元對該方程進行離散求解,正式創立了能量有限元方法。隨后,Wohlever,Bouthier,Cho等[3-5]先后建立了桿結構、梁、板及耦合結構等基本結構的能量有限元方法模型。在應用研究方面,密歇根大學Vlahopoulos教授研究團隊的Zhang、Wang、Lee、Yan等[6-9]先后建立了流體接觸結構、復合層合板等結構的能量有限元模型,并應用到艦船、汽車白車身、飛機旋翼等實際結構的高頻響應分析中。此外,游進等[10-11]、張文博等[12]、Kong等[13]分別對隨機能量有限元、熱環境下的能量有限元理論、能量有限元方法的有效性判據等問題進行了研究。

然而,截至目前,有關附加阻尼結構高頻振動響應預示的能量有限元模型的相關研究尚未見報道。現有的能量有限元理論在控制方程推導過程中對振動波的波數及群速度進行了小阻尼近似,當結構阻尼較小時,這種基于小阻尼假設的近似引起的誤差較小,然而,對于附加阻尼結構而言,為達到相應的減振降噪效果,結構的實際阻尼損耗因子一般很大,若仍采取小阻尼假設而直接使用現有的EFEA控制方程進行求解,勢必會帶來更大的誤差[14]。因此,有必要建立大阻尼條件下的能量有限元模型,對這類結構的高頻響應進行精確預示。

本文將復剛度理論與能量有限元相關理論結合起來,以附加阻尼結構中最為常用的附加自由阻尼結構為研究對象,對現有的能量有限元推導思路進行了大結構阻尼工況下的修正,推導了附加自由阻尼梁及局部附加自由阻尼梁高頻響應的能量密度控制方程,并采用有限單元法對控制方程進行離散求解。最后,通過與多個工況下的精確模態解析解的對比,證明了本文所建立的能量有限元方法模型可以準確的預示附加阻尼梁結構的高頻響應。

1附加自由阻尼梁的能量密度控制方程推導

1.1等效彎曲剛度的求取

圖1為附加自由阻尼梁的橫截面示意圖。在梁整個上表面粘貼了一層附加阻尼層,組成典型的附加自由阻尼結構。其中彈性基層(梁本身)與附加阻尼層的厚度分別為H與H1,梁的寬度為b。圖中點劃線代表組合梁的中性面位置。

圖1 附加自由阻尼梁的橫截面示意圖 Fig.1 Cross section of FLD beam

(1)

阻尼效應對振動的影響可以用結構阻尼損耗因子η來表述,對附加阻尼層而言,有

(2)

其中:η1為附加阻尼層的結構阻尼損耗因子。

一般而言,梁的彈性基層即梁本身的材質為金屬,與附加黏彈阻尼層相比,其阻尼系數要遠遠小于附加阻尼層,故彈性基層的阻尼系數可以忽略。因此有

E*=E(1+iη)≈E

(3)

式中:η為梁本身的結構阻尼損耗因子。

(4)

另一方面,附加自由阻尼梁的等效阻尼損耗因子ηeq可表示為

在求得了等效復彎曲剛度Deq和等效阻尼損耗因子ηeq后,附加自由阻尼梁的橫向振動方程便可表示為

(6)

其中:w為梁的橫向位移;m為附加自由阻尼梁的總質量,m=ρS+ρ1S1,ρ與ρ1分別為彈性基層及附加黏彈阻尼層的密度,S與S1分別為彈性基層及附加黏彈阻尼層的橫截面積;F(x-x0)為作用在x=x0位置上的激勵力。

1.2大阻尼條件下能量密度控制方程的推導

忽略掉近場項后,該附加自由阻尼梁的橫向振動方程式(6)的位移w的解可表示為

(7)

式中:A與B為由邊界條件確定的待定系數;k*為復彎曲波數。

聯立式(6)與(7),可得等效復彎曲波數的表達式為

(8)

其中:定義tanφ=ηeq;ω為振動的圓頻率。

對現有傳統的能量有限元理論而言,由于一般系統的結構阻尼較小,為了簡化推導,對復彎曲波數的表達式(8)進行了小阻尼假設基礎上的近似與簡化。將指數函數次冪項中的三角函數在小阻尼假設的基礎上,進行了等效近似。對傳統的能量有限元理論[3]而言,波數可近似表示為

(9)

一般地,定義keq=keq1+jkeq2,則有

(10)

對小阻尼結構而言,上述近似帶來的誤差較小[12]。然而,對于附加阻尼結構而言,由于黏彈阻尼層的存在,整個結構的阻尼一般都遠大于結構本身的阻尼,當系統的結構阻尼顯著增大,對波數進行小阻尼近似帶來的誤差將顯著增加。因此,對附加阻尼結構,為了確保預示精度,更準確的描述系統的響應特性,應該對現有的傳統能量有限元理論進行修正,使用未進行小阻尼近似的波數式(8)來進行能量控制方程的推導。

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

(11)

則等效群速度Cg-eq可以表示為

(12)

橫向振動的附加阻尼梁的能量密度為勢能密度和動能密度之和。時間平均的能量密度可以表示為[3]

式中:W為能量密度;“〈〉”表示時間平均;“*”表示共軛運算。

時間平均的能量流可以表示為[3]

(14)

式中:q為能量強度;Re表示取實部運算。

將位移的遠場解式(7)及其相應的等效復彎曲波數式(8)代入式(13)和(14),并在半個波長上進行如下式的本地空間平均,化簡后最終可得到

(15)

(16)

聯立式(15)和(16),可得能量密度與能量流的關系為

(17)

定義ηeq-EFEA為附加阻尼結構的等效能量有限元阻尼系數,且

ηeq-EFEA=

(18)

對整個振動系統而言,穩態下的能量平衡方程可以表示為

(19)

振動系統的能量耗散關系可以近似表示為

(20)

聯立能量密度與能量流的關系式(17)、能量平衡方程(19)及能量耗散關系式(20),可得到附加自由阻尼梁的EFEA能量密度控制方程為

(21)

式(21)即為大阻尼條件下附加自由阻尼梁高頻振動的能量密度控制方程。與傳統的梁結構的能量密度控制方程[3]相比,式(21)既考慮了附加阻尼層對梁的影響,又采用了完整的波數表達式,無需小阻尼假設,因此能更好的預示附加阻尼結構的高頻響應。

采用有限單元法對能量密度控制方程式(21)進行離散求解。離散后的最終的節點單元線性方程則可表示為

[Ke]{We}={Fe}+{Qe}

(22)

式中:上標e為表示該量為單元量;{We}為節點能量密度向量;[Ke] 為單元系數矩陣;{Fe}為外加輸入功率矩陣;{Qe}為單元邊界節點處內部功率列陣。

2局部附加自由阻尼梁能量模型的建立

上節的推導基于整個梁均進行了自由阻尼處理的情況,在工程實際中,對整個結構均進行附加阻尼處理會大幅提升處理成本。因此,通常僅對彈性基體的某些位置進行局部附加阻尼處理。因此,本節將建立進行了局部附加阻尼處理的梁結構高頻響應的能量有限元模型。

圖2為一個典型的進行了局部附加自由阻尼處理的梁結構,僅在其中部區域粘貼了黏彈阻尼材料,而其它位置未進行阻尼處理。當振動波傳遞到附加阻尼層交界面位置時,由于交界面處的剛度、厚度和橫截面積等參數發生了改變,振動波將產生部分反射與部分傳遞的現象。因此,可以將進行了局部附加自由阻尼處理的梁結構等效為數個具有不同物理參數的耦合梁結構的振動問題來處理。

圖2 進行了局部附加自由阻尼處理的梁結構示意圖 Fig.2 A beam with partial free layer damping treatment

圖3為交界面位置的節點示意圖,交界面節點位置的能量密度W和能量流q可以用傳播出去部分(上標用“t”表示)和反射回來部分(上標用“r”表示)來描述

(23)

qi=qit+qir

(24)

式中:下標i=1表示進行了附加阻尼處理的部分梁,i=2或3表示未進行附加阻尼處理的那部分梁。由于兩側未進行附加處理的那部分梁物理參數相同,故下述推導中用下標2、3或者其中之一(下標2)來統一表示。

圖3 進行了附加阻尼處理和未進行阻尼處理 的交界面位置的節點關系示意圖 Fig.3 The interface between the sector with or without partial free layer damping treatment

對于單個正在傳遞的振動波而言,能量流與能量密度存在如下關系

(25)

如圖3所示,流出交界面節點每一側的凈能量可以分別表示為

(26)

(27)

式中:τij為i部分梁到j部分梁之間的能量傳遞系數;rii為i部分梁的能量反射系數。對于保守耦合的一維梁結構而言,有τij=τji,rii=rjj同時τ+r=1。

將式(25)代入式(26)及(27)可得到

(28)

(29)

聯立式(23)~(29)最終可得到進行了附加阻尼處理和未進行阻尼處理的交界面位置的振動能量轉移關系為

(30)

(31)

另外,對于未進行阻尼處理的部分梁的高頻響應,可用傳統的梁結構彎曲振動的EFEA能量密度控制方程來表示

(32)

參考文獻能量傳遞系數與能量反射系數通過耦合處理交界面位置的位移、剪力彎矩等連續性條件求取,具體方法可[5]。

將進行了附加自由阻尼處理的梁部分的控制方程式(21)及未進行阻尼處理的梁部分的控制方程(32)與阻尼處理的交界面位置的振動能量轉移關系式(30)和(32)聯立,即可求取局部附加阻尼結構的下高頻響應。

3附加自由阻尼梁能量有限元模型的驗證

為了驗證本文推導的能量控制方程,本節將分別利用一個進行了完全附加自由阻尼處理的簡支梁和一個對局部位置進行了附加自由阻尼處理的簡支梁進行相關對比驗證。EFEA的計算結果將與相應的“精確”模態解析解進行對照驗證。為了保證計算結果的收斂性和精度,所有的模態解析解均截取前5000階模態的計算結果。

3.1附加自由阻尼梁的驗證

首先進行完全附加自由阻尼處理的簡支梁的能量有限元模型的驗證。梁的示意圖見圖4,簡諧激勵力作用在梁的中間位置。梁及附加黏彈阻尼層的相關參數如表1所示。

圖4 附加自由阻尼梁示意圖 Fig.4 A beam with free layer damping treatment

參數名稱參數值梁尺寸L×b(長×寬)/m1×0.01梁本身材料密度ρ/(kg·m-3)2700梁本身彈性模量E/Pa7.1×1010梁本身厚度H/mm2梁本身結構阻尼損耗因子η0.01附加阻尼層材料密度ρ1/(kg·m-3)980附加阻尼層結構阻尼損耗因子η10.8附加阻尼層彈性模量E/Pa3.3×108附加阻尼層厚度H1/mm2或4(按工況不同)

需要說明的是,在結構高頻響應預示領域,高頻的定義是一個相對的概念,并不是一個確定性頻率范圍,而是一個與結構材料及尺寸參數及激勵頻率有關的相對概念。一般認為,當分析頻率下對應的振動波長遠大于結構尺寸時,即可認為該振動頻率下的預示問題為高頻問題[15-16]。本部分將選取激勵頻率為2 000 Hz和8 000 Hz兩種頻率工況進行相關驗證,這兩種工況下各自對應的振動波波長(0.13 m及0.06 m)均遠小于梁的幾何長度,屬于高頻響應問題。

首先計算了激勵頻率為2 000 Hz,附加阻尼層厚度為4 mm時附加自由阻尼梁沿x方向的振動能量密度,并與模態解析解進行了對比,計算結果見圖5。圖中能量密度的參考值為10-12J/m(下同)。由圖中可以看出,EFEA的預示結果與模態解析解之間有著良好的一致性。EFEA的計算結果較為精確的反映了模態解析解的均值。而在激勵點附近,EFEA的預示結果與模態解析解之間存在著一定的誤差,其原因是在能量密度控制方程的推導中,忽略了代表漸逝波的近場項。

圖5 激勵頻率為2 000 Hz,附加阻尼層厚度為4 mm時, 沿著x方向上梁的能量密度分布圖 Fig.5The energy density distribution along the x-axis across the excitation point when frequency centered on 2 000 Hz with the thickness of damping layer equals 4 mm

圖6 激勵頻率為8 000 Hz,附加阻尼層厚度為4 mm時, 沿著x方向上梁的能量密度分布圖 Fig.6The energy density distribution along the x-axis across the excitation point when frequency centered on 8 000 Hz with the thickness of damping layer equals 4 mm

保持其他參數不變,將激勵頻率增加到8 000 Hz,EFEA與模態解析解的預示結果見圖6。可以看到,EFEA的預示結果依然與模態解析解吻合。并且,當激勵頻率增加后,振動能量的衰減更為明顯,與2 000 Hz相比,由激勵點到邊界位置共0.5 m的長度上,能量衰減量增加了約30 dB,因此,附加自由阻尼對高頻振動的減振效果更加的明顯。

若附加阻尼層厚度發生了變化,見圖7,當附加阻尼層的厚度由4 mm減小到2 mm時,整個梁結構的等效阻尼損耗因子也相應降低,附加阻尼層對振動能量的衰減作用變弱,與4 mm厚度時的預示結果相比,由激勵點到邊界位置共0.5 m的長度上,能量衰減量降低了將近40 dB。EFEA的預示結果準確的反映了這一現象。

圖7 激勵頻率為8 000 Hz,附加阻尼層厚度為2 mm時, 沿著x方向上梁的能量密度分布圖 Fig.7 The energy density distribution along the x-axis across the excitation point when frequency centered on 8 000 Hz with the thickness of damping layer equals 2 mm

圖8 激勵頻率為50 000 Hz,附加阻尼層厚度為2 mm時, 沿著x方向上梁的能量密度分布圖 Fig.7The energy density distribution along the x-axis across the excitation point when frequency centered on 50 000 Hz with the thickness of damping layer equals 2 mm

為了驗證本章所給出的能量有限元模型在更高的激勵頻率下對附加自由阻尼梁的預示情況,圖8計算了激勵頻率為50000 Hz、附加阻尼層的厚度為2 mm時,梁結構沿x方向的振動能量密度分布情況。由圖5~9可知,當激勵頻率增加到50 000 Hz后,振動能量更加快速地衰減,而本文所給出的EFEA模型的預示結果依然與模態解析解吻合,可以準確地預示更高頻率下的附加自由阻尼梁的響應。

3.2局部附加自由阻尼梁的驗證

下面對進行了局部附加自由阻尼處理的簡支梁結構進行計算和對比驗證。在圖2所示的簡支梁結構中,梁的中部區域(0.3 m~0.7 m的區段)粘貼了黏彈阻尼層,梁本身及附加黏彈阻尼層的相關參數均如表1所示,其中附加阻尼層厚度為4 mm。

簡諧激勵力作用在梁的中間位置。激勵頻率為8 000 Hz。分別采用EFEA的能量密度控制方程和模態解析解對振動能量響應進行計算,結果見圖9。

圖9 激勵頻率為8 000 Hz,局部附加自由 阻尼梁沿x方向上能量密度分布圖 Fig.9The energy density distribution along the x-axis across the excitation point when frequency centered on 8 000 Hz of the beam with partial FLD treatment

從圖中可以看出,對激勵點位置進行了附加阻尼處理后,振動能量發生了明顯的衰減,結構整體的振動大幅度降低,充分地顯示了恰當的局部附加阻尼處理能夠帶來對高頻結構振動的顯著減振效果。同樣的,EFEA的預示結果較好的反映了進行了附加阻尼處理的區段和未進行附加阻尼處理的區段的“精確”模態解析解的平滑處理后的變化趨勢。

通過以上的算例可以看出,本文所給出的能量有限元模型可以較為精確的預示附加自由阻尼結構梁的高頻響應。當結構整體的阻尼較大時,預示結果依然具有較高的準確度。

4結論

本文將等效復剛度理與能量有限方法的相關理論相結合,得到了附加自由阻尼梁及局部附加自由阻尼梁高頻響應的能量有限元模型。通過與多個工況下的精確模態解析解的對比,對模型進行了驗證。對推導出來的控制方程模型討論如下:

(1)考慮到附加自由阻尼結構的實際結構阻尼較大的實際,在推導的過程中,對現有能量有限元方法的推導思路進行了大阻尼條件下的修正,采用了完整的波數及群速度進行推導,因而可以更準確的預示大阻尼條件下的響應情況,相關驗證證明了該點。

(2)對進行了局部附加自由阻尼處理的梁結構而言,本文將其等價為數個具有不同物理參數的耦合梁結構的振動問題來處理,相關驗證也證明了這種處理方法是可行的。

[1]Oberst H, Becker G, Frankenfeld K. über die D?mpfung der Biegeschwingungen dünner Bleche durch fest haftende Bel?ge II[J]. Acta Acustica United with Acustica,1954,4(1):433-444.

[2]Nefske D, Sung S. Power flow finite element analysis of dynamic systems-Basic theory and application to beams[J]. ASME, Transactions, Journal of Vibration, Acoustics, Stress, and Reliability in Design, 1989, 111: 94-100.

[3]Wohlever J, Bernhard R J. Mechanical energy flow models of rods and beams[J]. Journal of Sound and Vibration, 1992, 153 (1): 1-19.

[4]Bouthier O M, Bernhard R J. Models of space averaged energetics of plates[J]. AIAA Journal,1992,30(3): 616-623.

[5]Cho P E. Energy flow analysis of coupled structures[D]. United States: Purdue University, 1993.

[6]Zhang W G, Wang A M, Vlahopoulos N, et al. High-frequency vibration analysis of thin elastic plates under heavy fluid loading by an energy finite element formulation[J]. Journal of Sound and Vibration, 2003, 263(1): 21-46.

[7]Wang A, Vlahopoulos N, Buehrle R, et al. Energy finite-element analysis of the NASA aluminum test-bed cylinder[J]. The Journal of the Acoustical Society of America, 2006, 119(5): 3389-3389.

[8]Lee S, Vlahopoulos N, Waas A. Analysis of wave propagation in a thin composite cylinder with periodic axial and ring stiffeners using periodic structure theory[J]. Journal of Sound and Vibration,2010,329(16):3304-3318.

[9]Yan X Y. Energy finite element analysis developments for high frequency vibration analysis of composite structures[D]. United States: The University of Michigan, 2008.

[10]游進, 李鴻光, 孟光. 耦合板結構隨機能量有限元分析[J]. 振動與沖擊, 2009, 28(11): 43-46.

YOU Jin, LI Hong-guang, MENG Guang.Random energy finite element analysis of coupled plate structures[J]. Journal of Vibration and Shock, 2009, 28(11): 43-46.

[11]You J, Li H G, Meng G. Validity investigation of random energy flow analysis for beam structures[J]. Shock and Vibration, 2011, 18(1-2): 269-280.

[12]Zhang W B, Chen H L, Zhu D H, et al. The thermal effects on high-frequency vibration of beams using energy flow analysis[J]. Journal of Sound and Vibration, 2014,333(9):2588-2600.

[13]Kong X J, Chen H L, Zhu D H, et al. Study on the validity region of Energy Finite Element Analysis[J]. Journal of Sound and Vibration, 2014, 333(9):2601-2616.

[14]Cremer L, Heckl M, Petersson B A T. Structure-borne sound: structural vibrations and sound radiation at audio frequencies[M]. Springer, 2005.

[15]Bitsie F. The structural-acoustic energy finite-element method and energy boundary-element method[D]. United States: Purdue University, 1996.

[16]Moens I, Vandepitte D, Sas P. A wavelength criterion for the validity of the Energy Finite Element Method for plates[C]// Proceedings of the International Conference on Noise and Vibration Engineering. KU Leuven; 1998, 2001, 2: 599-606.

主站蜘蛛池模板: 99精品热视频这里只有精品7| 美女免费黄网站| 免费人成视网站在线不卡| 国产精品夜夜嗨视频免费视频| 日本a级免费| 国产精品女在线观看| 国产精品第页| 伊人91视频| 久久人搡人人玩人妻精品一| 日本久久网站| 日本道中文字幕久久一区| 国产日韩欧美中文| 99视频在线免费| 香蕉国产精品视频| 国产欧美日韩在线一区| 一级毛片基地| 国产成人乱码一区二区三区在线| 成人一级免费视频| 亚洲国产成人精品无码区性色| 91国内外精品自在线播放| 亚洲有无码中文网| 在线亚洲小视频| 亚洲欧洲自拍拍偷午夜色| 国产爽妇精品| 在线国产欧美| 一本色道久久88综合日韩精品| 日韩在线欧美在线| 国产精品香蕉| 国产精品林美惠子在线观看| 熟妇丰满人妻av无码区| 青青热久麻豆精品视频在线观看| 无码啪啪精品天堂浪潮av| 伊人91视频| 麻豆国产在线不卡一区二区| AV无码无在线观看免费| 久久亚洲精少妇毛片午夜无码| 在线不卡免费视频| 国产精品主播| 亚洲九九视频| 性视频久久| 国产十八禁在线观看免费| 国产无遮挡猛进猛出免费软件| 一级一级特黄女人精品毛片| 精品久久香蕉国产线看观看gif| 欧美色视频网站| 亚洲精品国产首次亮相| 好久久免费视频高清| 这里只有精品免费视频| 亚洲综合欧美在线一区在线播放| 国产精品制服| 日韩无码视频播放| 在线欧美国产| 男人的天堂久久精品激情| 午夜日b视频| 激情亚洲天堂| 亚洲精品在线91| 欧美亚洲欧美区| 国产免费怡红院视频| 99视频只有精品| 黄色网址免费在线| 青青草a国产免费观看| 国产又爽又黄无遮挡免费观看| 亚洲欧美成人综合| 亚洲首页在线观看| 日本久久久久久免费网络| 亚洲欧美人成电影在线观看| 精品国产免费观看| 国产91精品调教在线播放| 亚洲精品波多野结衣| 日韩在线视频网站| 国产粉嫩粉嫩的18在线播放91| 性做久久久久久久免费看| 91偷拍一区| a毛片在线播放| 亚洲欧美日韩色图| 丝袜亚洲综合| 亚洲欧美日韩色图| 自拍中文字幕| 国产青榴视频| 91av成人日本不卡三区| www.av男人.com| 午夜日b视频|