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

基于彈丸激波的彈目偏差測試與反演研究

2017-04-11 06:46:56顧國華張良吳海兵崔遜學俞波
兵工學報 2017年3期

顧國華, 張良, 吳海兵, 崔遜學, 俞波

(陸軍軍官學院, 安徽 合肥 230031)

基于彈丸激波的彈目偏差測試與反演研究

顧國華, 張良, 吳海兵, 崔遜學, 俞波

(陸軍軍官學院, 安徽 合肥 230031)

由于空中超聲速彈丸激波的復雜性和易被干擾性而難以捕捉,其彈目偏差測試是一大難題。為研究基于彈丸激波特性的彈目偏差測試原理和方法,在描述彈目偏差測試場景的基礎上,分析激波形成機理和激波形成前后參數關系,建立彈目偏差測試模型和基于彈群激波環境簡化的線性及非線性干擾模型。綜合采用時間寬度測試原理進行外場試驗與采用激波壓力場累積效應的壓力測試進行仿真反演相結合,結果表明仿真反演與工程測試相符合,為超聲速運動目標聲學定位和跟蹤提供了一種新思路和方法。

兵器科學與技術; 彈道激波; 彈目偏差測試; 聲場反演; 聲壓

0 引言

早在20世紀五六十年代學者們就對空中飛行物產生的激波理論進行了相關探索[1],后來的研究[2-4]集中在對超音速飛行器環境和噪聲影響的評估上。然而,由于激波干擾復雜、厚度難以感知、湍流不確定性等因素,常用數值模擬方法來探討有關激波的空氣動力學問題和飛機設計[5]。

當前,空中彈丸射擊彈目偏差的測試方法主要有光電交匯、圖像判讀、雷達探測、聲學探測、人工觀測等[6]。這些方法各有優點和不足:1)光學交匯法的優點[7]是,實時性強、測量精度高、自動化程度高,不足[8]是視場交匯重疊區不能覆蓋彈著點則無法測量,安裝調試復雜,環境影響因素多;2)圖像(錄像)判讀法的優點[9]在于,測量、存儲方便直觀,缺點[10]在于彈丸來自不同方向對精度影響大,難以實現實時成績評定和報送;3)雷達探測法適用于大口徑彈丸,小口徑彈丸分辨率較差;4)聲學探測法優點[11]在于定位精度高、安全性好,能實現全天候工作,不受光線和能見度影響,缺點是在復雜聲場探測困難。前3種方法,還存在難以工作在彈群環境中,不能辨識重疊(彈丸相距很近)彈丸。

因而聲學探測法具有其獨特的優勢。然而水下聲探測研究成熟,空中聲探測研究少;空中聲標量探測研究較多,空中聲矢量探測研究甚少。為發揮聲學測試方法的優勢,本文結合所在團隊和前人多年在聲學測試研究的基礎上[12-16],進一步研究空中激波環境中超音速彈丸的彈目偏差測試和仿真反演。在研究和分析彈丸彈目偏差解算時,為方便建立理論數學模型,進行了如下假設與簡化:1)將“N”形波傳播圓錐面理想化為2個幾何上的圓錐面,未考慮激波厚度[17-18];2)將有一定幾何尺寸的傳感器陣簡化為一個點;3) 將一傳感器到圓錐面距離近似簡化為到其切平面距離;4) 仿真時將空氣定為理想氣體[19-20](絕熱,傳播均勻,未考慮風速、氣溫等影響);5) 拖靶飛行速度為勻速;6) 其他理想化的因素。這些因素值得繼續深入研究和細化。

1 彈道激波形成

彈道激波是彈丸以超音速飛行時產生的一種非線性大氣聲學現象。如圖1所示,彈丸飛出炮口后,以超音速飛行時,前方的空氣被擠壓來不及向外擴散而形成一定程度的高壓;而彈丸后方的空氣被排開而膨脹,形成較為突然的負壓。大氣中壓縮區(高壓區)和膨脹區(負壓區)結合起來,產生了空氣動力學中的凹角轉折和凸角轉折現象,便在彈頭、彈帶和彈尾部形成近似圓錐形的稠密空氣層,就形成彈道激波,同時,在彈尾底部后面呈現低真空的渦流區。該現象在空氣動力學中稱為激波(也稱聲爆),在彈道學中分別稱為彈頭波、彈帶波和彈尾波,統稱為彈道波[21]。

