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

整體平動(dòng)自由結(jié)構(gòu)載荷時(shí)域識(shí)別技術(shù)研究

2016-04-27 02:14:35馬慶鎮(zhèn)韋冰峰
振動(dòng)與沖擊 2016年6期

彭 凡, 馬慶鎮(zhèn), 肖 健, 韋冰峰, 劉 杰

(1.湖南大學(xué) 機(jī)械與運(yùn)載工程學(xué)院,長(zhǎng)沙 410082; 2.北京強(qiáng)度與環(huán)境研究所,北京 100076)

?

整體平動(dòng)自由結(jié)構(gòu)載荷時(shí)域識(shí)別技術(shù)研究

彭凡1, 馬慶鎮(zhèn)1, 肖健2, 韋冰峰2, 劉杰1

(1.湖南大學(xué) 機(jī)械與運(yùn)載工程學(xué)院,長(zhǎng)沙410082; 2.北京強(qiáng)度與環(huán)境研究所,北京100076)

摘要:研究整體平動(dòng)自由結(jié)構(gòu)載荷識(shí)別的Green函數(shù)法及應(yīng)用途徑。建立測(cè)點(diǎn)絕對(duì)運(yùn)動(dòng)加速度與動(dòng)態(tài)激勵(lì)力卷積關(guān)系,采用截?cái)嗥娈愔捣纸獾恼齽t化方法求解反卷積問題。以受軸向激勵(lì)自由梁為對(duì)象,考慮不變及變截面模量?jī)煞N情形,討論整體平動(dòng)及變形運(yùn)動(dòng)在核函數(shù)中的構(gòu)成比例與核函數(shù)矩陣病態(tài)特性關(guān)系,分析結(jié)構(gòu)整體剛度及測(cè)點(diǎn)部位局部剛度對(duì)載荷反演精度的影響機(jī)制及測(cè)點(diǎn)布置需考慮的剛度因素,提高正則化方法計(jì)算穩(wěn)健性;分析復(fù)雜的組合薄壁結(jié)構(gòu),以一點(diǎn)響應(yīng)反求其在牽制釋放實(shí)驗(yàn)中所受沖擊載荷進(jìn)行結(jié)構(gòu)瞬態(tài)響應(yīng)有限元仿真,比較多點(diǎn)仿真與實(shí)測(cè)結(jié)果,驗(yàn)證反求方法的有效性。

關(guān)鍵詞:自由結(jié)構(gòu);載荷識(shí)別;絕對(duì)運(yùn)動(dòng);格林函數(shù);正則化

航天飛行器為自由結(jié)構(gòu)系統(tǒng),受動(dòng)態(tài)載荷作用時(shí)會(huì)產(chǎn)生整體平動(dòng)改變及結(jié)構(gòu)振動(dòng)。較準(zhǔn)確確定此類系統(tǒng)的動(dòng)力學(xué)環(huán)境,實(shí)現(xiàn)減振及結(jié)構(gòu)優(yōu)化設(shè)計(jì),仍需深入研究,利用遙測(cè)響應(yīng)數(shù)據(jù)進(jìn)行動(dòng)載荷時(shí)域識(shí)別具有重要應(yīng)用價(jià)值。針對(duì)整體平動(dòng)被完全消除的約束結(jié)構(gòu),已有多種載荷時(shí)域識(shí)別方法[1],如通過模態(tài)疊加將問題變換到模態(tài)空間中求解的模態(tài)力法[2-4]及基于線性疊加原理將響應(yīng)與動(dòng)態(tài)力關(guān)系以卷積表示[5-6]由反卷積獲得動(dòng)載荷法。劉杰等[7]基于正則化方法求解反卷積方程,實(shí)現(xiàn)結(jié)構(gòu)多源動(dòng)態(tài)載荷反求。王曉軍等[8]在反卷積基礎(chǔ)上研究動(dòng)載荷區(qū)間識(shí)別。而對(duì)航天飛行器一類自由結(jié)構(gòu),動(dòng)載荷識(shí)別方法尚少見。Kreitinger等[9]提出計(jì)權(quán)加速度法(SWAT),由加速度測(cè)量值乘以計(jì)權(quán)系數(shù)獲得等效合力,應(yīng)用難點(diǎn)在于如何確定計(jì)權(quán)系數(shù)。朱斯巖等[10]將模態(tài)矢量按測(cè)點(diǎn)自由度歸一化,構(gòu)造運(yùn)載火箭動(dòng)態(tài)載荷識(shí)別模式,由計(jì)算所得加速度及測(cè)量值在頻域的差異修改識(shí)別的載荷頻域成份,通過反變換獲取時(shí)域載荷。毛玉明等[11]通過模態(tài)分解將響應(yīng)輸入的彈性運(yùn)動(dòng)與整體平動(dòng)分離,在模態(tài)空間中反求彈性運(yùn)動(dòng)及整體平動(dòng)模態(tài)力。

