亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        回轉體傾斜入水空泡試驗及六自由度數(shù)值計算研究

        2017-12-12 19:49:56孫鐵志張桂勇
        宇航總體技術 2017年4期
        關鍵詞:空泡加速度角度

        侯 昭,孫鐵志,張桂勇,2,3,宗 智,2,3

        (1.大連理工大學船舶工程學院,遼寧省深海浮動結構工程實驗室,大連116024; 2.大連理工大學工業(yè)裝備與結構分析國家重點實驗室,大連116024; 3.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,上海200240)

        回轉體傾斜入水空泡試驗及六自由度數(shù)值計算研究

        侯 昭1,孫鐵志1,張桂勇1,2,3,宗 智1,2,3

        (1.大連理工大學船舶工程學院,遼寧省深海浮動結構工程實驗室,大連116024; 2.大連理工大學工業(yè)裝備與結構分析國家重點實驗室,大連116024; 3.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,上海200240)

        為研究回轉體傾斜入水空泡及參數(shù)變化,基于低速入水試驗裝置和CFD軟件對該問題進行了試驗和數(shù)值仿真研究。通過試驗對回轉體傾斜入水的空泡演變進行研究,得到了不同時刻空泡形態(tài)圖。數(shù)值計算選用基于N-S方程的雷諾平均(RANS)方法和基于k-ω的SST二方程湍流模型,建立了六自由度數(shù)值仿真方法。結果表明:數(shù)值計算得到的空泡形態(tài)與試驗結果一致性較好,依次經歷了撞擊水面、空泡形成、頸縮、空泡斷裂、空泡閉合、表面紊亂和空泡潰滅的過程。進一步分析發(fā)現(xiàn)速度、加速度和壓力均在入水瞬間和空泡斷裂時刻發(fā)生波動;偏轉角度在回轉體尾部刺穿空泡后增幅明顯;隨著入水角度增加,入水瞬間速度衰減加快、壓力峰值增加,峰值出現(xiàn)越早,空泡閉合越難。

        傾斜入水;數(shù)值計算;六自由度;彈道

        0 引言

        回轉體傾斜入水不可避免地要經歷撞擊水面、空泡形成、頸縮、空泡斷裂、空泡閉合、表面紊亂和空泡潰滅的過程。入水瞬間產生的沖擊載荷可能對回轉體結構造成破壞[1];水下航行過程中產生的空泡影響其在水中的航行姿態(tài)及彈道穩(wěn)定性[2-3],甚至可能發(fā)生沉底、彈跳等問題。因此,研究入水空泡問題能夠對提高回轉體航行性能和結構穩(wěn)定提供重要參考。

        對入水問題的試驗研究始于19世紀末,攝影技術的高速發(fā)展使得空泡從形成到潰滅過程得以被記錄研究。Worthington等[4]利用高速攝影技術,首次記錄了剛性球垂直入水空泡演變過程。隨后,Watanabe[5]在前人的基礎上進行了不同頭型結構低速入水試驗,研究了瞬時沖擊特性、最大阻力與入水速度的關系。針對不同頭部形狀回轉體在不同角度下入水的工況,Abelson[6]在1970年對其所形成的空泡壁面內壓力進行研究分析。May[7-8]對彈體垂直以及傾斜入水時空泡演變過程進行了研究,同時分析了彈體及流體的受力情況。何春濤等[9]進行了不同速度、角度的單圓柱體以及多圓柱體串行、并行入水空泡試驗。此外,國內諸多學者[10-13]還針對不同物體入水空泡演變過程開展了大量研究。隨著計算機的發(fā)展和入水空泡理論的成熟,數(shù)值計算逐漸成為入水空泡問題的重要研究方法。石漢成等[14]通過MSC Dytran研究了頭部形狀對水雷入水載荷及水下彈道的影響。He等[2]對垂直入水空泡內部壓強分布進行了數(shù)值仿真。張國軍等[15]基于VOF模型對導彈低速入水進行仿真模擬。程文鑫等[16]對魚雷在小入水角條件下的入水空泡形成過程進行了計算。朱珠等[17]仿真分析了入水攻角和速度對入水彈道的影響規(guī)律。

        本文對低速傾斜入水空泡問題進行了試驗,同時利用CFD軟件對試驗過程進行了數(shù)值仿真。入水過程涉及多相流的非線性問題,同時屬于流固耦合范疇,數(shù)值計算過程采用RANS方法和基于k-ω的SST二方程湍流模型,利用VOF(Volume of Fluid)模型、歐拉多相模型來描述流體體積及模擬空泡形成至潰滅的過程,采用VOF波捕捉自由液面并利用六自由度體定義回轉體,通過重疊網格技術模擬回轉體入水過程,從而建立了回轉體傾斜入水空泡問題的數(shù)值仿真方法,并與試驗結果進行對比驗證。在此基礎上進一步分析了傾斜入水過程中彈道、流體動力以及入水空泡演變特性。

        1 試驗系統(tǒng)和數(shù)值計算模型

        1.1 試驗系統(tǒng)

        試驗系統(tǒng)設計如圖1(a)所示,系統(tǒng)主要由底部支架、試驗水箱、滑軌機構、燈光系統(tǒng)、高速攝像機和計算機組成。試驗水箱由鋼化玻璃板制成,尺寸為1.5m×0.8m×1.2m,箱內水深為1m,試驗選用20℃的自來水;支架置于水箱底部,底部裝有滾輪,達到支撐固定實驗水箱及方便移動的目的;滑軌機構是實現(xiàn)回轉體以一定速度及角度入射的基礎,滑軌一端固定在水槽上,另一端固定在可調節(jié)支架上,通過調整回轉體在滑軌上的位置來改變入水速度,調整伸縮桿高度可以改變入水角度;入水空泡演變過程采用Photron高速工業(yè)攝像機記錄,將其放置在試驗水箱正前方,聯(lián)合使用計算機和高速攝像機專用軟件PCC進行數(shù)據(jù)采集和處理;在高速攝像機兩邊分別設置兩個Starison CE-1500Ws透射式聚光燈,構成燈光系統(tǒng),保證圖片采集質量。用于試驗的回轉體模型為鋁制實心圓柱體,密度為2740kg/m3,長度為0.197m,質量為1.06kg,直徑為0.05m。

        圖1 試驗裝置及回轉體模型Fig.1 Schematic of the water-entry experiment and revolution model

        1.2 數(shù)值計算模型

        1.2.1 基本控制方程

        本文計算所涉及的控制方程主要包括質量守恒方程和動量守恒方程。

        質量守恒方程:

        其中,ρ為流體密度,本文涉及氣、液兩相,故密度根據(jù)各相所占體積分數(shù)來確定;t為時間;u i為沿x、y、z方向上的速度分量;x i為在x、y、z方向上的位置。

        動量守恒方程:

        1.2.2 湍流模型

        RANS方程中由于雷諾應力的出現(xiàn),增加了未知量數(shù)量,因此需要采用特定的湍流模型來封閉方程。本文選用渦黏模型中基于k-ω的SST二方程湍流模型:

        其中,k為湍動能,ω為湍流耗散率,S k和Sω是用戶指定的源項,k0和ω0是用來抵消湍流衰減的環(huán)境湍流值,μ為動力黏性系數(shù);γeff是有效間隔,定義γ′=min[max(γeff,0.1),1]。

        二方程湍流模型中Gω根據(jù)標準k-ω模型估算得到:

        其中,γ為模型的混合系數(shù),S為平均應變張量的模。

        Dω為交叉導數(shù)項,定義為:

        其中,T是湍流時間尺度,在Durbin約束下計算得到;系數(shù)α*=1;系數(shù)a1=0.31;函數(shù)F2=為到近壁面的距離。

        模型系數(shù)?由混合函數(shù)F1計算得到:

        其中,系數(shù)?1包括:

        混合函數(shù)F1定義為:

        1.2.3 計算模型及設置

        本文采用的物理模型和計算水域均與試驗相同,以便進行對比驗證。網格劃分是整個建模仿真的基礎,網格質量優(yōu)劣直接關系計算結果的準確性。本文采用網格生成質量比較好的切割體網格生成器,加以表面重構提高網格質量。在壁面邊界加設棱柱層網格生成器,設置棱柱層厚度為0.5mm,棱柱層數(shù)為5層。為保證仿真模擬的質量,在回轉體運動彈道、水面以及回轉體1.5倍半徑范圍內進行網格加密。網格劃分時保持重疊區(qū)域垂直,網格生成后再對重疊區(qū)域旋轉平移到指定工況。本文計算的網格總數(shù)為3609899,網格劃分如圖2所示。

        圖2 網格劃分Fig.2 Mesh generation

        針對低速回轉體入水問題,采用重疊網格技術[18]實現(xiàn)對物體運動的模擬。重疊網格就是嵌套于流體域網格中、隨物體運動的網格,具有網格易生成、質量好的特點。重疊網格通過搜索指定區(qū)域,進行網格劃分;在物體運動過程中的每個時間步都需要將重疊區(qū)域從流體域內逐步剔除重新定位,以此來模擬入水過程;重疊網格與流體網格之間通過數(shù)據(jù)的插值進行信息傳遞。重疊網格局部如圖3所示。

        圖3 重疊網格Fig.3 Overlapping grid

        由于入水空泡問題涉及多相流,本文采用基于歐拉-歐拉型方法的VOF模型進行計算。將多相流中的流體視為均勻物質,通過各相流體所占體積分數(shù)來描述不同流體。體積分數(shù)總和為1,體積分數(shù)為1時代表液相,為0時代表氣相,處于0到1之間表示氣液的混合物。依此來表征各相流體體積,模擬低速回轉體入水過程中的空泡演變。

        本文數(shù)值計算過程采用VOF波定義的平波來捕捉自由液面,模擬氣液交界面。數(shù)值計算邊界條件為:入口設置為速度入口;出口為壓力出口;根據(jù)厲尚書[19]研究,為避免入水過程壓力壁面條件發(fā)生反射,將壁面設置為速度入口;重疊區(qū)域設定為重疊網格邊界,回轉體表面設定為壁面條件。

        2 結果與分析

        試驗中回轉體以4.85m/s初速度、60°入射角(與水平面的夾角)的初始條件入水,得到的空泡演變過程如圖4(a)所示。從圖4中可以看出,入水空泡的演變分為以下幾個階段:撞擊水面、空泡形成、頸縮、空泡斷裂空泡閉合、空泡表面紊亂和空泡潰滅。隨著回轉體入水,其周圍開始形成空泡,如25ms時的狀態(tài);入水50ms,回轉體完全被空泡包裹,空泡隨回轉體運動,逐漸擴展;85ms時,在尾部形成回射流;110ms時,伴隨入水深度的增加,流場的壓力及空泡的表面張力使空泡直徑減小,導致出現(xiàn)頸縮現(xiàn)象;入水150ms時,空泡完全閉合,尾部空泡脫落并在液體表面處形成漏斗形噴濺冠,閉合空泡受慣性影響跟隨回轉體繼續(xù)運動,受流場沖擊作用,空泡表面逐漸紊亂;180ms時,空泡潰滅為氣泡沿回轉體上浮;在整個入水過程中,液體表面都存在噴濺現(xiàn)象。

        數(shù)值計算過程首先采用與上述試驗相同的工況,表征空泡外形的含氣率的量值設定為0.5,從而得到圖4(b)中空泡演變過程,并與不同時刻回轉體入水空泡形態(tài)的試驗結果進行對比。從圖4中可以看出,數(shù)值計算得到的空泡形態(tài)演變與試驗結果一致性較好,包括入水空泡的形成、頸縮現(xiàn)象、空泡閉合、空泡脫落及潰滅、水面形成噴濺,從而驗證了所建立的數(shù)值計算方法。

        圖4 試驗與數(shù)值計算入水空泡演變對比Fig.4 Comparison of cavity evolution between the experimental results and numerical predictions

        圖5給出了入水過程典型時刻回轉體對稱面壓力分布云圖。從圖5中可以看出:在入水瞬間,自由液面對回轉體頭部產生較大壓力??张蓍]合前,近水面處壓力較大,因為此時空泡頸縮加劇并產生回射流,導致大量氣體向近水面移動,故壓力增加;另外,受空泡形態(tài)影響回轉體迎流面的低壓區(qū)明顯比背流面大。在0.15s時空泡閉合,近水面高壓區(qū)隨著氣流的流出而消失;回轉體迎流面低壓區(qū)與背流面尺寸差別縮小??张蓍]合后,近水面與回轉體之間形成長條狀低壓區(qū),這由少量氣泡脫落引起;同時,外部水壓大于空泡內部壓力并擠壓空泡,導致迎流面低壓區(qū)與背流面大小基本相同,空泡內部壓力相較空泡閉合前增加。

        圖5 入水過程典型時刻壓力分布云圖Fig.5 Pressure diagram at typical time of water-entry process

        為進一步探究傾斜入水空泡演變特征和機理,對回轉體傾斜入水過程中速度、加速度、彈道以及受力情況進行深入分析。圖6給出了入水過程中速度、加速度變化曲線。從圖6中可以看出:從接觸水面至入水0.01s過程,速度下降比較大,加速度迅速降低至0以下,這是入水瞬間產生巨大沖擊而后沖擊減緩引起的,這個過程中回轉體受到阻力也在瞬間增加而后減小,所以在初始時刻速度、加速度變化劇烈;回轉體入水后的0.01s~0.11s期間,受重力、慣性力和阻力聯(lián)合作用,速度基本上呈直線下降、加速度緩慢減小;0.11s~0.13s期間,加速度減小趨勢略增,速度變化與之前相差不大;0.13s~0.15s期間,加速度明顯減小,速度基本不變;0.15s~0.26s時間內,速度緩慢減小,加速度緩慢變化。

        回轉體入水過程中速度、加速度變化與空泡形態(tài)密切相關。由圖7可知,0.11s~0.13s時間內,空泡頸縮加劇,彈體尾部刺穿空泡,此時空泡開始閉合,回轉體受阻力略增,加速度有所減小但幅度不大,因此速度減小趨勢基本不變; 0.13s~0.15s期間,頭部空泡大量向尾部移動,回轉體受阻力大幅增加,造成加速度明顯減小且趨于0,速度基本保持不變,在0.15s時發(fā)生空泡斷裂,阻力瞬間減小,加速度有增加趨勢;0.15s ~0.26s期間,空泡跟隨回轉體運動,有少量氣泡脫落,阻力變化不大,因而加速度緩慢變化,速度緩慢減小。

        圖6 速度、加速度曲線Fig.6 Time histories of velocity and acceleration

        圖7 0.11s~0.26s時間內空泡形態(tài)演變Fig.7 The cavity evolution from 0.11s to 0.26s

        為分析入水后回轉體運動軌跡的變化,對其運動的偏轉角度 (與水平面的夾角)進行了討論,偏轉角度變化如圖8所示。從回轉體入水直至0.11s的過程中偏轉角緩慢增大,0.11s之后偏轉角度增幅明顯,逐漸趨近于90°。這是由于入水初期速度比較大,回轉體運動受慣性影響大且空泡正處于發(fā)展期,從而保持原有運動軌跡,偏轉角變化較小;由圖7可以看出,從0.11s開始,回轉體刺穿空泡,頸縮加劇,頭部空泡逐漸向尾部移動直至斷裂,造成阻力增加,速度逐漸減小,慣性的影響逐漸減弱,回轉體受重力影響大,因而偏轉角度增幅變大,逐漸趨于與水平面垂直狀態(tài)。

        圖8 偏轉角隨時間變化Fig.8 Time history of deflection angle

        圖9給出了回轉體入水過程中頭部壓力變化。可以看出回轉體撞擊水面時壓力急劇增加后又迅速減小,可見入水瞬間產生較大的沖擊載荷;隨著入水空泡的形成壓力有所減小;在0.12s~ 0.15s期間發(fā)生較大起伏,主要原因是此時頭部空泡向尾部運動并發(fā)生脫落,從而影響回轉體頭部與水作用的受力情況;0.15s后空泡隨物體運動,中心點壓力受水深的增加而小幅度增大。

        圖9 頭部中心點壓力曲線Fig.9 The pressure at the center of the head

        為研究不同入水角度對回轉體彈道特性的影響?;诮⒌臄?shù)值方法,同時開展了45°和75°傾斜入水工況的計算,其中不同角度入水的初速度相同。圖10對比了不同入水角度下回轉體頭部中心點壓力變化。從圖10中可以看出:不同入水角度下,回轉體頭部中心點壓力變化曲線趨勢基本一致;隨著入水角度的增加,入水瞬間回轉體頭部所受沖擊壓力增大;入水角度越大,壓力峰值和波動出現(xiàn)時間越早,這是因為大角度入水時頭部與自由液面接觸面大,故壓力增加迅速;在0.13s附近,由于大角度入水過程空泡發(fā)生斷裂,頸縮時間較早,導致出現(xiàn)波動壓力峰值時刻靠前。

        圖10 不同入水角壓力變化曲線Fig.10 Time histories of pressure at different angles

        圖11對比了不同入水角度下速度變化。從圖11中可以得出:回轉體在不同入水角度下,速度變化曲線趨勢基本一致;入水角度越大,入水瞬間回轉體所受沖擊壓力越大,導致入水速度衰減加快;隨著入水角度的增加,在速度曲線中變化緩慢的區(qū)域 (圖中a區(qū)域)持續(xù)時間更久,即空泡從頸縮加劇到閉合的過程持續(xù)時間長。

        圖11 不同入水角速度變化曲線Fig.11 Time histories of velocity at different angles

        3 結論

        本文對回轉體傾斜入水空泡問題進行了試驗和數(shù)值計算研究,對入水空泡演變、彈道變化以及受力情況進行了深入探討,得到的主要結論如下:

        1)回轉體傾斜入水依次經過了:撞擊水面、空泡開始形成、頸縮、空泡完全閉合、空泡表面紊亂、空泡潰滅等過程。

        2)入水瞬間壓力、速度、加速度受沖擊載荷影響產生較大波動。隨后,在空泡敞開階段,受壓差影響外部空氣進入空泡,空泡外壁受水壓影響逐漸閉合,閉合空泡受高速氣流進入的影響導致內部壓力降低。由于閉合空泡內氣壓的突然下降,導致回轉體的速度、加速度和壓力出現(xiàn)波動。回轉體偏轉角則從尾部刺穿空泡后開始增幅變大。

        3)不同入水角度對回轉體傾斜入水的影響主要體現(xiàn)在:隨著入水角度的增加,回轉體入水瞬間速度衰減加劇,壓力峰值增大,空泡從頸縮加劇到閉合的過程持續(xù)時間增長,壓力峰值和波動出現(xiàn)的時間提前。

        [1] 鄭金偉,宗智.三維剛體橢圓頭結構高速傾斜入水沖擊模擬[J].船海工程,2012,41(3):7-9.

        [2] He C T,Wang C,Wei Y J,et al.Numerical simulation of pressure distribution in vertical water-entry cavity[J].Journal of Ship Mechanics,2011,15 (9):960-968.

        [3] 王聰,何春濤,權曉波,魏英杰.空氣壓強對垂直入水空泡影響的數(shù)值研究 [J].哈爾濱工業(yè)大學學報, 2012,44(5):14-19.

        [4] Worthington A M,Irwin K G.A study of splashes[J]. Science,1909,29(742):464-465.

        [5] Watanabe S.Resistance of impact on water surface[J]. Scientific Papers of the Institute of Physical and Chemical Research of Japan,1934,23(484): 202-208.

        [6] Abelson H I.Pressure measurements in the waterentry cavity[J].Journal of Fluid Mechanics,1970, 44(1):129-144.

        [7] May A.Review of water-entry theory and data[J]. Journal of Hydronautics,1970,4(4):140-142.

        [8] May A.Water entry and the cavity-running behavior of missiles[R].Nasa Sti/recon Technical Report N,1975.

        [9] 何春濤,王聰,何乾坤,等.圓柱體低速入水空泡試驗研究[J].物理學報,2012,61(13):281-228.

        [10] 顧建農,張志宏,王沖,等.旋轉彈頭水平入水空泡及彈道的試驗研究[J].兵工學報,2012,33(5): 540-544.

        [11] 楊衡,張阿漫,龔小超,等.不同頭型彈體低速入水空泡試驗研究[J].哈爾濱工程大學學報,2014, 35(9):1060-1066.

        [12] 路中磊,魏英杰,王聰,等.基于高速攝像試驗的開放腔體圓柱殼入水空泡流動研究 [J].物理學報, 2016,65(1):301-315.

        [13] 蔣運華,徐勝利,周杰.圓盤空化器航行體入水空泡試驗研究[J].工程力學,2017,34(3):241-246.

        [14] 石漢成,蔣培,程錦房.頭部形狀對水雷入水載荷及水下彈道影響的數(shù)值仿真分析 [J].艦船科學技術, 2010,32(10):104-107.

        [15] 張國軍,閆云聚.基于VOF模型的導彈低速入水數(shù)值模擬方法 [J].空軍工程大學學報 (自然科學版),2013,14(6):23-26.

        [16] 程文鑫,蔡衛(wèi)軍,楊春武.魚雷小角度入水過程仿真[J].魚雷技術,2014,22(3):161-164.

        [17] 朱珠,袁緒龍,劉維.柱體大攻角入水彈道建模與仿真[J].火力與指揮控制,2015,40(2):13-18.

        [18] 徐偉光,趙發(fā)明,何術龍.基于重疊網格的船舶運動計算方法研究[J].船舶力學,2016,20(7): 824-832.

        [19] 厲尚書.空投魚雷入水沖擊數(shù)值模擬的研究[D].太原:中北大學,2015.

        Experimental Investigation and 6-DOF Simulation of Oblique Water-entry Cavity of Revolution Body

        HOU Zhao1,SUN Tie-zhi1,ZHANG Gui-yong1,2,3,ZONG Zhi1,2,3
        (1.Liaoning Engineering Laboratory for Deep-Sea Floating Structures,School of Naval Architecture, Dalian University of Technology,Dalian 116024,China;2.State Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology,Dalian 116024,China;3.Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration,Shanghai 200240,China)

        In order to investigate the problem of oblique water entry cavity and its related parameters,both experiment and numerical simulation of oblique water-entry process have been carried out.In the experiment,the cavity evolution of oblique water entry was studied and the cavity morphology at different times was obtained.The numerical model was constructed using the RANS method and SST two-equation turbulence model based on k-ωmodel,which can effectively simulate the 6-DOF movement of the revolution body.The results show that the cavity evolution obtained by numerical simulation agrees well with the experimental results,which mainly includes the stages of impacting water surface,forming cavity,appearing pinch off phenomenon,cavity cracking,cavity closure,coming out surface disordered and cavity collapse.The change of waterentry parameters has been further analyzed based on the numerical simulation.It was found that the velocity,acceleration and pressure all change abruptly when the revolution body impacts water surface and the cavity cracks.The body's deflection angle shows an increasing trend,which becomes more significant after the tail of the revolution pierces the cavity.With the increase of incident angle,the velocity of the revolution body decreases rapidly and the peak pressure increases greatly when the revolution body impacts water surface.Besides,the peak pressure appears earlier and the cavity is more difficult to be closed.

        Oblique water entry;Numerical simulation;Six degrees of freedom;Ballistic trajectory

        O352

        A

        2096-4080(2017)04-0038-08

        2017-09-12;

        2017-10-28

        國家自然科學基金(51579042,51639003,51709042);青年千人項目(D1007001);中央高?;究蒲袠I(yè)務費專項資金(DUT16ZD218,DUT17ZD311,DUT16RC(3)085)

        侯昭(1993-),女,碩士研究生,主要研究方向為回轉體入水流體動力特性研究。E-mail:619256848@qq.com

        張桂勇 (1978-),男,博士,教授,主要研究方向為無網格方法、船舶與海洋工程流固耦合分析及極地船舶等高技術船舶開發(fā)技術。E-mail:gyzhang@dlut.edu.cn

        猜你喜歡
        空泡加速度角度
        “鱉”不住了!從26元/斤飆至38元/斤,2022年甲魚能否再跑出“加速度”?
        當代水產(2022年6期)2022-06-29 01:12:20
        神奇的角度
        水下航行體雙空泡相互作用數(shù)值模擬研究
        天際加速度
        汽車觀察(2018年12期)2018-12-26 01:05:42
        一個涉及角度和的幾何不等式鏈的改進
        創(chuàng)新,動能轉換的“加速度”
        金橋(2018年4期)2018-09-26 02:24:46
        死亡加速度
        勞動保護(2018年8期)2018-09-12 01:16:14
        角度不同
        37°女人(2017年8期)2017-08-12 11:20:48
        人啊
        滇池(2017年7期)2017-07-18 19:32:42
        基于LPV的超空泡航行體H∞抗飽和控制
        久久久一本精品久久久一本| 亚洲 暴爽 av人人爽日日碰| 国产91成人精品亚洲精品| 亚洲成AV人片无码不卡| 久久精品亚洲一区二区三区画质| 精品国产一区二区三区2021| 中日韩精品视频在线观看| 91性视频| 午夜视频手机在线免费观看| 日本污ww视频网站| 中文字幕人妻av一区二区| 三级国产女主播在线观看| 一区二区三区国产天堂| 体验区试看120秒啪啪免费| 好吊色欧美一区二区三区四区| 国产二区中文字幕在线观看| 日本一二三四高清在线| 色综合色狠狠天天综合色| 亚洲av色欲色欲www| 欧洲人妻丰满av无码久久不卡| 野外三级国产在线观看| 在线亚洲妇色中文色综合| 无码无套少妇毛多18p| 乱子伦视频在线看| jiZZ国产在线女人水多| 青青青免费在线视频亚洲视频| 久久无码av一区二区三区| 日韩成人免费一级毛片| 丝袜美腿诱惑一二三区| 亚洲国产av无码精品无广告| 午夜亚洲av永久无码精品| 久久精品国产乱子伦多人| 亚洲第一女人的天堂av| 欧美又大又色又爽aaaa片| 国产欧美日韩午夜在线观看| av在线不卡一区二区三区| 日韩精品无码熟人妻视频| av无码久久久久久不卡网站| 男女上床视频在线观看| 91久久精品色伊人6882| 被群cao的合不拢腿h纯肉视频|