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

多因素分離過程蒙特卡羅仿真平臺設計

2015-07-26 11:13:02李慧通
系統工程與電子技術 2015年9期

李慧通,趙 陽

(哈爾濱工業大學航天學院,黑龍江哈爾濱150001)

多因素分離過程蒙特卡羅仿真平臺設計

李慧通,趙 陽

(哈爾濱工業大學航天學院,黑龍江哈爾濱150001)

火箭和導彈等運載工具在飛行過程中涉及多種分離,分離過程出現問題會導致整個任務的失敗。分離過程中分離體受到多種干擾因素和誤差的影響,而靶場試驗次數有限,無法得到多因素對分離過程的綜合影響。設計了通用型的分離模型,建立包含多種力和干擾的分離載荷庫,基于蒙特卡羅隨機打靶技術,能夠對多種分離過程進行打靶分析,得到多因素干擾的影響情況。對某導彈頭罩分離過程進行了蒙特卡羅打靶仿真,得到了在多種偏差因素下分離體的運動范圍。

仿真;蒙特卡羅方法;導彈頭罩分離

0 引 言

在火箭和導彈發射過程中,將已完成預定工作而且在繼續飛行中無用的部分分離并拋掉,可以改善飛行器后續飛行的質量特性,極大地提高運載能力。發射過程中飛行器一般經歷若干次分離,包括頭罩分離、級間分離、底罩分離、有效載荷分離等,分離方案變化大,分離過程受力復雜,如果針對每種分離情況開發仿真軟件則周期長、費用高。不同分離方案的分離過程不確定因素多,分離體受力和受干擾復雜,干擾對分離過程的綜合影響規律難以確定,很難通過幾次實驗就能找到分離的極限狀態,過去一般都是靠科研人員的經驗來進行估計和估算,這對分離機構的設計以及分離故障的排除帶來了不便。因此亟待于構建一個通用型的分離仿真平臺,能夠針對現階段兩體縱向分離過程進行仿真,并且針對多種干擾因素進行打靶仿真,供科研設計人員對仿真結果進行統計學分析。

目前,國內外學者對偏差對分離的影響進行了一系列研究[18],分離過程中,各種影響因素耦合程度高,無法通過單因素仿真實驗確切得到多種影響因素下分離的極端情況和分離體運動范圍,通過蒙特卡羅隨機多次打靶則能較好地得到分離體運動分布范圍。文獻[9]用蒙特卡羅打靶方法對大氣層內高超音速級間分離進行了研究,分析了隨機偏差對分離運動的影響。文獻[10]對微型導彈命中精度進行了蒙特卡羅打靶仿真研究,分析了多種干擾因素對導彈命中精度的影響。

對于飛行器動力學仿真和分離仿真平臺的開發國內外學者也進行了較多研究[11-17],文獻[18]對火箭飛行過程的運動和受力進行仿真分析中,得到火箭飛行過程中位置、速度、加速度和火箭各子結構相互間的作用力,以及級間力在火箭飛行過程中的變化規律。文獻[19]則對運載火箭在運輸和發射過程中所受到的載荷環境進行了研究,開發了具有一定通用性的大型運載火箭載荷仿真軟件。

現階段所設計的分離仿真軟件對于不同分離方案的通用性不強,對于多種偏差因素考慮不夠全面。本文基于蒙特卡羅打靶技術設計了高通用性的分離仿真軟件,能夠對多種分離方案進行仿真分析。本文最后對頭罩分離過程多種影響因素進行隨機蒙特卡羅打靶仿真,確定需要打靶的某些參數及數值上下限后,通過多次隨機打靶得到頭罩分離極限分布情況和打靶的數據結果,供科研人員參考分析。

1 分離動力學模型

通用的兩體縱向分離仿真平臺必須保證實現對多種分離方案的仿真支持,包括底罩分離、頭罩分離、級間分離、有效載荷分離等。這就要求分離模型具有通用性,并且需要通用的力庫供模型調用。

1.1 坐標系設計

分離過程中涉及到的坐標系包括慣性坐標系、彈體坐標系、速度坐標系和質心坐標系。