本文研究Green核函數(shù)法用于整體平動(dòng)自由結(jié)構(gòu)動(dòng)載荷反演的技術(shù)途徑。研究軸向平動(dòng)自由梁的載荷反求問題,分析結(jié)構(gòu)整體剛度及測(cè)點(diǎn)部位局部剛度等因素對(duì)核函數(shù)矩陣病態(tài)性影響特征,指導(dǎo)測(cè)點(diǎn)布置;利用測(cè)點(diǎn)實(shí)測(cè)響應(yīng)反求平動(dòng)組合薄壁結(jié)構(gòu)在牽制釋放實(shí)驗(yàn)中的沖擊載荷,以識(shí)別的載荷為激勵(lì)力計(jì)算結(jié)構(gòu)瞬態(tài)響應(yīng),比較結(jié)構(gòu)多點(diǎn)實(shí)測(cè)結(jié)果與基于載荷反求的響應(yīng)值,驗(yàn)證方法的有效性。

1整體平動(dòng)自由結(jié)構(gòu)的沖擊載荷反求

忽略剛-柔耦合效應(yīng),整體平動(dòng)自由結(jié)構(gòu)動(dòng)力學(xué)方程可表示為

(1)

式中:M,C,K為質(zhì)量、阻尼及剛度矩陣。

利用模態(tài)疊加將系統(tǒng)絕對(duì)運(yùn)動(dòng)分解為整體平動(dòng)與彈性振動(dòng),即

y=Φeqe+Φrqr

(2)

式中:Φr,Φe為整體平動(dòng)及彈性振動(dòng)模態(tài)矩陣;qr,qe為整體平動(dòng)及彈性運(yùn)動(dòng)模態(tài)坐標(biāo)。

將結(jié)構(gòu)在點(diǎn)x所受任意動(dòng)態(tài)載荷歷程在時(shí)域內(nèi)表示為

(3)

式中:δ為作用于結(jié)構(gòu)的單位沖激載荷。

在載荷作用點(diǎn)施加單位沖激載荷,在模態(tài)空間載荷產(chǎn)生的彈性振動(dòng)及整體平動(dòng)響可描述為

(4)

式中:ωi,ξi為彈性振動(dòng)各階頻率及模態(tài)阻尼比;Qe,Qr為對(duì)應(yīng)單位沖激力的彈性振動(dòng)及整體平動(dòng)模態(tài)力。

則載荷作用點(diǎn)到響應(yīng)測(cè)點(diǎn)x的加速度脈沖響應(yīng)函數(shù)為

g(x,t)=ge(x,t)+gr(x,t)

(5)

式中:g(x,t)為Green核函數(shù);ge(x,t),gr(x,t)為整體平動(dòng)、彈性振動(dòng)加速度脈沖響應(yīng)。

需要指出的是,以上模態(tài)分解及轉(zhuǎn)換即為便于闡述核函數(shù)矩陣的奇異特性。實(shí)際應(yīng)用時(shí)核函數(shù)可由有限元計(jì)算直接確定。

(6)

在時(shí)間域內(nèi)用m個(gè)等間隔△t采樣點(diǎn)離散式(6),并簡(jiǎn)單表示為