圖1 彈丸飛行產生的彈道波Fig.1 Ballistic shock wave of flying projectile

空氣的壓縮或膨脹使空氣流波形變陡,其過程實際上氣體的壓力、密度、溫度和速度等空氣參量的微擾動迅速累積而產生突變。顯然,激波具有一定厚度,可用空氣分子的自由行程加以計量,此處略。為更方便分析激波產生前后空氣參數的變化,認為激波是一種不連續面,忽略其厚度[22]。

彈丸產生的激波在空氣中向四周傳播,由于空氣介質的非均勻性(聲阻、氣溫、風),使得非線性因素產生累積效應,激波中的高壓區傳播速度快于聲速,而激波的低壓區傳播則低于聲速。傳播時使得彈丸激波高壓區和低壓區容易形成一個漸進的外形,形似英文字母“N”,所以通常又將其稱為“N”形波。

彈丸激波的形成大體上可由4個過程來描述[11],如圖2所示。彈丸頭部錐體形成的凹角轉折使壓力從常壓p0急速上升至p0+p1,錐柱界面形成的凸角轉折使壓力從p0+p1降至p0,柱面膨脹又使壓力從p0漸降至p0-p2,彈丸底端的界面變化又使壓力從p0-p2急速回升至p0. 圖2中,TF為激波時間寬度,ti為砰發點時刻,τ1、τ2為前后沿上升時間。

圖2 “N”形波Fig.2 “N”shock wave

經進一步研究表明,理想“N”形波前后沿上升時間τ1=τ2,幅值p1=p2.

2 彈目偏差測試場景與激波聲場

2.1 彈目偏差測試場景

當彈丸從靶標附近飛過,其彈道距靶心的最短距離稱為彈目偏差(也稱脫靶量)。其基本場景和主要參數關系如圖3所示。可知理論上彈目標偏差可表示為dcosμ(d為“N”形波傳播距離),而現實中由于炮彈、拖靶、探測器陣都在運動,以及環境、氣象等因素影響,實時情況非常復雜,而給理論模型和工程實踐帶來了難度。

圖3 彈目偏差向量測試原理Fig.3 Test principle of vector miss distance

當彈丸向拖靶發射時,通常瞄向其靶心的未來點。當彈丸運到某一位置P點時,“N”形波前沿4以平行于線d的方向傳播。當彈丸到達3′位置時“N”形波前沿4′到達靶前部探測器5′,此時靶在位置2′,在P點產生的“N”形波將遇見在位置5′的探測器,將P點稱作砰發點。

由于“N”形波傳輸到傳感器期間自身的運動和靶的運動,因而P到“N”形波接收點的距離并不是彈目偏差距離。圖3中用實線表示的拖靶是彈丸在位置3時距其最近時所處的位置,此時彈丸3到靶心T的距離(位置6)dT即為彈目偏差距離。

由于實踐測試的困難性,工程上是通過測量的“N”形波時間寬度來計算dT的。

2.2 運動聲場激波前后關系分析

為進行試驗反演,首先對激波前后參數關系進行建模[22-24]分析,描述彈丸激波前后關系如圖4所示。將超音速彈丸(速度vb)運動轉化為超音速來流狀態。產生激波前,來流為vb,氣流速度沿彈丸切線方向分量vbt,氣流速度沿彈丸法線的方向分量vbn;產生激波時,激波傾角為β;產生激波后,氣流的切線分量為vat,法向分量為van,激波轉折角為ω.

圖4 激波前后參數關系Fig.4 Parameter relations before and after shock wave forming

由克拉珀龍方程得

(1)