慣性坐標系(O-XYZ):以發射點為坐標原點O,位于大地或者大海表面,Y軸方向指向正上方,X軸為水平方向,指向與彈體坐標系X軸成銳角,Z軸與X軸、Y軸構成右手坐標系。慣性系主要用于輸入分離初始運動狀態。

彈體坐標系(Ot-XtYtZt):原點Ot設為飛行器頂端理論頂點處,理論頂點為整流罩頂端椎體的切線形成的曲面在空間的交匯點。Xt軸沿彈體縱軸指向彈頭的正方向,也就是導彈的對稱軸,Yt軸在導彈主對稱平面內,即在發射的瞬時與慣性坐標系XOY平面重合,Zt軸與Xt軸、Yt軸構成右手坐標系。彈體坐標系原點位置不會隨燃料的消耗而改變,主要用于分離初始受力參數的輸入。

質心坐標系(Oc-XcYcZc):原點為分離體的質心,Xc軸平行于彈體的對稱軸,Yc軸指向在分離前與彈體坐標系Yt軸的方向相同,Zc軸與Xc軸、Yc軸構成右手坐標系。主要用于分離后分離體運動計算。

速度坐標系(Ov-XvYvZv):原點位于分離體質心,Xv軸沿分離體的飛行速度方向,Yv在分離體縱對稱平面內,垂直于Xv軸,當速度矢量為水平時指向上方,Zv軸與Xv軸、Yv軸構成右手坐標系。主要用于分離氣動力的計算。

1.2 通用模型設計

設分離體所受外載荷為F(Fx、Fy、Fz),外力矩為M(Mx、My、Mz),分離體質量為m,轉動慣量與慣量積為J(Jx,Jy,Jz,Jxy,Jxz,Jyz),分離體在慣性坐標系下位移變化為r(rx,ry,rz),轉動角速度為ω(ωx,ωy,ωz)。6自由度歐拉方程表示為

本文采用的歐拉角按照Z-Y-X順序旋轉,按照歐拉角的定義和歐拉角與角速度的相互關系可以得到體角速度(ωxωyωz)與歐拉角速度()關系為

由于分離過程中,受到多種外力及干擾因素的影響,為了通用性地表示多種力和干擾,本軟件平臺設置了力庫供使用者調用,如圖1所示。

圖1 分離力庫設計

力庫包含3種類型的力,包括分離體所受單獨力、分離體相互作用力和附加力。每種類型的力最多可添加100個,能滿足現階段幾乎所有二體縱向分離仿真需求。

(1)分離體所受單獨力

分離體單獨受力包括氣動力、地球引力、主火箭推力和分離火箭推力,這類力僅單獨分離體受力,不考慮其對另外一個分離體的運動的影響。

(2)分離體相互作用力

分離體相互作用力包括導向摩擦力、分插拔脫力、彈簧力和推沖器力,這類力對上下級分離體都有作用力,大小相等,方向相反。

(3)附加力

附加力包括附加恒定力、附加時變力和附加隨距離變化的力。附加恒定力的大小和彈體系下力作用方向在分離過程中不發生變化;附加時變力的大小和方向隨分離時間變化而變化,通過對關鍵點力的值進行一階拉格朗日插值得到分離任意時刻力的大小和方向;附加隨距離變化的力是隨著分離體之間距離變化而變化的力,同樣通過一階插值得到任意分離距離力的大小和方向。對于未知的干擾和分離力,附加力可以很好地進行模擬,比如分離體之間的空氣負壓力就可以簡化為隨分離距離變化的力,這樣對于將來新的分離方案可以進行仿真。

1.3 干擾因素分析

分離過程包括底罩分離、頭罩分離、整流罩分離、級間分離、有效載荷分離等,所涉及到的高度范圍從海拔十幾米到幾百千米,分離環境復雜,需要考慮多種干擾因素。所有的干擾因素和分離體參數都要在打靶過程中考慮。