(7)

浙江某精密儀器實(shí)驗(yàn)室使用普通恒溫恒濕空調(diào),室內(nèi)溫度始終偏低,濕度始終過高,無法達(dá)到設(shè)定的需求,在引入高精度恒溫恒濕智能空調(diào)系統(tǒng)后,通過調(diào)整PID算法參數(shù)和溫濕度偏差閾值等,即使室外為陰雨連綿的天氣,室內(nèi)濕度也可控制在±2%,溫度控制在±0.3℃,建設(shè)了長(zhǎng)期以來所期望的理想實(shí)驗(yàn)室環(huán)境。

對(duì)核函數(shù)矩陣G作奇異值分解,即

G=UDiag(σi)VT

(8)

式中:σi(i=1,2, …,m)為奇異值;U=(u1,u2,…,um),V=(v1,v2,…, vm)為左、右奇異值向量矩陣,滿足UTU =VTV=I,I為單位矩陣。

考慮測(cè)量噪聲,由式(7)、(8)得動(dòng)態(tài)外載荷為

(9)

式(9)表明,由于核函數(shù)矩陣存在小奇異值σi,會(huì)放大測(cè)量誤差甚至導(dǎo)致反求失敗。為此引入正則化算子f(α,σi),小奇異值σi趨于零時(shí)f(α,σi)/σi也趨于零,得實(shí)際載荷穩(wěn)定的近似估計(jì)為

(10)

f(α,σ)=1, (σ2≥α)0, (σ2<α){

(11)

由式(11)可得TSVD算法[12]。采用廣義交叉驗(yàn)證準(zhǔn)則[13](GCV)獲得α,構(gòu)造函數(shù)為

(12)

式中:E=I-G[GTG+αI]-1GT;tr(E)表示E的跡;Ω(α)極小值點(diǎn)即GCV準(zhǔn)則確定的正則化參數(shù)。

由式(10)可見,核函數(shù)矩陣的小奇異值對(duì)正則化反演穩(wěn)健性較關(guān)鍵??疾焓?4)、(5),整體平動(dòng)及振動(dòng)加速度核函數(shù)隨時(shí)間變化特征不同,前者具有瞬時(shí)性且與載荷同步,而后者具有延續(xù)性,各自構(gòu)成的核函數(shù)矩陣特性也不同,整體平動(dòng)gr(x,t)的核函數(shù)矩陣為對(duì)角陣,利于減弱反求的病態(tài)性,故由式(6)構(gòu)成的核函數(shù)矩陣特性與兩種運(yùn)動(dòng)加速度占比有關(guān),需分析剛度因素對(duì)病態(tài)特性的影響。

2結(jié)構(gòu)整體與測(cè)點(diǎn)局部剛度對(duì)反求精度的影響特征

文獻(xiàn)[10]的飛行器簡(jiǎn)化模型-自由梁見圖1,由長(zhǎng)50 m、直徑3 m、壁厚0.005 m的鋁質(zhì)圓筒構(gòu)成。材料彈性模量為70 GPa,泊松比0.33。梁有限元模型共101個(gè)節(jié)點(diǎn)、100個(gè)等長(zhǎng)單元,在各節(jié)點(diǎn)加1 600 kg等效集中質(zhì)量。自由梁所受載荷為F(t)=1.414× 105te-5t。

圖1 自由運(yùn)行飛行器簡(jiǎn)化模型Fig.1 Simplified model of space vehicle in unconstrained state

考慮材料與幾何特性沿軸向分布一致的均勻梁,分別由節(jié)點(diǎn)20、80加速度響應(yīng)反求載荷,在載荷作用點(diǎn)1(圖1)施加單位沖激脈沖力計(jì)算20、80節(jié)點(diǎn)的絕對(duì)加速度響應(yīng)獲得Green核函數(shù),再進(jìn)行F(t)作用的正問題響應(yīng)計(jì)算,將20、80節(jié)點(diǎn)絕對(duì)加速度響應(yīng)迭加30%的隨機(jī)噪聲模擬實(shí)測(cè)響應(yīng)。分別據(jù)節(jié)點(diǎn)20、80響應(yīng)識(shí)別的載荷見圖2。由圖2可見,雖疊加的噪聲水平較高,但識(shí)別結(jié)果與實(shí)際載荷仍吻合良好,且識(shí)別精度受測(cè)點(diǎn)位置影響不明顯。