式中:k為熱容;ρ表示氣體密度,ρb為激波產生前空氣密度,ρa為激波產生后空氣密度;p為單位體積空氣壓力,pb為激波產生前空氣壓力,pa為激波產生后空氣壓力;m為單位空氣質量;M為氣體摩爾質量;R為氣體常量;vs為考慮熱消耗時激波產生后的速度;c0為理想狀態空氣聲速。

由能量守恒,激波產生前后認為空氣是絕熱的,以及由激波前后的質量流量相等原理,求解(1)式,可得

(2)

(3)

(4)

(5)

式中:Mab為激波產生前氣流的瞬時馬赫數。

根據能量方程及郎金- 雨貢紐關系式,有

(6)

(7)

在空氣中,氣流通過激波時可認為總溫度保持不變,有

(8)

至此,可建立激波產生前后的迭代關系,為解算彈目偏差提供試驗理論依據。

3 彈目偏差測試基本模型

所謂彈目偏差是指火炮在瞄準靶射擊過程中,彈丸飛經靶附近時距靶心的最短距離[12]。為提高脫靶量測試精度和判別射彈散布區域,工程實現是采用了4個聲傳感器Mi(i=1,2,3,4)并組成一個陣列,布置在一個平面上,每個傳感器代表一個象限,方位信息則通過首遇“N”形波的傳感器得到,如圖5所示,圖中vt為拖靶運動方向。

圖5 傳感器布局示意圖Fig.5 Layout diagram of sensors

傳感器M1和M2、M3和M4分別為一組,用于誤差修正。用4個傳感器檢測和計算脫靶距離R,彈目偏差各參量關系如圖6所示。

圖6 彈目偏差各參量關系圖Fig.6 Parameter relations of miss distance

根據sinμ=c0/vp=1/Ma,有

(9)

式中:vp為彈丸速度。

第1節中指出,彈丸激波為一圓錐面,為更好分析激波傳播情況,確定任意時刻激波參數,根據Whitham理論在文獻[16]基礎上建立激波面模型:

(10)

式中:A(r)為激波某區域的面積,r為區域半徑;uf為氣流速度。由(10)式可推導出:

1.1.2 父本來源。父本由皖芝2164 Co60輻射誘變系選,單稈型,一般株高160 cm。葉綠色,白花,一葉三蒴,蒴果四棱,單株蒴果數82.5個。始蒴部位較低,中長果型,結蒴較密,每蒴粒數60.8粒左右,千粒重3.02 g,種皮白色。耐旱性較強,抗枯萎病。正常夏播全生育期90 d左右。2013年參加安徽省區試,平均產量為1 479.0 kg/hm2,比豫芝4號增產0.68%。2014年通過安徽省鑒定。

(11)

(12)

根據速度勢[25]

(13)

求出推遲時間可表示為

(14)

式中:R為有效輻射半徑;Q0為聲源原始強度;Ma0為聲源原始馬赫數;ct為聲源傳播速度;x為接受點處位置,y為聲源在傳播方向所處位置;x1、x2、x3分別為某3個位置接收點。

激波面上有

(15)

式中:ρ0為聲源處空氣原始密度。

因而,超聲速聲源傳播(輻射)距離可表示為(根號前取正號)

(16)

(17)

至此,可解算出任意時刻,激所處的位置極其壓力。

4 彈群激波干擾分析

(18)

式中:U、V為激波在x、y方向上的速度;E為單位質量氣體總能量;u、v為彈丸在x、y方向上的速度;Θx、Θz為x、z方向上張力瞬態值;τxx、τzz、τxz、τzx為黏性應力張量在不同方向上的分量。

假設在擾動情況下,流場矢量q的擾動前為q0,擾動量為q′,則q=q0+q′,擾動的通量包含了非線性的擾動項。在忽略了擾動黏性項后,可對(18)式進行分析,將擾動分為線性和非線性擾動。設ρ0、p0、u0、v0、e0分別是密度、壓強、速度以及總能量的各方向平均量,其對應的擾動量分別為ρ′、p′、u′、v′、e′;E0為總能量,對應的擾動量為E′.

1)當為線性擾動時,有