質量幾何偏差:在飛行器制造過程中會出現質量和轉動慣量偏差,質心位置也可能有一定的偏移,在飛行過程中,隨著燃料消耗,相關質心質量參數也會與理論值出現偏差。

初始運動偏差:分離體分離時刻的初始運動參數和轉動參數可能存在一定的偏差。特別是考慮氣動力時,分離時組合體初始姿態角度誤差對分離體后續運動影響很大。

氣動系數偏差:分離過程流場變化比較復雜,氣動系數具有高度非線性,很難得到準確的氣動系數模型,平臺采用一階插值方法得到氣動系數曲線,并且軟件平臺設置氣動系數偏差百分比輸入接口。

主火箭偏差:主火箭的推力作用點可能存在橫移,推力線方向和導彈縱軸可能存在一定夾角,這對分離過程有較大影響,推力值大小也有一定偏差,需設置主火箭推力偏差量接口,在軟件平臺中進行詳細考慮。

分離火箭偏差:分離火箭推力實際值與理論值存在偏差,開機關機時間也有一定的不確定性,開機時推力值無法立即達到理論值,關機后也會存在后效推力,這些在打靶過程中都需要考慮。

2 仿真平臺設計

2.1 平臺總體結構

平臺需要良好的人際交互界面,方便科研人員使用。平臺界面采用Visio Studio分區架構形式,界面直觀,易于操作和顯示。主界面如圖2所示。

圖2 平臺界面

動力學軟件的總體結構對軟件系統的總體功能和系統的可維護性等方面有著至關重要的作用。分離過程蒙特卡羅打靶仿真平臺采用現階段動力學分析軟件常用的前處理、中間分析以及后處理結果輸出3大模塊劃分形式,如此劃分使平臺功能清晰,易于用戶使用和進行軟件維護。其平臺架構如圖3所示。

圖3 打靶仿真平臺架構圖

用戶通過前處理模塊輸入分離建模相關參數、分離體受力參數、相關參數偏差量以及打靶仿真控制參數,平臺將數據傳輸到打靶處理模塊進行多次循環打靶計算,將得到的結果送到后處理模塊進行結果輸出,根據用戶的命令進行數據分析以及繪圖顯示。

2.2 打靶設計方案

在打靶仿真中,需要生成隨機數,因為現階段無法得到真正的隨機數,一般采用取中法、位移法、乘同余法、混合同余法等獲取偽隨機數代替。本文采用混合同余法生成偽均勻分布的隨機數,該方法對初值的依賴較小,產生的隨機數列的性質較好,穩定性也比較好[20]。

用戶設置隨機數生成范圍邊界m和n,則在m~n的范圍內采用混合同余法生成某隨機數p。假設某參數的上邊界值為aup,下邊界值為adown,則生成的隨機抽樣值為

生成隨機抽樣值后代入分離動力學仿真部分進行仿真求解,并多次重復這個過程。利用蒙特卡羅方法進行分離過程仿真的基本步驟:

步驟1 根據用戶輸入參數,建立經過簡化的分離過程動力學模型;

步驟2 生成均服從均勻分布的隨機抽樣值,這樣得到的結果更為保守;

步驟3 將抽樣值加載到分離動力學仿真模型,并進行仿真計算;

步驟4 重復進行步驟2和步驟3,多次進行仿真,即可獲得分離運動過程的子樣集;

步驟5 對多次隨機仿真結果進行分析和輸出。

根據以上基本步驟可以得到打靶仿真流程如圖4所示。打靶仿真流程可分為前處理、打靶計算和后處理3部分,對應平臺架構的前處理模塊、打靶處理模塊和后處理模塊。

圖4 打靶仿真流程圖

2.3 平臺描述及實現

多因素分離過程蒙特卡羅仿真平臺可劃分為3大模塊,包括前處理模塊、打靶處理模塊和后處理模塊。

2.3.1 前處理模塊

前處理模塊功能包括仿真管理參數、模型參數配置、干擾參數配置以及打靶參數控制4大部分,如圖5所示。主要功能為輸入仿真參數和打靶控制參數。

圖5 前處理模塊功能劃分