圖2 針對(duì)均勻梁反求出的載荷Fig.2 Identified loads acting on uniform free beam

考慮非均勻自由梁,將梁(圖1)右半部分彈性模量降低至0.14 GPa,其余參數(shù)不變,左半部分剛度較右半部分大500倍,仍分別由節(jié)點(diǎn)20、80的響應(yīng)反求載荷,節(jié)點(diǎn)20位于剛度大的左半部分,而節(jié)點(diǎn)80則在剛度小的右半部分。在兩節(jié)點(diǎn)響應(yīng)中分別疊加10%、1%的隨機(jī)噪聲,反求結(jié)果見圖3、圖4。由兩圖可見,以80作為響應(yīng)測(cè)點(diǎn)時(shí)1%的噪聲導(dǎo)致較大反求誤差,而以20的響應(yīng)所得反求載荷具有一定精度,但較均勻梁(圖2)結(jié)果差。

圖3 非均勻梁中以節(jié)點(diǎn)20的響應(yīng)反求出的載荷Fig.3 Identified loads by the response of node 20 in non-uniform free beam

圖4 非均勻梁中以節(jié)點(diǎn)80的響應(yīng)反求出的載荷Fig.4 Identified loads by the response of node 80 in non-uniform free beam

圖2~圖4反映的不同識(shí)別精度可由Green核函數(shù)矩陣的奇異值差別解釋。節(jié)點(diǎn)20、80分別作為響應(yīng)測(cè)點(diǎn)在均勻、非均勻梁的位置均相同。圖2對(duì)應(yīng)均勻梁,核函數(shù)矩陣后4個(gè)最小奇異值見表1右側(cè),圖3、圖4對(duì)應(yīng)非均勻梁,核函數(shù)矩陣后4個(gè)最小奇異值見表1左側(cè)。由表1可見,圖2對(duì)應(yīng)的奇異值高于圖3與圖4若干量級(jí),此因均勻梁整體剛度大,核函數(shù)中整體平動(dòng)加速度占比大,核函數(shù)矩陣病態(tài)性明顯減弱,故識(shí)別結(jié)果精度較高。圖4以80為測(cè)點(diǎn),核函數(shù)矩陣最小奇異值較小,故反演結(jié)果對(duì)噪聲高度敏感,因非均勻梁整體剛度較均勻梁小,且測(cè)點(diǎn)處于剛度最小的右半梁,變形運(yùn)動(dòng)加速度在整個(gè)絕對(duì)運(yùn)動(dòng)加速度中占比高,核函數(shù)矩陣病態(tài)性加劇。圖3的20測(cè)點(diǎn)位于非均勻梁局部剛度大的左半部分,較80節(jié)點(diǎn)核函數(shù)矩陣病態(tài)性得以改善,使識(shí)別效果較好。

以上討論表明,用Green函數(shù)法反求自由結(jié)構(gòu)動(dòng)載荷時(shí)測(cè)點(diǎn)布置需在綜合分析基礎(chǔ)上,考慮其位置局部剛度,以提高反演法的穩(wěn)健性及抗干擾能力。

表1 核函數(shù)矩陣4個(gè)最小奇異值

3組合結(jié)構(gòu)牽制釋放載荷反求及驗(yàn)證