(19)

2)當為非線性擾動時,有

(20)

由完全氣體的狀態方程(γ為比熱比),E可表示為

5 外場試驗

為驗證模型的有效性和系統功能,在某試驗場進行了野外試驗。工況為:地面立靶2 m×2 m,靶中心高出地面2.4 m,靶面下邊緣離開地面1.4 m. 現場布置示意圖如圖7所示。

圖7 外場試驗布置Fig.7 Field test arrangement

5.1 現場布置

試驗儀器包括,指示器(2傳感器系統和4傳感器系統各一套)、地面站主機、示波器及打印機等。

設置A型炮射擊距離為525 m、B型炮射擊距離為600 m,指示器先后放在距靶心5 m、10 m、

15 m、20 m、25 m、30 m、35 m的位置上,在每一處接收3發能同時測到彈丸存速的原始數據,并測量出M1到靶洞的脫靶距離。

調整指示器支架使傳感器M1、M2、M3和M4在同一平面上,且上面兩傳感器與下面兩傳感所在直線平行。調整指示器使指示器平面與立靶平面在同一平面上。

試驗所采用的彈目偏差的經驗公式[12,15],即

(21)

式中:C1、C2為校準系數。

5.2 試驗數據

I型高炮某試驗場的部分試驗數據見表1.

表1 A型高炮部分試驗數據Tab.1 Some experimental data of antiaircraft gun A

注:指示器采用4傳感器系統,炮目距離為550 m,指示器位于靶平面前側3 m.

采用Levenberg-Marquardt法和通用全局優化法進行擬合,得出:均方根誤差(RMSE)為0.666 346 546,殘差平方和(SSE)為6.660 265 800,相關系數R為0.997 007 629,相關系數之平方R2為0.994 024 213,決定系數(DC)為0.993 966 89,Chi-Square系數為0.207 317 434, F-Statistic統計為2 162.445 651. 結果得出C1=253.902,C2=30.884. 實測試驗數據與擬合的理論數據線曲線與如圖8所示,誤差如表1所示。精度較高,符合系統測試需求。

圖8 外場試驗數據點分布擬合圖Fig.8 Field test data points and their fitted curve

6 仿真反演與分析

為了驗證模型和方法的有效性,對2發超聲速彈丸鄰近干擾情況進行了仿真模擬和計算。初始條件來流為理想空氣,Ma=2.12,靜壓為101 325 Pa,來流靜溫為103.2 K,彈丸攻角為3°,通過計算流體力學(CFD)軟件計算結果如圖9~圖13所示。

殘差結果如圖9所示,空氣密度分布等值圖如圖10所示,彈丸I、彈丸II壓力分布如圖11所示,圖12為上下邊界所受的壓力圖。

圖9 殘差曲線圖Fig.9 Residual error plot

圖10 壓力分布等值圖Fig.10 Pressure distribution contour map

圖11 彈丸Ⅰ與彈丸Ⅱ壓力分布圖Fig.11 Pressure profile of Projectile Ⅰ and Projectile Ⅱ

圖12 邊界壓力分布圖Fig.12 Pressure distribution at boundary

圖13 仿真擬合數據Fig.13 Simulated data and its fitted curves

通過計算仿真,根據彈丸飛行參數,對指示器在距靶心5 m、10 m、15 m、20 m、30 m的位置上抽取每組3項彈丸的壓力幅值如表2所示。

通過壓力解算得出反演的彈目偏差,采用粒子群算法,擬合結果如圖13所示。得出:RMSE為0.302 236 05,SSE為1.370 199 46,R為0.999 458 7,R2為0.998 917 66,DC為0.998 758 82,Chi-Square系數為0.056 691 626,F-Statistic統計為1 997.941 786,將彈目偏差結果與實測比較,得出誤差小于5%,壓力幅值反演誤差小于4%.

表2 仿真反演數據及誤差Tab.2 Inversion data and its errors

注:仿真反演計算時忽略彈速、氣溫的微小差別。

7 結論