仿真管理參數功能為輸入積分步長和積分時間等。模型參數配置包括質量幾何參數以及初始運動參數,主要功能根據用戶輸入的參數構建分離動力學通用模型。

干擾參數配置包括分離體受力參數、分離體干擾參數以及受力類型配置。主要功能為根據用戶的輸入構建分離體受力和受外界干擾的數學仿真模型;包括氣動力、主火箭推力、分離火箭推力、分插拔脫力、空氣負壓力、推沖器力等,根據用戶的選擇和輸入的參數生成分離體受力模型。

打靶參數控制包括打靶次數頻率、上下限配置和打靶方式配置。主要功能為控制打靶次數和流程,由用戶選定需要打靶的某些參數和限定相關參數的打靶范圍,并進對隨機抽樣值的生成進行選擇。

2.3.2 打靶處理模塊

打靶處理模塊包括模型求解模塊、動力學模型、隨機數生成和打靶流程控制4部分,如圖6所示。動力學模型根據前處理模塊傳遞的數據自動生成,并送入模型求解模塊進行積分求解運算。隨機數生成模塊根據用戶的選擇生成均勻分布隨機數。打靶流程控制監控打靶計算流程,根據需要停止打靶仿真,并將結果數據打包傳遞給后處理模塊。

圖6 打靶處理模塊功能劃分

2.3.3 后處理模塊

后處理模塊包括統計學分析、結果顯示和結果存儲3部分,如圖7所示。主要功能為分析打靶數據和向用戶直觀的展現打靶結果,便于用戶分析。

圖7 后處理模塊功能劃分

統計學分析模塊對用戶選定的需要輸出的參數進行均值、方差和標準差的處理。結果顯示模塊根據用戶的選擇輸出打靶結果的云圖、餅形圖和柱形圖。結果存儲模塊將打靶結果進行存儲,以txt文檔形式保存,便于閱讀,供以后調用分析。

3 頭罩分離算例分析

潛射導彈一般帶有外頭罩保護彈頭以避免海水沖擊的損壞,當導彈發射后距離海面一定高度后,外頭罩工作任務結束,需要將外頭罩拋離導彈,以減少后續飛行過程中的彈體質量。

分離體上面級為頭罩,下面級為彈體,彈體下部安裝主推火箭,頭罩側面安裝4枚分離火箭,頭罩和彈體之間用爆炸螺栓連接,分離體之間有導向機構和控制線纜插頭。由于分離是在低海拔稠密大氣層中分離,所以需要考慮氣動力對分離過程的影響。頭罩分離示意圖如圖8所示。

圖8 頭罩分離示意圖

頭罩主要受氣動力、分離火箭力、空氣負壓力、分插拔脫力和分離機構沖擊力的影響;彈體主要受氣動力、空氣負壓力、分插拔脫力、主火箭力和分離機構的沖擊力影響。

頭罩上安裝4枚分離火箭,以90°夾角沿頭罩圓周均勻分布,推力線與導彈中軸線呈35°夾角。彈體坐標系下分離火箭安裝形式如圖9所示。

圖9 分離火箭安裝示意圖

頭罩在分離階段受力情況為

下面級導彈在分離階段受力為

式中,G為地球引力;R為氣動力;Fi(i=1,2,3,4)為4個分離火箭的推力矢量;Ff為分插拔脫力;Fc為爆炸螺栓沖擊力;Fn為空氣負壓力;Fz為下面級主火箭推力。MR為氣動力矩;Mi(i=1,2,3,4)為4個分離火箭的推力力矩;Mf為分插拔脫力矩;Mc為爆炸螺栓沖擊力矩;Mn為空氣負壓力矩;Mz為下面級主火箭推力矩。

在實際飛行過程中,導彈的飛行狀態比較復雜,受到的影響因素很多,對分離體進行打靶分析時必須首先進行一定得簡化和模型假設。

(1)分離體看作剛體,不考慮其形變對氣動參數的影響;

(2)由于分離時間短暫,不考慮分離時橫風對分離過程的影響;