考慮的組合結(jié)構(gòu)由聯(lián)接底部的內(nèi)外加筋薄壁構(gòu)件組成,內(nèi)部結(jié)構(gòu)含沿z向的薄壁圓筒、縱隔板及垂直z向橫隔板,結(jié)構(gòu)形狀較復(fù)雜,不同部位剛度差別較大。組合結(jié)構(gòu)牽制釋放實(shí)驗(yàn)描述:將其與工裝固連,工裝上端由多股彈性繩張拉,下部連接觸發(fā)機(jī)構(gòu)接地,實(shí)驗(yàn)開始觸發(fā)機(jī)構(gòu)與工裝下端突然脫鉤,組合結(jié)構(gòu)連同工裝在彈性繩的張力作用下突然啟動(dòng)。力學(xué)模型簡(jiǎn)化見圖5(a),底部受均勻分布的沖擊力作用引發(fā)結(jié)構(gòu)z向自由運(yùn)動(dòng)。對(duì)組合結(jié)構(gòu)用板殼單元、梁?jiǎn)卧?、集中質(zhì)量及耦合單元進(jìn)行有限元網(wǎng)格劃分,見圖5(b)??疾靸?nèi)部結(jié)構(gòu)上、中、下部三點(diǎn)A、B、C,其中C點(diǎn)位于橫隔板上,B點(diǎn)位于z向剛度大的縱向板中部。以B點(diǎn)加速度傳感器測(cè)量的z向數(shù)據(jù)為輸入,反求施加于結(jié)構(gòu)的沖擊載荷,再以識(shí)別的載荷重新加載進(jìn)行動(dòng)力學(xué)正演計(jì)算,獲得整個(gè)結(jié)構(gòu)的瞬態(tài)響應(yīng)。比較A、B、C三點(diǎn)結(jié)果,實(shí)測(cè)的軸向加速度響應(yīng)、基于載荷反求的數(shù)值仿真結(jié)果見圖6~圖8。由三圖看出,基于反求的響應(yīng)與實(shí)測(cè)結(jié)果變化趨勢(shì)相同,吻合較好,雖兩者在某些峰值點(diǎn)存在差別,但工程上仍在可接受范圍內(nèi)??紤]組合結(jié)構(gòu)的整體平動(dòng)簡(jiǎn)化與實(shí)際運(yùn)動(dòng)差異、模型復(fù)雜引起的有限元建模難度大及局部模態(tài)阻尼信息不完備等因素產(chǎn)生的誤差,可認(rèn)為基于反求的數(shù)值仿真能較好預(yù)示實(shí)測(cè)響應(yīng),表明本文反求技術(shù)對(duì)復(fù)雜的整體平動(dòng)自由結(jié)構(gòu)載荷反演有效。

圖5 組合薄壁結(jié)構(gòu)與分析模型Fig.5 Composite thin-walled structure and its analysis model

圖6 組合結(jié)構(gòu)B點(diǎn)軸向加速度的比較Fig.6 The comparison of axial acceleration at point B in composite structure

圖7 組合結(jié)構(gòu)C點(diǎn)軸向加速度的比較Fig.7 The comparison of axial acceleration at point C in composite structure

圖8 組合結(jié)構(gòu)A點(diǎn)軸向加速度的比較Fig.8 The comparison of axial acceleration at point A in composite structure

4結(jié)論

通過基于Green函數(shù)法,以測(cè)點(diǎn)絕對(duì)運(yùn)動(dòng)加速度響應(yīng)為輸入直接反演整體平動(dòng)自由結(jié)構(gòu)的動(dòng)載荷,對(duì)均勻、非均勻自由梁軸向激勵(lì)反求問題分析研究,結(jié)論如下:

(1)因整體平動(dòng)與彈性振動(dòng)組成比例不同會(huì)使核函數(shù)矩陣病態(tài)特性呈巨大差異,為提高正則化方法的計(jì)算穩(wěn)定性及抗干擾能力,宜將測(cè)點(diǎn)布置在剛度較大處。

(2)通過將形成方法被用于組合薄壁結(jié)構(gòu)在牽制釋放時(shí)所受沖擊載荷反求,并對(duì)比結(jié)構(gòu)瞬態(tài)響應(yīng)的實(shí)測(cè)及基于載荷反求的數(shù)值仿真結(jié)果,驗(yàn)證反演技術(shù)的有效性。亦為復(fù)雜平動(dòng)自由結(jié)構(gòu)動(dòng)載荷反演提供實(shí)施新途徑。