仿真反演是氣動聲學研究的重要手段,通過數值仿真和外場試驗的結合能更有效地檢驗理論研究,優化工程應用效果。

1)實測試驗數據與擬合的理論數據表明基于彈丸激波的彈目偏差測試精度較高,誤差為7%以內,符合實際需求。

2)對2發超聲速彈丸鄰近干擾情況進行了仿真模擬和計算表明,采用壓力解算仿真反演得出的彈目偏差與外場試驗接近,得出距離誤差小于5%,壓力幅值反演誤差小于4%.

研究結果表明,仿真反演與工程測試相符,方法可行,為超聲速運動目標聲學定位和跟蹤提供了一種新思路和方法。

References)

[1] DV Gaitonde. Progress in shock wave/boundary layer interactions[J]. Progress in Aerospace Sciences,2013(72):80-99

[2] Libal U, Spyra K. Wavelet based shock wave and muzzle blast classification for different supersonic projectiles[J]. Expert Systems with Applications, 2014, 41(11):5097-5104.

[3] Heimbs S, Ritzer J, Markmiller J. A numerical method for blast shock wave analysis of missile launch from aircraft[J]. International Journal of Aerospace Engineering, 2015, 2015:1-8.

[4] Aswin G, Chakraborty D. Numerical simulation of tansverse side jet interaction with supersonic free stream[J]. Aerospace Science and Technology, 2010, 14(5):295-301.

[5] 楊訓仁, 陳宇. 大氣聲學[M]. 北京:科學出版社, 2007:256-257. YANG Xun-ren, CHEN Yu. Atmospheric sound[M]. Beijing:Science Press, 2007:256-257.(in Chinese)

[6] 國蓉, 何鎮安, 王偉. 被動聲探測技術與彈著點定位方法綜述[J]. 電聲技術, 2010, 34(11):48-49. GUO Rong, HE Zhen-an, WANG Wei. Review of passive acoustic localization technique and impact point location methods[J]. Elementary Electroacoustics, 2010, 34(11):48-49.(in Chinese)

[7] 于雪媛, 左乾縣, 李永鋒. 一種測試彈丸落炸點的數學模型[J]. 火炮發射與控制學報, 2008(4):123-126. YU Xue-yuan, ZUO Qian-xian, LI Yong-feng. A mathematical model of gun ball fall explosion point measurement[J]. Journal of Gun Launch & Control,2008(4):123-126.(in Chinese)

[8] 王英, 曾光宇. 雙線陣CCD交匯測量立靶精度系統研究[J]. 光電工程, 2011, 38(10):33-38. WANG ying, ZENG Guang-yu. Intersection measuring system of erecting target with dual liner CCD[J]. Opto-Electronic Engineering, 2011, 38(10):33-38.(in Chinese)

[9] 張先葉, 高俊釵, 李亞強. 彈丸立靶坐標的圖像處理算法研究[J]. 科學技術與工程, 2008, 8(11):3018-3021. ZHANG Xian-ye, GAO Jun-cha, LI ya-qiang. Image processing algorithm of the pill coordinates[J]. Seience Technology and Engineering, 2008, 8(11):3018-3021.(in Chinese)

[10] 田會, 倪晉平, 魯倩, 等. 光幕與圖像融合式自動報靶系統[J]. 西安工業大學學報, 2012, 32(8):617-621. TIAN hui, NI jin-ping, LU qian, et al. Automatic target scoring system based on light screen and image porcessing[J]. Journal of Xi’an Technological University, 2012, 32(8):617-621.(in Chinese)

[11] 顧國華, 陳棟, 王彩, 等. 基于聲學靶傳感器的彈著點測試研究與實現[J]. 電子測量技術, 2007, 30(2): 154-156. GU Guo-hua, CHEN Dong, WANG Cai, et al. Research and implement of acoustics testing impact points based on target sensor[J]. Electronic Measurement Technology, 2007, 30(2):154-156.(in Chinese)