(3)不考慮分離火箭的噴流對下面級分離體流場的影響;

(4)分離火箭推力假定為恒定力,將火箭推力上升段和下降段的變力和中間的恒定推力在保證總沖量一致的前提下整合,簡化為一恒定推力;

(5)分離機構的沖擊力簡化為一短時間的恒定力;

(6)分離運動仿真不考慮地球自轉和地球曲率的影響。

3.1 計算參數配置

本文研究的火箭頭罩分離由于不同的設計方案,質量特性變化較大,在1 000~1 500 kg范圍內浮動,同時轉動慣量也同時出現變化,質心位置由于制造誤差,也會產生一定范圍內的偏差。頭罩質量幾何偏差如表1所示。

表1 頭罩質量幾何特性表

分離火箭在實際情況下關機時間會出現誤差,不可能在指令發出的精確時間點關機,將關機分離火箭推力下降段簡化為一個短時間的恒定力,這個恒定力和推力下降段沖量一致。分離火箭的推力為30 k N,火箭的開關機時刻表如表2所示。

表2 分離火箭開關機時刻表

3.2 仿真計算結果

500次打靶計算后得到大量結果數據,進行分析并且繪圖輸出,由于本文的篇幅有限,本文列舉部分具有代表性的結果圖。本文取分離后0.6 s時間節點進行分析,結果云圖中黑色球形塊為打靶結果值。

在慣性坐標系下,頭罩X方向、Y方向和Z方向位置分布圖如圖10所示,X方向位置在-2~2 m之間,Y方向高度在1 072~1 080 m之間,Z方向在-1~3 m之間。在Y方向上頭罩的散布范圍比X方向和Z方向要大。

根據打靶結果分析,頭罩和彈體質心的距離與頭罩質量變化相關性較大,如圖11所示,隨著上面級質量的增加,分離距離相應減少,質心距離在12~18 m之間。

質心相對距離與分離火箭3以及分離火箭4的關機時間關系如圖12所示,分離后0.6 s時質心相對距離變化范圍為12~18 m之間。此時分離火箭1和分離火箭2仍在工作。

圖10 頭罩位置打靶結果圖

圖11 質心相對距離與上面級質量關系圖

圖12 質心相對距離與分離火箭工作時間關系圖

圖13 分離體相對角度分布情況

由圖13可知,頭罩分離到0.6 s時,頭罩和下面級導彈的相對俯仰角在-30°~40°之間,相對滾轉角在-15°~15°之間,相對偏航角在-40°~40°之間。

由圖14可知,相對俯仰角度與質心Y方向的位置關聯性較大,隨著質心位置在彈體系下由-0.1 m過渡到0.1 m,相對俯仰角由-30°增加到40°。

圖14 質心位置Y方向與相對俯仰角關系圖

由于本文篇幅有限,無法給出更多的結果圖,通過打靶計算可以得到分離后頭罩和下面級火箭的分布范圍,得到分離結果和偏差量的對應關系,供給科研人員進行分析和驗證。

4 結 論

本文設計的通用分離蒙特卡羅打靶平臺通用性強,能夠仿真多種分離方案,考慮的偏差和干擾因素較多,界面直觀易用。平臺能夠幫助科研人員對多偏差條件下分離過程分離體的運動范圍、潛在風險以及失敗概率進行預測,并且輔助科研人員進行分離機構的設計開發和驗證,模擬多種情況下分離過程,研究導致分離出現問題的原因,找到分離失敗的可能情況。

[1]Singaravelu J,Jeyakumar D,Nageswara R B.Taguchi's approach for reliability and safety assessments in the stage separation process of a multistage launch vehicle[J].Reliability Engineering&System Safety,2009,94(10):1526-1541.

[2]Duprey K E,Saucier E R.Separation systems comparison for ares I launch Vehicle[C]∥Proc.of the 44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit,2008:1-6.

[3]Zhu X,Li H,Yu T,et al.Research on reliability analysis for low-altitude and high-speed payload fairing separation[C]∥Proc.of the International Conference on Quality,Reliability,Risk,Maintenance,and Safety Engineering,2013:90-94.

