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

        ?

        基于CFD方法的對拍撲翼推進特性分析

        2014-09-12 11:22:24張偉偉蔣躍文葉正寅
        空氣動力學學報 2014年4期
        關鍵詞:振動效率研究

        張偉偉,徐 楊,蔣躍文,葉正寅

        (西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室,陜西西安 710072)

        基于CFD方法的對拍撲翼推進特性分析

        張偉偉,徐 楊,蔣躍文,葉正寅

        (西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室,陜西西安 710072)

        雙翼對拍推進相對于單撲翼推進具有效率高,機械設計容易,升力和慣性力平穩(wěn)的特點。論文采用雷諾平均 NS方程和非結構混合運動網(wǎng)格方法求解非定常流場,以 NACA0012翼型為例,研究了低速狀態(tài)下雙翼對拍時沉浮運動的振幅、頻率、雙翼的間距以及沉浮俯仰兩自由度相角差等參數(shù)對推進效率和推進系數(shù)的影響。研究發(fā)現(xiàn),相對單撲翼推進,雙翼對拍能夠得到更大的推力和更高的推進效率。附加俯仰自由度后,選擇合適的沉浮和俯仰相角差,兩自由度組合撲動系統(tǒng)的推進性能相比單純沉浮撲動系統(tǒng)又將進一步得到提高。

        撲翼推進;雙翼對拍;數(shù)值模擬;非定常流場

        0 引 言

        撲翼推進是一種有別于傳統(tǒng)旋轉機械(螺旋槳、渦輪發(fā)動機)和噴氣推進的又一種推進方法。隨著近年來微型飛行器和水下機器人研究的興起,這種推進方式顯示出來越來越多的優(yōu)勢。這主要是由于在低雷諾數(shù)下,傳統(tǒng)螺旋槳推進的效率較低,而這種撲動式非定常運動不僅能夠實現(xiàn)固定翼無法達到的高升力,而且其推進效率相比傳統(tǒng)方式也具有一定的優(yōu)越性。天空飛行的鳥類和水下游動的魚類都是采用這類非定常運動方式實現(xiàn)推進的。也正是這些生物卓越的飛行能力和游動能力啟發(fā)了人們對撲動推進的研究興趣。

        很早以前就有科學家對這個領域的研究進行了嘗試。Knoller[1]和Betz[2]分別在1909年和1912年發(fā)現(xiàn)了撲動的機翼可以產(chǎn)生當量迎角,而這個迎角產(chǎn)生的法向力在來流方向進行投影,便可產(chǎn)生推力。因此這種撲動機翼產(chǎn)生有效推力的效應就被稱為Knoller-Betz效應。在1922年,科學家 Katzmayr[3]通過實驗驗證了 Knoller-Betz效應。接下來的研究者通過豐富的手段對撲動機翼(包含單獨沉浮、俯仰、沉?。┭鲴詈希┑耐七M機理作進一步的研究。1935年,Von Karman和 Burgers[4]通過對尾跡旋渦的位置和方向的研究,對阻力和推力的產(chǎn)生提出了理論的解釋。幾乎在同一時間,Garricks[5]采用Theodorsen的無粘、不可壓縮振蕩平板理論解釋了撲動機翼的推力產(chǎn)生機理。Jones[6]研制出一個具有上下兩對均可產(chǎn)生拍打和俯仰運動的翼展為1270mm 的撲翼機構,并進行了風洞試驗。中國科學技術大學的童秉綱、陸夕云等人[7-8]采用活體觀測、實驗測量、數(shù)值模擬和理論分析多種手段,以及力學和生物學相結合的綜合研究方法,系統(tǒng)地開展了關于飛行和游動的生物運動力學以及相關的仿生技術研究。

        如今,CFD方法已經(jīng)廣泛用來模擬撲翼推進問題。國內外不少團隊在這方面作了較多的工作,進行各種魚類推進模式和鳥類飛行模式的機理研究,并且建立了多種仿生模型來模擬撲翼推進。楊亮,蘇玉民[9]等應用CFD技術開展了仿金槍魚擺動尾鰭的水動力性能與推進機理的研究。楊樹池[10]等人以 NACA0012翼型作為模型,驗證了其在低速情況下做沉浮及俯仰運動歐拉方程的適用性。大連理工大學的王濤[11]等利用非定常雙時間NS方程的預處理方法對雙翼地效推進器進行了推力和推進性能的分析,數(shù)值結果顯示了沉浮振動的推力性能。推力和效率與速度的變化關系表現(xiàn)出一致性,固定頻率,改變速度,發(fā)現(xiàn)推力和效率存 在極大值點 。余 春錦 ,昂海 松等[12]對柔性膜微型飛行器的氣動特性進行了數(shù)值模擬,考察了柔性對撲翼的推力和升力的影響。西北工業(yè)大學楊文青等人[13]研究了展向折疊撲動模式對提升微型飛行器撲動性能的影響。錢靖[14]等數(shù)值模擬后緣帶彈性薄板的翼型在沉浮運動時翼型繞流和薄板彈性振動之間的流固耦合問題及撲動運動的推力效應。

        近年來一種新型撲翼推進引起了人們新的關注,即雙翼對拍推進方式。這種推進主要利用了類似流體力學中的地面效應原理。利用這種鏡像效應,實現(xiàn)了一加一大于二的效果。另外這種機構的另一些優(yōu)勢在于通過雙翼的對稱運動特性可以極大程度抵消單個撲翼一個周期內結構慣性力和升力的脈動,不僅能夠有效改善微型飛行器的飛行品質,也為動力機構和系統(tǒng)結構的設計帶來很多方便。然而目前針對雙翼對拍推進的研究相 對不 多,僅 見Liu和王濤[14-16]的研 究,主要研究了振動頻率和來流速度對推進特性的影響。但在上述研究中,僅考慮了沉浮方向的一個自由度,也未對雙翼平衡間距對推進特性的影響進行分析。

        本論文在上述工作的基礎上,系統(tǒng)對比了單翼和雙翼撲動的推進特性,進一步研究了雙翼平衡間距對推進效率和推力系數(shù)的影響,并考慮沉浮和俯仰兩自由度系統(tǒng),研究了俯仰沉浮自由度相角差對推進特性的影響。

        1 數(shù)值仿真方法

        采用課題組自行開發(fā)的NWind2D計算程序,以二維可壓縮雷諾平均 N-S方程為控制方程??臻g離散采用非結構混合網(wǎng)格技術。繞翼型壁面附近采用結構化網(wǎng)格,用于描述附面層流動,其它部分采用非結構化網(wǎng)格。由于流動涉及邊界的運動,故而采用運動網(wǎng)格技術來實現(xiàn)網(wǎng)格的變形,避免了網(wǎng)格的重新生成和流場參數(shù)的插值。壁面附近的結構化網(wǎng)格運動與壁面運動同步,而非結構網(wǎng)格采用彈簧方法實現(xiàn)網(wǎng)格的變形。這樣的處理方式可以保證即使壁面大幅運動也不會出現(xiàn)附面層附近的網(wǎng)格破壞。流場求解的空間離散格式采用Jameson中心有限體積法,時間推進格式采用“雙時間”方法,湍流模型采用S-A一方程模型。關于該程序的具體實現(xiàn)過程和可靠性可參考文獻[17]。

        2 參數(shù)定義

        研究算例采用NACA0012翼型,研究中參數(shù)都進行無量綱化。定義雙翼布局模型的翼間距L=l/c,沉浮運動振幅斯特羅哈爾數(shù)為后緣點的最大運動路程,非定常運動的減縮頻率為翼型弦長。 定義一個運動周期內平均推力一個運動周期內平均輸入功率P應的推力系數(shù)輸入功率系數(shù)

        上下兩個翼型分別繞平衡位置作簡諧沉浮振動和繞軸點c/4處的簡諧俯仰運動,初始來流迎角始終為0,兩個翼型運動方向完全相反,故每個時刻兩翼與x軸都呈對稱分布(圖1)。運動定義如下:h(t)= hmsin(ωt-φ),α(t)=αmsin(ωt),針對上一片翼型,沉浮運動定義上翼向上運動為正方向,俯仰運動定義上翼抬頭為正方向。

        圖1 翼型非定常運動示意圖Fig.1 Sketch map for unsteady motion of airfoils

        遠場網(wǎng)格直徑20倍翼型弦長,附面層15層,附面層網(wǎng)格第一層高度為4×10-4。單翼布局網(wǎng)格的總網(wǎng)格個數(shù)是8398,附面層部分結構網(wǎng)格數(shù)目是3060;雙翼布局網(wǎng)格的總網(wǎng)格個數(shù)大約是15000,附面層部分結構網(wǎng)格數(shù)目是6120(圖2)。

        3 算例與分析

        3.1 程序驗證

        為了驗證程序的正確性,本文首先計算了兩個算例,計算模型均采用NACA0012翼型,計算結果和實驗值或參考文獻值的比較見圖3和圖4。

        (1)計算狀態(tài)為馬赫數(shù)M=0.6,雷諾數(shù)Re=9× 106,平均迎角α0=2.89°,最大振幅迎角αm=2.41°,減縮頻率k=0.0808,力矩的積分點為27.3%弦長處。

        圖2 單翼及雙翼模型網(wǎng)格Fig.2 Grids of single airfoil model and dual airfoils model

        圖3 升力系數(shù)和力矩系數(shù)隨迎角的變化Fig.3 Variation of lift and drag coefficient as a function of angle of attack

        圖4 hm=0.75和0.25時推力系數(shù)隨st的變化Fig.4 Variation of thrust force coefficient as a function of St number for hm=0.75 and 0.25

        (2)計算狀態(tài)為M=0.1,Re=40000,hm=0.25 和0.75,φ=90°,俯仰運動繞前緣1/3弦長點處轉動。圖4為翼型非定常運動隨St數(shù)變化產(chǎn)生的推進系數(shù)值同文獻及實驗值的比較。

        從上述的算例比較中可見,該程序不論是針對升力、力矩還是推力的計算,都能保證較高的數(shù)值精度,可用于接下來的推進性能研究。

        3.2 不同參數(shù)對撲翼推進特性的影響

        論文研究了包括沉浮運動振幅,沉浮運動減縮頻率,雙翼翼間距,沉浮俯仰運動相位差在內的多種參數(shù)對撲翼推進性能的影響,并對比了單翼和雙翼撲動在推進效果上的優(yōu)劣。

        流場參數(shù)為 Ma=0.1,Re=40000。為方便對比,雙翼運動產(chǎn)生的推力和所需要的輸入功率均除以2之后與單翼比較。故實際的雙翼推力和所需輸入功率均為給出對應結果的2倍。

        (1)沉浮運動振幅hm對撲翼推進的影響

        圖5為撲翼在不同沉浮運動振幅下推進性能的參數(shù)變化,其運動參數(shù)為k=1.5,L=0.8。圖5(a)可以看出當沉浮運動的振幅比較小時,無論是雙翼還是單翼振動都不會產(chǎn)生正推力,隨著振幅的增大,雙翼和單翼振動產(chǎn)生的推力都會增大,而且分別到達一個臨界點后開始由負推力轉變成為正推力,但雙翼出現(xiàn)正推力的幅值比單翼略小。由圖5(b)可以看出,隨著振幅的增大,撲翼所需要的輸入功率也越大,同時雙翼作動需要的功率要比單翼作動大。但由圖5(c)可以看出雙翼振動的推進效率絕大多數(shù)情況下都優(yōu)于單翼,而且效率的最大值比單翼大。

        圖5 沉浮振幅hm對撲翼推進特性的影響Fig.5 The effect of plunging amplitude hmon flapping airfoil propulsion

        (2)雙翼間距L對推進效果的影響

        圖6為雙翼翼間對推進特性的影響,其運動參數(shù)為k=1.5,hm=0.15。由于翼間距只是針對雙翼模型,所以這個計算條件下單翼產(chǎn)生的推進特性僅起到與雙翼對比的作用。由圖6(a、b)可以看出雙翼間距越小,振動產(chǎn)生推力越大,但需要的輸入功率也越大。由圖6(c)可以看出該狀態(tài)翼間距L在0.7附近推進效率最高。另外無論是推力,輸入功率還是推進效率,在同樣計算條件下雙翼模型都要大于單翼模型。

        圖6 沉浮運動雙翼翼間距L對推進特性的影響Fig.6 The effect of pacing between the dual airfoils L on flapping airfoil propulsion

        (3)減縮頻率k對推進特性的影響

        圖7為撲翼在不同減縮頻率k下作沉浮運動時的推進特性變化,其運動參數(shù)為L=0.8,hm=0.2。減縮頻率對推進的影響和沉浮振幅對推進的影響類似,圖7(a)可以看出當沉浮運動的頻率比較小時,無論是雙翼還是單翼振動不足以產(chǎn)生正推力,隨著頻率的增大,雙翼和單翼振動產(chǎn)生的推力都會增大,而且分別到達一個臨界頻率后開始由負推力轉變成為正推力。由圖7(b)可以看出,隨著頻率的增大,撲翼作動所需要的輸入功率也越大,同時雙翼作動需要的功率要比單翼大。由圖7(c)可以看出雙翼振動的推進效率基本上都要優(yōu)于單翼振動,而且效率的最大值比單翼要大。另外隨著減縮頻率的增大雙翼振動會比單翼振動更早到達效率最大點。

        圖7 減縮頻率k對沉浮推進特性的影響Fig.7 The effect of plunging reduced frequency k on flapping airfoil propulsion

        (4)相位差φ對沉浮俯仰組合推進效果的影響

        圖8、圖9為雙翼模型做沉浮俯仰運動,不同的相位差對推進效果的影響。選取的兩個計算狀態(tài)分別為L=0.8,k=1.5,hm=0.15,αm=5°和L=0.8,k =0.8,hm=0.15,αm=5°,圖中虛線給出的是相同狀態(tài)下雙翼僅作沉浮推進的結果。從圖8(a、c)和圖9 (a、c)中可以看出加入俯仰模態(tài)后,推力以及推進效率并不一定要優(yōu)于單獨的沉浮模態(tài),然而在某一的俯仰沉浮相角差下,推力系數(shù)和推進效率都可以到達一個比較高的值。另外,從圖8和圖9可見兩種不同減縮頻率下的運動達到最高推力和最高推進效率的相位角相差較大。

        圖8 k=1.5時沉浮俯仰運動相位差φ對推進特性的影響Fig.8 The effect of phase differenceφon flapping airfoil propulsion for k=1.5

        圖9 k=0.8時沉浮俯仰運動相位差φ對推進效果的影響Fig.9 The effect of phase differenceφon flapping airfoil propulsion for k=0.8

        3.3 對拍撲翼高效性能的流場機理分析

        圖10 撲翼的升力、阻力和力矩系數(shù)隨時間的變化Fig.10 Lift,drag and moment coefficient of flapping airfoils vs.time

        圖11 單翼和雙翼非定常運動的尾渦結構和流線關系Fig.11 Vortical pattern and streamline of single flapping airfoil and dual flapping airfoils

        為了從本質上說明雙翼對拍的優(yōu)越性,本文從非定常運動的流場特征進一步進行分析。圖10給出了一典型狀態(tài)下單翼和雙翼中上翼的氣動載荷響應結果。非定常運動狀態(tài)的參數(shù)為L=0.8,hm=0.2,k= 1.5。可以看出在雙翼對拍不僅產(chǎn)生的平均推力要大于單翼撲動,而且一個周期內,雙翼推力系數(shù)除了在上下兩個極限位置附件略微低于單翼外,其它行程中都大于單翼推力系數(shù)。在平衡位置附近,雙翼的推力明顯大于單翼的推力,其中通過平衡位置向上的瞬時點更明顯,此時雙翼的推力系數(shù)為0.598,約為單翼推進系數(shù)的1.27倍(單翼為0.472)。為了進一步揭示這一現(xiàn)象的物理機理,圖11給出了該瞬時(對應單翼和雙翼中的上翼恰好向上通過平衡位置)的流場特征。在此瞬時,雙翼之間的后緣尾跡區(qū)形成了一對明顯的旋渦(單翼沒有),該旋渦對是由于上下兩翼迅速遠離形成的低壓效應造成的尾緣回流。從圖12的壓力表面分布來看,雙翼對應的壓力分布較單翼更加飽滿(壓差更大),不僅增加了雙翼的升力幅值,而且也使得雙翼的推力較單翼大。

        圖12 撲翼瞬態(tài)的表面壓力分布Fig.12 Cpdistributions of flapping airfoils

        4 結 論

        本文以NACA0012翼型為例,研究了雙翼對拍時沉浮運動的振幅、頻率、雙翼的間距以及沉浮俯仰組合運動相角差等參數(shù)對推進性能的影響,得到了如下結論:

        (1)隨著沉浮振幅的增大,雙翼對拍產(chǎn)生的推力也會增大,逐漸會由負推力轉變成為正推力。同時隨著振幅的增大,撲翼所需要的輸入功率也越大,但總體推進效率較單翼振動更高。

        (2)雙翼翼間距越小,振動產(chǎn)生推力越大,需要的輸入功率也越大,但推進效率隨間距的減小不呈單調變化,在某一間距下推進效率達到最高。

        (3)雙翼和單翼的推力和所需輸入功率都隨沉浮頻率的增大而單調增加,但推進效率不呈單調變化,在某一頻率下達到最優(yōu)。雙翼的最大推進效率較單翼更高。

        (4)在合適的相位差下,俯仰和沉浮組合運動模式的推進特性優(yōu)于單純沉浮運動模式。另外,不同減縮頻率下的運動,其達到最大推力和最高推進效率的相位角也不相同。

        (5)對翼型非定常運動產(chǎn)生的流場特征進行分析,闡明了雙翼對拍推進性能較高的機理。

        [1] KNOLLER R.Die Gesetze des luftwiderstandes[J].Flug und Motortechnik(Wien),1909,3(21):1-7.

        [2] BETZ A.Ein beitrag zur erklaerung des segelfluges[J].Zeitschrift fuer Flugtechnik und Motorluftschiffahrt,1912,3:269-272.

        [3] KATZMAYR R.Effect of periodic changes of angle of attack on behavior of airfoils[R].NACA Report.Oct,1922.

        [4] Von KARMAN T,BURGERS J M.General aerodynamic theory-perfect fluids[M].Division E,Vol.II,Aerodynamic Theory,Ed.Durand,W.F.,1943.308.

        [5] GARRICK I E.Propulsion of a flapping and oscillating airfoil [R].NACA Technical Report,1936,(567).

        [6] JONES K D,PLATZER M F.Numerical computation of flapping-wing propulsion and power extraction[R].AIAA 97-0826,1997.

        [7] TONG B G,LU X Y.A review on biomechanics of animal flight and swimming[J].Advances in Mechanics,2004,34(1):1-7.(in Chinese)童秉綱,陸夕云.關于飛行和游動的生物力學研究[J].力學進展,2004,34(1):1-7.

        [8] LU X Y,YANG J M,YIN X Z.Some studies on biofluiddynamics of flight and swimming animals and biomimetic technology[J].Journal of University of Science and Technology of China,2007,37(10):1159-1163.(in Chinese)陸夕云,楊基明,尹協(xié)振.飛行和游動的生物運動力學和仿生技術研究[J].中國科學技術大學學報,2007,37(10):1159-1163.

        [9] YANG L.Research on hydrodynamic characteristics and propulsive mechanism of oscillating tuna-tail[D].Harbi:Harbin Engineering University,2009.(in Chinese)楊亮.仿金槍魚擺動尾鰭的水動力性能與推進機理研究[D].哈爾濱:哈爾濱工程大學,2009.

        [10]YANG S,LUO S,LIU F.Computation of the flows over flapping airfoil by the euler equations[A].The AIAA 43rd Aerospace Sciences Meeting,Reno[C].NV,10-13,2005.

        [11]WANG T.Performance investigation of symmetrical plunging dual-foils propulsor[J].Journal of Hydrodynamics,2009,A24(5):647-653.(in Chinese)王濤.對稱雙翼沉浮振動推進性能研究[J].水動力學研究與進展,2009,A24(5):647-653.

        [12]YU C J,ANG H S.Numerical study of aerodynamics for flexible membrane flapping-wing MAV[J].Journal of University of Science and Technology of China,2009,39(12):1305-1310.(in Chinese)余春錦,昂海松.柔性膜微型撲翼飛行器氣動力的數(shù)值研究[J].中國科學技術大學學報,2009,39(12):1305-1310.

        [13]YANG W Q,SONG B F,SONG W P.Aerodynamic performance research of micro flapping-wing in low Reynolds number flow[J].ACTA Aerodynamica Sinica,2011,29(1):32-38.(in Chinese)楊文青,宋筆鋒,宋文萍,等.微型撲翼低雷諾數(shù)繞流氣動特性研究[J].空氣動力學學報,2011,29(1):32-38.

        [14]QIAN J,ZHANG Z,LUO S,et al.Numerical study of flow characteristics of a plunging rigid airfoil with elastic trailing-edge plate[J].ACTA Aerodynamica Sinica,2012,30(1):113-119.(in Chinese)錢靖,張正科,羅時鈞,等.二維彈性撲翼沉浮運動流動特性的數(shù)值研究[J].空氣動力學學報,2012,30(1):113-119.

        [15]LIU P.Propulsive performance of a twin-rectangular-foil propulsor in a counter-phase oscillation[J].Journal of Ship Research,2005,49(3):207-215.

        [16]LIU P,WANG L,HUANG G,et al.Propulsion characteristics of wing-in-ground effect dual-foil propulsors[J].Applied Ocean Research,2010,32(1):103-112.

        [17]JIANG Y W,YE Z Y.A cell-centered finite volume method for arbitrary grid type[J].Chinese Journal of Theoretical and Applied Mechanics,2010,42(5):830-837.(in Chinese)蔣躍文,葉正寅.適用于任意網(wǎng)格拓撲和質量的格心有限體積法[J].力學學報,2010,42(5):830-837.

        Analysis of propulsive performance for dual flapping airfoils using CFD method

        ZHANG Weiwei,XU Yang,JIANG Yuewen,YE Zhengyin

        (National Key Laboratory of Aerodynamic Design and Research,Northwestern Polytechnical University,Xi′an 710072,China)

        Dual flapping airfoils propulsion has many advantages than a single flapping airfoil,such as higher propulsive efficiency,easier mechanical design,more stationary lift and inertial force.The unsteady flow fields around flapping airfoils are simulated by solving Reynolds Averaged Navier-Stokes equations based on dynamically deformable hybrid grids.Some flapping parameters,such as plunging amplitude,plunging reduced frequency,spacing between the dual airfoils,phase difference between plunge and pitch of composite motion on thrust force coefficient,and propulsive efficiency,are studied by simulating the unsteady low speeds flow over dual flapping NACA0012 airfoils.It is found that dual flapping airfoils achieve a higher thrust force coefient and propulsive efficiency compared with a single flapping airfoil.It is also shown that there is an optimum phase angle of composite motion for a given flapping motion mode.The thrust force coefficient and propulsive efficiency of the composite motion are higher than those of the pure plunging motion.

        flapping airfoil propulsion;dual flapping airfoil;numerical simulation;unsteady flow field

        V211.3

        Adoi:10.7638/kqdlxxb-2012.0148

        0258-1825(2014)04-0446-07

        2012-09-10;

        2013-02-05

        航空基金(20121353014);西北工業(yè)大學基礎研究基金

        張偉偉(1979-),男,博士,教授,博士生導師,主要研究方向:氣動彈性力學、非定??諝鈩恿W、流動控制、葉輪機流固耦合力學.E-mail:aeroelastic@nwpu.edu.cn

        張偉偉,徐楊,蔣躍文,等.基于CFD方法的對拍撲翼推進特性分析[J].空氣動力學學報,2014,32(4):446-452.

        10.7638/kqdlxxb-2012.0148. ZHANG W W,XU Y,JIANG Y W,et al.Analysis of propulsive performance for dual flapping airfoils using CFD method[J].ACTA Aerodynamica Sinica,2014,32(4):446-452.

        猜你喜歡
        振動效率研究
        振動的思考
        科學大眾(2023年17期)2023-10-26 07:39:14
        FMS與YBT相關性的實證研究
        遼代千人邑研究述論
        提升朗讀教學效率的幾點思考
        甘肅教育(2020年14期)2020-09-11 07:57:42
        振動與頻率
        天天愛科學(2020年6期)2020-09-10 07:22:44
        視錯覺在平面設計中的應用與研究
        科技傳播(2019年22期)2020-01-14 03:06:54
        EMA伺服控制系統(tǒng)研究
        中立型Emden-Fowler微分方程的振動性
        跟蹤導練(一)2
        “錢”、“事”脫節(jié)效率低
        久久无码人妻一区二区三区午夜| 91精品国产综合成人| 国产亚洲sss在线观看| 日本黑人人妻一区二区水多多| 偷拍一区二区三区高清视频| 香港三日本三级少妇三级视频| 国模无码视频一区| 欧美精品高清在线xxxx| 亚洲一区免费视频看看| 边添小泬边狠狠躁视频| 日本公妇在线观看中文版| 狠狠躁天天躁无码中文字幕图| 国产精品一区二区三区黄片视频 | 欧美在线专区| 中文字幕人妻少妇久久| 男人的天堂一区二av| 国产肥熟女视频一区二区三区 | 亚洲精品国产美女久久久| 一区二区av日韩免费| 国产精品白浆一区二区免费看| 国产精品亚洲lv粉色| 在线精品免费观看| 日本中文字幕一区二区在线观看 | 岛国熟女一区二区三区| 国产日产韩国级片网站| 比较有韵味的熟妇无码| 国产日韩欧美亚洲精品中字| 夫妻一起自拍内射小视频| 在线视频色系中文字幕| 一区二区三区人妻无码| 亚州精品无码人妻久久| 一区二区在线观看视频亚洲| 婷婷色综合视频在线观看| 日本老熟欧美老熟妇| 国产精品无码久久AⅤ人妖| 在线播放草猛免费视频| 国产精品美女久久久久久| 99久久久精品免费香蕉| 一区二区三区在线观看精品视频| 欧美最猛性xxxx| 亚洲精品欧美二区三区中文字幕|