[12] 張飛猛, 馬春茂. 對空射擊聲學靶脫靶量測試系統的精度分析[J]. 兵工學報, 2000, 21(1):23-26. ZHANG Fei-meng, MA Chun-mao. Accuracy analysis for the measuring system of target deviation for projectiles shooting an acoustic target[J]. Acta Armamentarii, 2000, 21(1):23-26.(in Chinese)

[13] 張飛猛, 徐長根. 高炮射擊效果實時評價系統測量誤差修正[J]. 測試技術學報, 2007, 21(6):505-509. ZHANG Fei-meng, XU Chang-gen. Test error revision of real time evaluation system for anti-aircraft firing effect[J]. Journal of Test and Measurement Technology, 2007, 21(6):505-509.(in Chinese)

[14] Levin D, Gannot S, Habet S. Direction-of-arrival estimation using acoustic vector sensors in the presence of noise[C]∥2011 IEEE International Conference on Acoustics, Speech and Signal Processing. Prague, CZ:IEEE, 2011:105-108.

[15] 張飛猛. 火炮對懸吊靶射擊的脫靶向量測試系統[J]. 測試技術學報, 1999, 13(1):28-31. ZHANG Fei-meng. The missing target vector measurement system of the hang targets fired by the gun[J]. Journal of Test and Measurement Technology, 1999,13(1):28-31.(in Chinese)

[16] 吳壽榮. 球面激波對平面的馬赫反射的簡化處理[J]. 空氣動力學報, 1990(4):480-481. WU Shou-rong. Simple disposal of mach reflection of spherical shock wave on a plane[J]. Acta Aerodynamica Sinica, 1990(4):480-481.(in Chinese)

[17] Wagner B, Schmidt W. Theoretical investigations on shock wave-boundary layer interaction in cryogenic nitrogen[J]. Springer International Publishing, 2015, 365(3):341-350.

[18] 陶如意, 王浩, 趙潤祥. 超音速子母彈分離激波干擾特性研究[J]. 兵工學報, 2011, 32(10):1206-1211. TAO Ru-yi, WANG Hao, ZHAO Run-xiang, et al. Research on shock/shock wave disturbance characteristics in separation of supersonic cluster munition[J]. Acta Armamentarii, 2011, 32(10):1206-1211.(in Chinese)

[19] Kremeyer K . Energy deposition systems, equipment and method for modifying and controlling shock waves and supersonic flow: US, 8960596B2[P].2015-2-24.

[20] Lombard B, Matignon D, Gorrec Y L. A fractional Burgers equation arising in nonlinear acoustics: theory and numerics[J]. IFAC Proceedings Volumes, 2013, 46(23):406-411.

[21] 王東生. 炮口波、彈道波、爆炸波[J]. 兵器知識, 1985(3):23. WANG Dong-sheng. Muzzle wave, ballistic wave and blast wave[J]. Ordnance knowledge, 1985(3):23.(in Chinese)

[22] 王南炎. 彈丸空氣動力學[M]. 南京:中國人民解放軍總字150部隊, 1962:148-155. WANG Nan-yan. Projectile aerodynamics lecture [M]. Nanjing:Unit 150 of PLA, 1962:148-155.(in Chinese)

[23] 王保國. 空氣動力學基礎[M]. 北京:國防工業出版社, 2014:122-123. WANG Bao-guo. Aerodynamics basis[M]. Beijing:National Defense Industry Press, 2014:122-123.(in Chinese)

[24] Gnani F, Lo K H, Zare-Behtash H, et al. Shock wave diffraction in the presence of a supersonic co-flow jet[J]. Shock Waves, 2016, 26(3):1-10.

[25] 張強. 氣動聲學基礎[M]. 北京:國防工業出版社,2012:145-147. ZHANG Qiang. Aeroacoustics foundation[M]. Beijing:National Defense Industry Press, 2012:145-147.(in Chinese)

[26] 張金明. 氣動聲學數值計算的算法研究[D]. 南京:南京航空航天大學, 2009:6-8. ZHANG Jin-ming. Research of numerical method for computational aeroacoustics[D]. Nanjing:Nanjing University of Aeronautics and Astronautics, 2009:6-8.(in Chinese)