[4]Singaravelu J,Jeyakumar D,Nageswara R B.Reliability and safety assessments of the satellite separation process of a typical launch vehicle[J].The Journal of Defense Modeling and Simulation:Applications,Methodology,Technology,2012,9(4):369-382.

[5]Oh C S,Sun B C,Park Y K,et al.Payload fairing separation analysis using constraint force equation[C]∥Proc.of the Inter-national Conference on Control,Automation and Systems,2010:1134-1138.

[6]Roshanian J,Talebi M.Monte Carlo simulation of stage separation dynamics of a multistage launch vehicle[J].Applied Mathematics and Mechanics,2008,29(11):1411-1426.

[7]Pamadi B N,Tartabini P V,Toniolo M D,et al.Application of constraint force equation methodology for launch vehicle stage separation[J].Journal of Spacecraft and Rockets,2013,50(1):191-205.

[8]Gusman M R,Barad M F.Aerodynamic database generation for SRB separation from a heavy lift launch vehicle[C]∥Proc.of the 29th AIAA Applied Aerodynamics Conference,2011:1-10.

[9]Jia R Y,Jiang Z Y,Zhang W H.Simulation of off-nominal parameters disturbance of hypersonic vehicle stage separation[J].Journal of Rocket Technology,2012,35(5):578-582.(賈如巖,江振宇,張為華.高超聲速飛行器級間分離偏差干擾仿真[J].固體火箭技術,2012,35(5):578-582.)

[10]Luo Q,Zhang W,Li W.Simulation research of a miniature missile with Monte-Carlo method[J].Flight Dynamics,2013,1(3):265-268.(羅俏,張偉,李偉.微型導彈蒙特卡羅打靶仿真研究[J].飛行力學,2013,1(3):265-268.)

[11]Schneider S A,Chen V W,Pardo-Castellote G,et al.Control-Shell:a software architecture for complex electromechanical systems[J].The International Journal of Robotics Research,1998,17(4):360-380.

[12]Liu X T,Liu L,Song K,Important considerations of building pattern to complicated system[J].Journal of System Simulation,2007,19(13):3073-3075.(劉興堂,劉力,宋坤,等.對復雜系統建模與仿真的幾點重要思考[J].系統仿真學報,2007,19(13):3073-3075.)

[13]Maier M W.System and software architecture reconciliation[J].

Systems Engineering,2006,9(2):146-159.

[14]Garcia J G,Ortega JG,Garcia A S,et al.Robotic software architecture for multisensor fusion system[J].IEEE Trans.on Industrial Electronics,2009,56(3):766-777.

[15]Wu X G,Yang Y J.The design of integrated simulation for underwater weapon system[J].Systems Engineering and Electronics,1996,18(7):61-68.(吳旭光,楊益軍.水下武器系統一體化仿真技術研究與設計[J].系統工程與電子技術,1996,18(7):61-68.)

[16]Bunzel S.AUTOSAR-the standardized software architecture[J].Informatik-Spektrum,2011,34(1):79-83.

[17]Miyachi C.Agile software architecture[J].Association for Computing Machinery(ACM)Special Interest Group on Software Engineering(SIGSOFT)Software Engineering Notes,2011,36(2):1-3.

[18]LüM L.Research on the rocket flight process simulation based on multi-body dynamic method[D].Beijing:Beijing Jiaotong University,2013.(呂明亮.基于多體動力學的火箭飛行仿真研究[D].北京:北京交通大學,2013.)

[19]Wang F.The load calculation of launch vehicle and common software implementation[D].Changsha:National University of Defense Technology,2001.(王鋒.運載火箭載荷計算及通用軟件實現[D].長沙:國防科學技術大學,2001.)

[20]Zheng L,Song Z Y.Algorithms to generate pseudo random numbers and comparison[J].Journal of Hubei University of Technology,2008,23(5):65-68.(鄭列,宋正義.偽隨機數生成算法及比較[J].湖北工業大學學報,2008,23(5):65-68.)

Design of separation process Monte-Carlo simulation platform considering multiple factors