參 考 文 獻(xiàn)

[1] Sanchez J, Benaroya H. Review of force reconstruction techniques[J]. Journal of Sound and Vibration,2014,333(14): 2999-3018.

[2] 張運(yùn)良,林皋,王永學(xué),等. 一種改進(jìn)的動(dòng)態(tài)載荷時(shí)域識(shí)別方法[J]. 計(jì)算力學(xué)學(xué)報(bào),2004,21(2):209-215.

ZHANG Yun-liang, LIN Gao, WANG Yong-xue, et al. An improved method of dynamic load identification in time domain[J]. Chinese Journal of Computational Mechanics, 2004, 21(2):209-215.

[3] 蔡元奇,朱以文. 基于逆向?yàn)V波器的動(dòng)態(tài)載荷識(shí)別方法[J]. 振動(dòng)工程學(xué)報(bào),2006,19(2):200-205.

CAI Yuan-qi, ZHU Yi-wen. Time domain identification method for dynamic loads based on inverse direction filter[J].Journal of Vibration Engineering,2006,19(2):200-205.

[4] 張方,朱德懋,張福祥. 動(dòng)載荷識(shí)別的時(shí)間有限元模型理論及其應(yīng)用[J]. 振動(dòng)與沖擊, 1998,17(12):1-5.

ZHANG Fan, ZHU De-mao, ZHANG Fu-xiang. The theory and application of time element model for dynamic load identification[J].Journal of Vibration and Shock,1998,17(2): 1-5.

[5] Inoue H, Harrigan J J, Reid S R. Review of inverse analysis for indirect measurement of impact force[J]. Applied Mechanics Review, 2001, 54(6):503-525.

[6] Liu G R,Ma W B,Han X. An inverse procedure for identification of loads on composite laminates[J]. Composites Part B:Engineering,2002,33(6):425-432.

[7] 韓旭,劉杰,李偉杰,等. 時(shí)域內(nèi)多源動(dòng)態(tài)載荷的一種計(jì)算反求技術(shù)[J]. 力學(xué)學(xué)報(bào), 2009, 41(4): 595-602.

HAN Xu, LIU Jie, LI Wei-jie,et al. A computational inverse technique for reconstruction of multisource loads in time domain[J]. Chinese Journal of Theoretical and Applied Mechanics, 2009, 41(4):595-602.

[8] 王曉軍,楊海峰,邱志平,等.基于Green 函數(shù)的動(dòng)態(tài)載荷區(qū)間識(shí)別方法研究[J]. 固體力學(xué)學(xué)報(bào),2011,32(1):95-101.

WANG Xiao-jun, YANG Hai-feng,QIU Zhi-ping,et al. Research on interval identification method for dynamic loads based on Green’s function[J]. Chinese Journal of Solid Mechanics, 2011, 32(1): 95-101.

[9] Kreitinger T, Luo H L. Force identification from structural response[C]//Proceedings of the 1987 Sem Spring Conference on Experimental Mechanics.Houston,1987:851-855.

[10] 朱斯巖,朱禮文. 運(yùn)載火箭動(dòng)態(tài)載荷識(shí)別研究[J]. 振動(dòng)工程學(xué)報(bào), 2008, 21(2):135-139.

ZHU Si-yan, ZHU Li-wen. Dynamic load identification on lunch vehicle[J]. Journal of Vibration Engineering, 2008, 21(2):135-139.

[11] 毛玉明,郭杏林,趙巖,等.自由-自由運(yùn)行體系動(dòng)態(tài)載荷反演問題研究[J]. 計(jì)算力學(xué)學(xué)報(bào),2010,27(1): 35-39.

MAO Yu-ming, GUO Xing-lin, ZHAO Yan,et al. Study of dynamic force identification for free-free structural system [J]. Chinese Journal of Computational Mechanics, 2010,27(1):35-39.

[12] Hansen P C. Truncated SVD solutions to discrete ill-posed problems with ill-determined numerical rank.[J]. SIAM Journal on Scientific and Statistical Computing,1990,11(3):503-518.