[27] Sofronov I L. Non-reflecting inflow and outflow in a wind tunnel for transonic time-accurate simulation[J]. Journal of Mathematical Analysis & Applications, 2010, 221(1):92-115.

Research on Miss Distance Test Based on Projectile Shock Wave and Its Inversion

GU Guo-hua, ZHANG Liang, WU Hai-bing, CUI Xun-xue, YU Bo

(Army Officer Academy, Hefei 230031, Anhui, China)

Since the shock wave of supersonic projectile is difficult to be captured due to its complexity and susceptible interference, the miss distance test is a difficult problem. The testing scene of miss distance is described, the formation mechanism of shock wave is analyzed, and the parameter relations before and after shock wave forming are identified. The test model of miss distance and the linear and nonlinear interference models based on simplified shock wave environment are established. The time duration test principle is used for field experiment, and the cumulative effect of pressure of shock wave field is used for simulating inversion in order to verify the accuracy and effectiveness of test methods and engineering test.

ordnance science and technology; ballistic shock wave;miss distance;acoustic field inversion; acoustic pressure

2016-06-03

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

顧國華(1981—),男,博士研究生。E-mail:Lmonkey2010@sohu.com

崔遜學(1969—),男,教授,碩士生導師。E-mail:xxcui@tsinghua.org.cn

TJ011+.2

A

1000-1093(2017)03-0576-09

10.3969/j.issn.1000-1093.2017.03.022

主站蜘蛛池模板: 毛片在线播放a| 欧美午夜在线视频| 中文字幕久久亚洲一区| 欧美无专区| 欧美丝袜高跟鞋一区二区 | 国产大片黄在线观看| 欧美精品在线看| 久久久亚洲国产美女国产盗摄| 久久综合色天堂av| 54pao国产成人免费视频| 天堂在线www网亚洲| 久久99精品久久久久久不卡| 精品国产毛片| 色哟哟精品无码网站在线播放视频| 欧美日韩国产成人高清视频| 国产色图在线观看| 亚洲精选高清无码| 国产一国产一有一级毛片视频| 91人妻日韩人妻无码专区精品| 国产精品青青| 成人午夜视频免费看欧美| 91青青视频| 亚洲色图欧美在线| 色综合中文| 无码高潮喷水专区久久| 亚洲精品日产精品乱码不卡| 成人在线视频一区| 中文字幕无码制服中字| 尤物成AV人片在线观看| 女人18毛片一级毛片在线| 欧美黄网在线| 国产欧美日韩综合在线第一| 亚洲成人网在线播放| 亚洲成年人网| 婷婷午夜天| 欧美久久网| 国产成人精品日本亚洲| 九色视频线上播放| 日韩精品一区二区三区swag| 国产亚洲精品精品精品| 成人福利在线免费观看| 91精选国产大片| 日韩天堂在线观看| 国产第一福利影院| 在线免费a视频| 谁有在线观看日韩亚洲最新视频| 日a本亚洲中文在线观看| 免费看美女毛片| 亚洲热线99精品视频| 女同国产精品一区二区| 亚洲男人的天堂在线| 久久精品人人做人人爽| 国产欧美日韩精品第二区| 国产不卡网| 欧美在线视频不卡| 精品国产自在现线看久久| 国产青榴视频| 国产经典免费播放视频| 精品视频一区二区观看| 日本成人一区| 亚洲AV无码乱码在线观看代蜜桃| 性69交片免费看| 日韩色图区| 内射人妻无套中出无码| 久久精品视频亚洲| 色视频国产| 亚洲专区一区二区在线观看| 四虎精品免费久久| 国产黄色爱视频| 欧美视频在线播放观看免费福利资源 | 国产一区二区三区在线精品专区 | 国产黑丝一区| 日韩天堂网| 四虎综合网| 中文精品久久久久国产网址| 国产永久无码观看在线| 中文字幕永久视频| 亚洲天堂网在线视频| av天堂最新版在线| 久久综合五月| 在线观看精品自拍视频| 亚洲人成网18禁|