LI Hui-tong,ZHAO Yang
(School of Astronautics,Harbin Institute of Technology,Harbin 150001,China)

Launch vehicles such as rocket and missile always experience a variety of separations in flights,and any problem in the process of separations may result in failure of the whole task.The separation process of detached body is influenced by many kinds of interferences and errors,and the trial number of target range is generally limited,thus it is difficult to precisely determine the influence of multiple factors in the separation process.A general separation model is designed,and a separation load library containing many forces and interferences is established.The simulation platform is on the base of Monte-Carlo method which can be used into simulation analysis for many different kinds of separation process,and then obtain the outcome considering the influence of multiple factors.The Monte-Carlo trajectory simulation for separation process of some missile hood is carried out and the detached body movement range is acquired.

simulation;Monte-Carlo method;separation of missile hood

V 475 文獻標志碼:A DOI:10.3969/j.issn.1001-506X.2015.09.32

李慧通(1988-),男,博士研究生,主要研究方向為飛行器仿真、多體動力學。

E-mail:lihuitongyx@126.com

趙 陽(1968 ),男,教授,博士,博士研究生導師,主要研究方向為飛行器仿真、振動與沖擊。

E-mail:yangzhao@hit.edu.cn

1001-506X(2015)09-2169-07

2014-08-18;

2015-01-04;網絡優先出版日期:2015-03-23。

網絡優先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20150323.1706.003.html

主站蜘蛛池模板: 综合久久五月天| 亚洲精品图区| 一本大道无码日韩精品影视| 成人精品视频一区二区在线| 人妻熟妇日韩AV在线播放| jizz亚洲高清在线观看| 在线观看91香蕉国产免费| 97se亚洲| 亚洲成人在线网| 色呦呦手机在线精品| 91破解版在线亚洲| 日韩精品一区二区深田咏美| 国产成人喷潮在线观看| 免费A级毛片无码免费视频| 久久久久国产精品免费免费不卡| 国产一区二区福利| 制服丝袜 91视频| 伊伊人成亚洲综合人网7777| 亚洲男人天堂网址| 久久青草精品一区二区三区| 欧美日韩一区二区三| 欧美亚洲一二三区| www.亚洲一区| 一本一道波多野结衣av黑人在线| 精品视频在线观看你懂的一区| 91色老久久精品偷偷蜜臀| 国产高潮视频在线观看| 国产激情第一页| 国产精品入口麻豆| 国产经典免费播放视频| 日韩专区第一页| 中文字幕亚洲精品2页| 深爱婷婷激情网| 91精品国产一区| 久久国产精品波多野结衣| 日本不卡在线播放| 欧美日韩国产在线人成app| 日韩无码白| 国产毛片基地| 国产精品一区二区不卡的视频| 亚洲欧美在线综合一区二区三区| h视频在线观看网站| 亚洲视频在线网| 日本欧美一二三区色视频| 成人在线观看不卡| 久久99国产综合精品1| 国产精品自在线拍国产电影| 熟妇无码人妻| 久草视频一区| 米奇精品一区二区三区| 91无码人妻精品一区| www.youjizz.com久久| 亚洲第一区精品日韩在线播放| 91精品国产无线乱码在线| 欧美成人午夜视频免看| 国产精品人成在线播放| 欧美69视频在线| 欧美精品xx| 色成人综合| 亚洲伊人天堂| 老司国产精品视频| 国产午夜不卡| 中国毛片网| 亚洲精品无码av中文字幕| 国产欧美日韩在线一区| 高清色本在线www| 日韩中文字幕亚洲无线码| 动漫精品中文字幕无码| 日本高清有码人妻| 综合人妻久久一区二区精品 | 日本欧美视频在线观看| 婷婷丁香色| 国产黄色片在线看| 国产尤物在线播放| 无套av在线| 无码专区第一页| 韩日午夜在线资源一区二区| 国产免费久久精品99re不卡| 67194在线午夜亚洲| 免费看a毛片| 成人一区在线| 天堂岛国av无码免费无禁网站|