Load identification technique in time domain for free structures with overall translation

PENGFan1,MAQing-zhen1,XIAOJian2,WEIBing-feng2,LIUJie1

(1. College of Mechanical and Mehicle Engineering, Hunan University, Changsha 410082, China;2. Research Institute of Beijing Structure and Environment Engineering, Beijing 100076, China)

Abstract:The green kernel function method to identify dynamic loads on structures with overall translation was studied. The convolution integral relationship between dynamic loads and absolute acceleration at measuring points was derived, and the regularization algorithm of truncated singular value decomposition was applied to solve corresponding de-convolution equations. The free beams with constant and variable stiffness distribution were analyzed. The effects of overall stiffness and local stiffness in the area of measuring point on the minimum singular values were examined. It is concluded that the large stiffness near measuring points can improve the computing robustness and noise resistance of the regularization method. The method proposed was used to identify the shock load on a composed thin-walled structure in hold down and release experiment. The response in numerical simulation based on load identification was compared with the measured result, and the validity and accuracy of the identification technique were verified.

Key words:free structure; load identification; absolute motion; Green function; regularization

中圖分類號(hào):O326; O347.1

文獻(xiàn)標(biāo)志碼:A

DOI:10.13465/j.cnki.jvs.2016.06.016

收稿日期:2014-06-24修改稿收到日期:2015-03-31

第一作者 彭凡 男,博士,教授,1963年9月生

主站蜘蛛池模板: 中文字幕天无码久久精品视频免费 | 青草视频免费在线观看| 美女无遮挡拍拍拍免费视频| 制服丝袜亚洲| 色精品视频| 成人av专区精品无码国产| 国产精品永久在线| 色悠久久综合| 大香网伊人久久综合网2020| 久草中文网| 91丝袜在线观看| 凹凸国产熟女精品视频| 欧美性天天| 亚洲自偷自拍另类小说| 亚洲成a人片77777在线播放| 久久免费精品琪琪| 精品超清无码视频在线观看| 大陆精大陆国产国语精品1024| 亚洲国内精品自在自线官| 色AV色 综合网站| 亚洲最大看欧美片网站地址| 亚洲美女AV免费一区| 亚洲 成人国产| 国产91精品调教在线播放| 国产欧美日韩在线一区| 国产精品亚洲片在线va| 久久男人资源站| 久久精品国产精品一区二区| 国产精品视频导航| 亚洲美女一区| 国产精品熟女亚洲AV麻豆| 亚洲天堂777| 综合久久五月天| 狠狠ⅴ日韩v欧美v天堂| 国产人成在线观看| 蜜桃视频一区| 日韩成人午夜| 国产精品手机在线观看你懂的| 天天色综合4| 久久99国产综合精品1| 欧美日韩国产成人高清视频| 国产精品浪潮Av| 国产后式a一视频| 成人在线第一页| 国产一级特黄aa级特黄裸毛片| 最新国产精品第1页| 色精品视频| 好紧好深好大乳无码中文字幕| 欧美性久久久久| 欧美一区二区精品久久久| 69精品在线观看| 国产xx在线观看| 国内嫩模私拍精品视频| 91精品久久久久久无码人妻| 91视频首页| 亚洲伊人天堂| 中文字幕av一区二区三区欲色| 久久精品国产亚洲麻豆| 激情无码字幕综合| 久久黄色毛片| 日韩欧美视频第一区在线观看| 一区二区三区国产精品视频| 日韩在线2020专区| 一区二区三区高清视频国产女人| 91成人在线观看| 欧美色综合久久| 无码国产伊人| 亚洲精品视频在线观看视频| 爱爱影院18禁免费| 欧洲免费精品视频在线| www欧美在线观看| 欧美综合中文字幕久久| 国产女人18毛片水真多1| 亚洲精品在线91| 免费AV在线播放观看18禁强制| 欧美成人精品一级在线观看| 一区二区三区四区精品视频 | 亚洲成A人V欧美综合| 少妇精品在线| 91蝌蚪视频在线观看| 四虎AV麻豆| 成人国产一区二区三区|