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

無人機果樹施藥旋翼下洗氣流場分布特征研究

2019-11-08 01:07:20祁力鈞吳亞壘程湞湞劉婠婠ElizabethMusiu楊澤鵬
農業工程學報 2019年18期
關鍵詞:風速

張 豪,祁力鈞,吳亞壘,程湞湞,劉婠婠,Elizabeth Musiu,肖 雨,楊澤鵬

無人機果樹施藥旋翼下洗氣流場分布特征研究

張 豪,祁力鈞※,吳亞壘,程湞湞,劉婠婠,Elizabeth Musiu,肖 雨,楊澤鵬

(中國農業大學工學院,北京 100083)

植保無人機懸停果樹施藥時的旋翼下洗氣流場分布對霧滴空間運動和在冠層內部的附著、穿透有重要影響。該文基于計算流體動力學(computational fluid dynamic,CFD)方法,結合RNG κ–ε湍流模型、多孔介質模型和滑移網格技術,通過構建虛擬果園,對六旋翼植保無人機懸停果樹施藥時的下洗氣流流場進行數值模擬,分析在無人機不同懸停高度、不同果樹生長階段和不同自然風速下的氣流場分布特征,并進行標記點下洗氣流速度測試試驗。研究結果表明:1)自然風速大于3 m/s時,旋翼下洗氣流速度已淹沒于環境自然風速中,不再滿足植保無人機懸停施藥作業條件;2)自然風破壞了旋翼下洗氣流的中心對稱狀態,向下風方向出現后揚,且隨著自然風速和懸停高度的增大,后揚距離隨之增大;3)與無自然風狀態比較,果樹生長時期對其噴頭處速度分布影響不顯著,主要受自然風影響,且豎直向下的向氣流占主體地位,對霧滴的對靶運輸起主導作用,應將噴頭安裝于可使霧滴獲得較大向速度的旋翼正下方0.2 m處附近;4)無人機懸停位置沿逆風方向調整后,冠層內部上、中、下層氣流平均速度較調整前分別由1.36、0.80、0.81 m/s增大至3.04、2.37、1.63 m/s;上、下層速度分布變異系數分別由74.26%、35.80%降至45.39%和22.70%,中層略有增大,總體利于實現對靶噴霧。試驗結果表明,標記點下洗氣流速度測量值和模擬值之間具有較好的一致性。該文可為動態環境條件下植保無人機懸停果樹施藥的對靶噴霧自適應控制技術研究提供參考。

無人機;農藥;模型;下洗氣流;果樹冠層;多孔介質

0 引 言

農業航空精準施藥技術是減少農藥殘留、提高農藥效益及降低農藥對環境負面影響的有效手段[1]。作為現代農業的重要組成部分,農業航空植保具有效率高、成本低、作業不受地形條件和作物長勢限制等優勢[2–6]。但在旋翼植保無人機懸停狀態下進行果樹施藥時,果樹冠層對旋翼下洗氣流的影響不可忽略[7],且旋翼下洗氣流裹挾霧滴運動,影響霧滴的空間運動及在冠層內部的附著和穿透。為研究旋翼下洗氣流與作物冠層之間的互作關系[8],有必要明確果樹冠層參與下的旋翼下洗氣流場分布特征。

目前,圍繞旋翼植保無人機作業參數對霧滴飄移[9-10]、沉積[11-12]和噴幅[13-14]的影響已有較多研究。另外,Bruno等[15-16]利用無人機接受無線傳感器網絡(wireless sensor networks,WSNs)的回傳信息,實現無人機航線在動態環境中自主調整,達到精準對靶施藥的目的。而在旋翼下洗氣流方面,Ramasamy等[17]、Wall等[18]采用粒子圖像測速(particle image velocimetry,PIV)方法對翼尖渦流形態特征進行了研究;Yeo等[19]建立了旋翼無人機下洗氣流的經驗模型,并用3種機型對模型可靠性進行了評估;汪沛等[20]、胡煉等[21]、李繼宇等[22]采用風場測量系統分別針對不同型號植保無人機旋翼下洗氣流進行田間測量,通過分析采集的風速數據得到風場分布規律,但都未能直觀反映出作物冠層參與下的旋翼下洗氣流三維形態。陳盛德等[23]則用該系統測得四旋翼無人機下洗氣流分布情況,研究了不同方向風速對霧滴沉積分布的影響;Yang等[24-25]采用計算流體動力學(computational fluid dynamic,CFD)方法,對六旋翼植保無人機的空載懸停下洗氣流進行數值模擬,揭示了無冠層的旋翼下洗氣流分布規律,并在載荷3 kg數值求解的下洗氣流場中添加霧滴離散相,初步研究了旋翼下洗氣流與霧滴運動的內在聯系;文晟等[26]采用格子玻爾茲曼方法(lattice-boltzman method,LBM)對單旋翼植保無人機旋翼流場進行數值求解,探尋了不同飛行速度下翼尖渦流對霧滴飄移的影響規律。數值模擬過程中對果樹冠層的處理方法有2種:一是創建完整的3D果樹模型(包括樹干、樹枝及樹葉等);二是用多孔介質代替冠層[27]。前者需要較大計算機資源,不宜開展數值模擬研究,后者在果園風送式噴霧機數值模擬方面已有較多成熟研究,如Endalew等[28-30]、Hong等[31-32]、Duga等[33-34]、Salcedo等[35-36]在模擬過程中引入多孔介質模型處理果樹冠層,廣泛開展有冠層參與下的噴霧機氣流場和霧滴分布規律研究。Hong等[37]將果樹冠層用多孔介質代替進行CFD模擬,利用霧滴飄移模擬數據開發了風送式噴霧機模擬(simulation of air–assisted sprayers,SAAS)應用軟件,該軟件可根據用戶輸入條件對噴霧機霧滴飄移做出預測。雖然植物冠層內氣流的不均勻和湍流性質、冠層葉面積密度參數誤差及考慮霧滴時氣流和靶標冠層之間相互作用建模困難等因素,限制了采用多孔介質代替冠層后的預測精度[27]。但上述研究已充分證明利用多孔介質模型處理果樹冠層,進行作物冠層參與下的氣流場數值模擬研究是準確可行的。總的來看,基于多孔介質模型對果樹冠層參與下的旋翼下洗氣流場分布特征的研究報道較少。

為明晰多旋翼植保無人機懸停果樹施藥時的氣流場分布特征,本文通過構建虛擬果園,采用多孔介質模型處理果樹冠層,基于商用軟件ANSYS Fluent16.0對六旋翼植保無人機懸停果樹施藥氣流場進行數值模擬。分析不同懸停高度、不同果樹生長階段和不同自然風速下的互作流場分布特征,并搭建試驗平臺開展標記點下洗氣流速度驗證試驗。

1 物理模型

六旋翼植保無人機定制于深圳市金銘睿電子有限公司(圖1),因初步只對旋翼下洗氣流進行研究,所以該植保無人機未設計噴霧系統。無刷電機由外接電源提供動力,旋翼轉速由無線遙控器進行調控,另外加裝有RC41轉速測量儀,反饋旋翼轉速用于數值模擬,連接桿用于將無人機懸掛固定于支架上。無人機基本參數為:軸距0.8 m,旋翼型號1 555,無刷電機型號X4114 KV370,最大載藥量5 kg。

1.轉速測量儀 2.連接桿 3.無刷電機 4.旋翼

植保無人機旋翼是核心旋轉部件,其外表面為復雜曲面,直接進行正向精準建模比較困難,且六旋翼植保無人機的正、反旋翼相間分布,需分別對二者進行建模。為獲得精確旋翼三維模型用于數值模擬,本文采用MCS五四軸全自動三維掃描系統分別對正、反旋翼進行掃描,獲取了旋翼表面點云數據,而后采用Geomagic Studio軟件對點云數據進行后處理,實現旋翼模型的逆向重建,得到精確的旋翼三維模型(圖2)。

1.旋翼托架 2.三維掃描儀 3.旋翼 4.旋轉托盤 5.點云數據 6.計算機 7.旋翼重構模型

數值模擬過程中將無人機模型進行合理簡化,忽略無人機起落架、機臂及噴霧系統等其他部件影響,只考慮各旋翼的存在,以使計算機資源與數值求解準確性達到平衡。噴霧系統的忽略可能對原藥箱所在處的氣流場造成影響,但旋翼距藥箱的水平距離一般較遠,不影響氣流的向下發展,且已有研究發現這種簡化對旋翼下洗氣流的整體發展演變影響不大[24],簡化后模型可用于旋翼氣流場的數值模擬研究。簡化后模型如圖3所示,定義坐標系軸正方向為無人機下降方向,旋翼1、3、5為反旋翼(逆時針旋轉),旋翼2、4、6為正旋翼(順時針旋轉),平面與6個旋翼所在平面重合,坐標系原點位于6個旋翼旋轉軸所在圓的圓心處。

圖3 無人機簡化模型

2 數值模擬

2.1 主控方程與湍流模型

納維-斯托克斯(Navier-Stokes,N-S)方程是描述流體運動的控制方程,但運用三維非定常N-S方程對旋翼下洗氣流進行直接數值求解尚難以實現。應用雷諾平均法在N-S方程中引入湍流模型將方程組(RANS方程)封閉求解,計算量小,更適用于工程數值計算[38],在旋轉坐標系下以絕對速度為變量的主控方程具體形式可表示為[39]

式中分別為控制體邊界面上的對流通量和黏性通量;為狀態變量;為流體的密度,kg/m3;、、分別表示速度矢量在3個坐標方向(、、)上的分量,m/s;為單位體積的總能量;W為控制體的體積,m3;為控制體邊界法向面;為旋轉角速度矢量,rad/s;為、、的矢量和;為邊界的外法向3個分量(、、)的矢量和,m/s;W為邊界面上網格運動速度,m/s;為壓力,Pa;為時間,s;為黏性剪切應力張量項(為取、、時的組合);為熱傳導項。

本文選擇RNG-湍流模型對主控方程進行封閉。該模型對湍動黏度進行了修正,并且考慮了平均流動中的旋轉及旋流流動情況[40],具體執行公式為

式中為湍動能,m2/s2;為湍動耗散率,m2/s3;為有效黏度,Pa×s;為流體動力黏度,Pa×s;uu為速度分量時均值,m/s;G為由于平均速度梯度引起的湍動能的產生項,m2/s2;1和2為模型常數,1=1.42,2=1.68;分別為和的有效湍流Prandtl數的倒數,==1.39;為平均應變率;0、均為默認常數。

2.2 冠層多孔介質模型

果樹冠層的阻力作用會造成氣流動量損失,本文在總結前人研究成果[27-37]基礎上選擇多孔介質模型處理果樹冠層,將冠層所在區域用多孔介質代替,通過在多孔介質區域中增加動量損失源項來模擬冠層對氣流的阻力作用。動量、湍動能和湍動耗散率方程的附加源項可描述為[31]

已知冠層壓力損失系數和葉面積密度之間的關系為C=2CL,可通過定義不同壓力損失系數表征冠層疏密程度。基本假設有:多孔介質為各向同性;多孔介質流動阻力忽略黏性項影響,只保留慣性損失項。

2.3 計算區域與邊界條件

由于針對旋翼下洗氣流進行田間試驗研究不僅耗時費力、精確測量設備成本高,而且自然風速等試驗條件極不可控,因此,本文通過構建虛擬果園,對六旋翼植保無人機懸停果樹施藥時的氣流場進行數值模擬研究。構建的虛擬果園(圖4)長20.0 m,寬18.0 m,行間距和株間距分別被設定為3.0和2.0 m。果樹模型被簡化為球形冠層和圓柱樹干的組合[31]。考慮標準化果園里果樹結構已生長完全,且冠層常年被修剪為固定大小形狀,所以對不同生長階段(發芽期、半葉期和全葉期)的果樹只考慮冠層壓力損失系數的變化,不考慮結構尺寸的變化。結構尺寸由測量實驗室內的試驗果樹尺寸得到。如圖5所示,果樹冠層球直徑為1.2 m,樹高1.9 m,圓柱樹干直徑為0.15 m。果樹發芽期、半葉期和全葉期的葉面積密度分別為1.6、3.0、8.0 m–1,對應冠層壓力損失系數分別為0.8、1.5、4.0[32]。若對整個果園進行模擬,需較高計算機資源,求解耗費時間長。為提高計算效率,在不影響旋翼下洗氣流充分發展擴散的前提下,選擇虛擬果園的一部分作為數值求解計算區域。

圖4 虛擬果園示意圖

圖5 果樹模型參數

如圖6所示,計算區域長10 m,寬9 m,考慮本文無人機最高懸停高度為4 m,計算區域高度設定為6 m,植保無人機懸停于靶標果樹正上方進行懸停施藥作業。包含無人機旋翼的6個旋轉子域、果樹的15個冠層子域和15個樹干子域、1個空氣子域。利用滑移網格技術處理旋翼的旋轉(3 000 r/min),在旋轉子域和空氣子域交界面處采用動態搭接面邊界條件;由于樹干子域不需求解,在前處理過程中采用布爾差運算將其從計算區域中除去,只保留樹干外邊界,設置為壁面邊界。沿軸正方向,自左向右依次為自然風進口和出口,進口邊界為速度入口,出口邊界為壓力出口,出口壓力為0。地面設置為壁面邊界,其余邊界為對稱邊界。計算區域網格劃分的質量好壞直接影響數值求解的準確性和計算效率,本文采用ANSYS Meshing中適用于復雜實體的非結構四面體網格對計算區域進行劃分。通過控制旋翼壁面尺寸、旋轉子域與冠層子域邊界面尺寸對該2個子域的網格進行加密,提高計算精度。劃分后計算區域網格節點總數約192萬,網格單元總數約1 032萬。對于四面體網格,可采用偏斜率(Skewness)衡量網格質量,查看偏斜率最大值為0.85,平均值為0.23,符合四面體網格最大偏斜率應低于0.95,平均偏斜率應低于0.33的要求[32],可用于數值模擬。

圖6 計算區域示意圖

2.4 計算方法

基于有限體積法將控制方程離散化求解瞬態氣流場,壓力–速度耦合采用適用于非定常計算并對壓力進行2次修正的PISO算法,空間域上選擇二階迎風離散格式對動量、湍動能和湍動耗散率進行離散,用該離散格式對非結構四面體網格進行求解可達到較高求解精度,時間域上選擇一階隱式離散格式,壓力插值格式選擇適合高速旋流及多孔介質的PRESTO!格式[34]。

3 模擬結果與分析

3.1 旋翼下洗氣流場分布特征

為明晰多旋翼植保無人機懸停果樹施藥時氣流場的分布情況,首先選取無人機懸停高度為3 m,自然風速為0 m/s情況下的數值模擬結果進行分析。由圖7可知,在不考慮自然風的情況下,該結果與單一果樹的研究結果一致[7],具體表現為靶標果樹冠層周圍氣流從冠層上半部區域開始呈“圓錐形”向下發展,以一傾斜角發展到地面之后形成地面鋪展,并且隨著果樹葉面積密度增大,冠層對旋翼下洗氣流向四周的擴散作用增強。該擴散作用將增加靶標果樹冠層周圍霧滴的橫向飄移,最終造成霧滴的地面無效沉積,不利于無人機懸停對靶精準噴霧。

圖7 各時期XOZ截面總速度分布云圖

植保無人機進行田間施藥作業,噴霧質量受噴灑系統、農藥理化特性及自然風等諸多因素影響[5],其中自然風是造成霧滴發生遠距離飄移的主要原因之一。為研究自然風動態環境條件下多旋翼植保無人機懸停果樹施藥時氣流場的具體形態,以果樹半葉期為例,對不同自然風速下的氣流場進行數值模擬。圖8給出了無人機懸停高度為3 m,自然風速為1~5 m/s狀態下的總速度三維分布云圖,三維形態代表渦的分布情況,顏色代表總速度分布情況。從圖8知,當自然風(指向軸正向)存在時,旋翼下洗氣流不再是豎直向下發展,而是在自然風的作用下向下風方向出現后揚。此時,旋翼下洗氣流已不能將靶標果樹冠層完全包裹。隨著自然風速的增大,旋翼下洗氣流后揚距離隨之增大,這將增加霧滴沿下風方向的飄移距離,也將造成靶標果樹的漏噴以及相鄰已噴果樹的重噴。當自然風速小于等于3 m/s時,仍可觀察到部分后揚的旋翼下洗氣流發展至地面,且有完整的渦分布。但隨著自然風速的增大,旋翼氣流場渦形態的三維分布被破壞。當自然風速大于3 m/s時,旋翼下洗氣流速度已淹沒于環境自然風速中,此時已不再適宜開展施藥作業,且無人機飛行安全也將受到威脅。

圖8 不同自然風速狀態下總速度三維分布云圖(半葉期)

為觀察自然風對近地面旋翼下洗氣流的影響,圖9給出了果樹半葉期,無人機懸停高度為3 m,自然風速為1~5 m/s狀態下的截面總速度分布云圖。無自然風時,旋翼下洗氣流發展至地面形成對稱的地面鋪展(圖7),當自然風存在時,其不再呈對稱分布狀態(圖9)。自然風上行側近地面鋪展氣流被削弱,同時被逆向近地面自然風吹起,在自然風來流與上行側地面鋪展氣流的共同作用下形成了較大順時針卷揚氣流,而自然風下行側的近地面鋪展氣流則被加強(圖9)。當自然風速較小時,近地面卷揚氣流只出現在地面鋪展氣流末端(圖9a)。當自然風速增大到2 m/s時,卷揚氣流被向前推至靶標果樹上行側邊緣,且卷揚強度增大,該卷揚氣流增大了靶標果樹上行側氣流的湍流程度(圖9b)。當自然風速大于等于3 m/s時,已不再出現上行側地面鋪展氣流,旋翼下洗氣流隨自然風向下行側發展,且上述卷揚現象消失。

圖9 不同自然風速狀態下XOZ截面總速度分布云圖(半葉期)

3.2 氣流后揚距離因素分析

為探究自然風和無人機懸停高度對旋翼下洗氣流后揚距離的影響,以果樹全葉期為例,進行2個案例的數值模擬。案例1:懸停高度3 m,自然風速1 m/s、2 m/s、3 m/s;案例2:自然風速2 m/s,懸停高度2.5 m、3 m、3.5 m、4 m。因案例存在1組交叉組,故只需模擬6組。

對比圖9和圖10可知,有自然風時,全葉期和半葉期的旋翼下洗氣流后揚形態基本相同。懸停高度3 m狀態下,自然風速1和2 m/s時旋翼下洗氣流后揚距離分別達到1和2 m;自然風速為3 m/s時,后揚距離已超過2 m,氣流拍打到下行側相鄰果樹冠層(圖10);其中,自然風速為1 m/s時,對比圖9a和圖10a可知,在自然風下行側的靶標果樹冠層下部氣流形態存在不同,這是因為全葉期果樹冠層葉面積密度較半葉期大,導致氣流動量損失也較大。該現象在圖11a中更為明顯。自然風速2 m/s狀態下,懸停高度2.5 m、3 m和3.5 m時旋翼下洗氣流后揚距離相差不大,都為2 m。但懸停高度為2.5 m時,到靶氣流速度最大達8.4 m/s,易對枝葉造成損傷且不利于霧滴在冠層頂端葉面上的沉積;懸停高度3 m和3.5 m時,前者與靶標果樹冠層存在接觸,后者并不接觸,這是因為懸停高度越高,旋翼下洗氣流到達冠層高度時的氣流速度越小,抵抗自然風來流能力變弱,更易被自然風吹向下行側,后揚旋翼下洗氣流變得更細長;懸停高度4 m時,后揚距離已超過2 m(圖11)。

圖10 不同自然風速狀態下XOZ截面總速度分布云圖(全葉期)

圖11 不同高度狀態下XOZ截面總速度分布云圖(全葉期)

綜合比較,懸停高度為3 m(距離冠層高度在1 m左右)時開展植保無人機懸停果樹施藥作業效果更好。另外,考慮到在自然風動態環境條件下存在旋翼下洗氣流后揚現象,可根據自然風速的大小對植保無人機懸停位置向逆風方向做出適當調整,以適應旋翼下洗氣流包裹靶標果樹冠層,迫使霧滴盡可能多的到達冠層,從而實現精準對靶噴霧。

3.3 噴頭處氣流速度分布

為進一步研究旋翼下洗氣流后揚情況下噴頭不同安裝距離范圍內氣流的速度分布,針對植保無人機懸停高度3 m、自然風速0~2 m/s狀態下進行施藥作業時,在旋翼2、旋翼3、旋翼5和旋翼6的正下方安裝噴頭,并分別放置速度檢測線。檢測線按旋翼編號命名。考慮無刷電機和噴頭固定裝置需要一定距離及噴頭安裝不可能無限制靠下,選擇旋翼正下方=0.1~0.6 m范圍進行分析。

分析果樹各時期不同自然風速狀態下檢測線L2、L3、L5和L6的、、方向速度可知,同一自然風速狀態下,不同時期(發芽期、半葉期和全葉期)的速度變化趨勢基本一致,表明果樹生長時期對同一自然風速下的檢測線各方向速度分布無顯著影響。而在同一時期,不同自然風速狀態下檢測線的、、方向速度發生變化,以下行側檢測線L2、L3的向速度變化最為明顯,表明檢測線各方向速度變化主要受自然風影響。

以果樹發芽期為例,從圖12知,自然風速為0 m/s時,4條檢測線的、、方向速度變化趨勢基本相同,表明各檢測線處氣流速度分布較為一致,呈中心對稱狀態,且向速度較、向速度大,說明向(豎直向下)氣流占主體地位,這與已有研究結論一致[24]。表明向氣流將對霧滴的對靶運輸起主導作用,最大向速度近9 m/s。當存在自然風時,由于氣流出現后揚,造成檢測線L2、L3與檢測線L5、L6的向速度變化趨勢不同。L2、L3的向速度與無自然風時較為一致,而L5、L6的向速度同步減小。與自然風速1 m/s時比較,自然風速為2 m/s時L5、L6的向速度衰減更快。隨著自然風速的增大,各檢測線向速度值均逐漸呈現正值,表明自然風的存在破壞了旋翼下洗氣流的中心對稱狀態,旋翼下洗氣流向下行側的分速度增大,這將迫使霧滴向下行側運動,造成靶標果樹冠層上行側與下行側的霧滴沉積分布不均勻。但由于自然風為單一風向,、、方向速度仍以截面為對稱面呈對稱分布。綜合比較,各狀態下檢測線L2、L3的最大向速度分布在旋翼正下放0.2~0.3 m區域內,當安裝距離超過0.2 m時,L5、L6的向速度減小明顯,為使各噴頭霧滴均獲得較大向速度,應將噴頭安裝于旋翼正下方0.2 m附近。

注:“La-Vb”表示檢測線La的b方向速度(a取2、3、5或6;b取x、y或z)

3.4 靶標果樹冠層氣流速度分布

圖13為無人機懸停高度3 m,自然風速0~2 m/s狀態下到靶氣流總速度分布云圖。從圖13知,無自然風時,旋翼下洗氣流以中心對稱形態到達冠層表面,對冠層實現包裹,且存在6個均布的高速區。該到靶特征有助于對靶噴霧的實現,可將霧滴均勻的吹送至冠層內部。由于自然風的存在,到靶氣流不再呈中心對稱形態,風速越大,到靶區域越偏向下行側,且到靶區域面積越小,這將顯著降低霧滴到達冠層的幾率。

從圖14知,無自然風時,冠層內部氣流速度呈對稱分布且相對較為均勻,該分布狀態有助于霧滴在冠層內部的均勻擴散。自然風的存在,造成冠層內部氣流向下行側出現偏斜,且隨著風速的增大氣流偏斜越嚴重,該不均勻分布狀態將造成冠層內部霧滴的不均勻分布,呈現為冠層下行側沉積相對較多的趨勢。

圖13 不同自然風速狀態下到靶氣流總速度分布云圖(全葉期)

圖14 不同自然風速狀態下冠層內XOZ截面總速度分布云圖(全葉期)

為進一步探究自然風動態環境條件下,調整無人機懸停位置對實現對靶噴霧的效果,將圖13c無人機懸停的水平位置沿逆風方向調整1 m,模擬結果見圖15。從圖15知,調整后旋翼下洗氣流從冠層上行側開始對冠層實現包裹,且到靶區域面積、冠層內部氣流分布范圍及氣流速度大小較調整前均有所增大,這將有助于霧滴向冠層下行側穿透,增加靶標果樹冠層內部的有效沉積,降低空中飄移和地面無效沉積。

為量化比較無人機調整前后靶標果樹冠層內部速度分布的均勻性,將冠層分為上(=1.4 m)、中(=1.7 m)、下(=2.0 m)3層,在各層分別均勻取9個樣本點,以速度分布變異系數作為評價指標,計算公式為[41]

由表1知,調整前冠層內部上、中、下層的平均速度分別為1.36、0.80和0.81 m/s,調整后冠層內部上、中、下層的平均速度分別為3.04、2.37和1.63 m/s,較調整前增大明顯,這有助于氣流攜帶霧滴進入冠層內部,且有助于樹葉的翻動,促使葉片背面附著霧滴。調整前后冠層內部上、中、下層的速度分布變異系數發生變化。上層速度分布變異系數顯著降低,從74.26%變為45.39%;中層從17.50%增大至28.69%,此處氣流速度的不均勻程度雖有增大,但相比于調整前中層氣流平均速度很小,在冠層內部達不到運輸霧滴的效果,調整后中層的平均速度增大明顯,反而有助于運輸霧滴;下層速度分布變異系數從35.80%降至22.70%。

表1 無人機位置調整前后冠層內部各層速度分布對比

4 驗證試驗

為驗證數值模擬的可靠性,于中國農業大學植保機械實驗室搭建試驗平臺,將數值模擬所用物理模型六旋翼植保無人機固定于支架上,果樹放置于無人機正下方,進行標記點向下洗氣流速度測試試驗,室內無自然風。在檢測線L2、L3、L5和L6上各選取4個標記點(=0.2 m、=0.3 m、=0.4 m、=0.5 m),采用testo 405i熱線式風速測量儀對標記點進行測量。風速儀在手機端配有Smart Probe APP,每2 s記錄一次數據,通過藍牙實現數據的無線傳輸及實時顯示。調整旋翼轉速至3 000 r/min,待旋翼下洗氣流發展穩定后開始測量,各標記點持續測量30 s,導出數據,求脈動速度的平均值為標記點速度值。測量現場如圖16所示,測量結果及相對誤差見表2。

圖16 氣流速度測量現場

表2 標記點模擬值和測量值對比

以均方根誤差(root mean square error, RMSE)和平均相對誤差(average relative error)來比較測量值和模擬值,計算公式為[25]

由表2知,各標記點測量值和模擬值的相對誤差均保持在10%以內,計算得到均方根誤差和平均相對誤差分別為0.48 m/s和6.24%,證明測量值和模擬值之間具有較好的一致性。某些標記點的相對誤差接近10%,可能是因為實際測量過程中,不能完全保證風速儀的熱線探頭準確置于旋翼正下方與檢測線重合,測量點位置存在一定偏差導致。

作為動態環境條件下無人機懸停果樹施藥對靶噴霧自適應控制技術的前期研究,本文尚未涉及霧滴離散相,自然風向僅為單一方向,且只對無自然風狀態下的氣流速度進行驗證試驗,研究仍存在不足。后續將以無人機果樹施藥田間作業環境為基礎,增加霧滴因素開展進一步研究,為自適應控制技術的研究提供參考。

本文后續擬建立動態環境條件下霧滴飄移預測模型,開展動態環境條件下無人機懸停果樹施藥的對靶噴霧自適應控制技術研究。

5 結 論

1)多旋翼植保無人機懸停狀態下進行果樹施藥,無自然風時,隨果樹葉面積密度增大,冠層對旋翼下洗氣流向四周的擴散作用增強;自然風存在時,旋翼下洗氣流中心對稱狀態被破壞,氣流向下風方向出現后揚,且隨自然風速和懸停高度的增大,氣流后揚距離隨之增大。各標記點下洗氣流速度測量值和模擬值相對誤差均在10%以內,均方根誤差和平均相對誤差分別為0.48 m/s和6.24%,測量值和模擬值之間具有較好的一致性。

2)當自然風速小于等于2 m/s時,靶標果樹上行側出現較大順時針卷揚氣流;當自然風速大于3 m/s時,旋翼下洗氣流速度已淹沒于環境自然風速中,不再適宜開展植保無人機懸停果樹施藥作業。

3)通過在噴頭安裝范圍內布置速度檢測線,分析表明,與無自然風狀態相比,果樹生長時期對噴頭處氣流速度的分布無顯著影響,自然風影響較大。并且,豎直向下的向氣流占主體地位,對霧滴的對靶運輸起主導作用,應將噴頭安裝于可使霧滴獲得較大向速度的旋翼正下方0.2 m處附近。

4)無人機懸停位置沿逆風方向調整后,旋翼下洗氣流從靶標果樹冠層上行側開始對冠層實現包裹,且冠層內部上、中、下層氣流平均速度較調整前增大明顯。分別由1.36、0.80、0.81 m/s增大至3.04、2.37、1.63 m/s。上、下層速度分布變異系數分別由74.26%、35.80%降至45.39%和22.70%,中層略有增大,但總體利于實現對靶噴霧。

[1] Lan Yubin, Chen Shengde, Bradley K Fritz. Current status and future trends of precision agricultural aviation technologies[J]. International Journal of Agricultural and Biological Engineering, 2017, 10(3): 1-17.

[2] 張東彥,蘭玉彬,陳立平,等. 中國農業航空施藥技術研究進展與展望[J]. 農業機械學報,2014,45(10):53-59.

Zhang Dongyan, Lan Yubin, Chen Liping, et al. Current status and future trends of agricultural aerial spraying technology in China[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(10): 53-59. (in Chinese with English abstract)

[3] 周志艷,臧英,羅錫文,等. 中國農業航空植保產業技術創新發展戰略[J]. 農業工程學報,2013,29(24):1-10.

Zhou Zhiyan, Zang Ying, Luo Xiwen, et al. Technology innovation development strategy on agricultural aviation industry for plant protection in China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(24): 1-10. (in Chinese with English abstract)

[4] 薛新宇,蘭玉彬. 美國農業航空技術現狀和發展趨勢分析[J]. 農業機械學報,2013,44(5):194-201.

Xue Xinyu, Lan Yubin. Agricultural aviation application in USA[J]. Transactions of the Chinese Society for Agricultural Machinery, 2013, 44(5): 194-201. (in Chinese with English abstract)

[5] 周志艷,明銳,臧禹,等. 中國農業航空發展現狀及對策建議[J]. 農業工程學報,2017,33(20):1-13.

Zhou Zhiyan, Ming Rui, Zang Yu, et al. Development status and countermeasures of agricultural aviation in China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(20): 1-13. (in Chinese with English abstract)

[6] He Xiongkui, Bonds J, Herbst A, et al. Recent development of unmanned aerial vehicle for plant protection in East Asia[J]. International Journal of Agricultural and Biological Engineering, 2017, 10(3): 18-30.

[7] 張豪,祁力鈞,吳亞壘,等. 基于Porous模型的多旋翼植保無人機下洗氣流分布研究[J]. 農業機械學報,2019,50(2):112-122.

Zhang Hao, Qi Lijun, Wu Yalei, et al. Spatio-temporal distribution of down–wash airflow for multi-rotor plant protection UAV based on porous model[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(2): 112-122. (in Chinese with English abstract)

[8] 李繼宇,蘭玉彬,施葉茵. 旋翼無人機氣流特征及大田施藥作業研究進展[J]. 農業工程學報,2013,34(12):104-118.

Li Jiyu, Lan Yubin, Shi Yeyin. Research progress on airflow characteristics and field pesticide application system of rotary-wing UAV[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(12): 104-118. (in Chinese with English abstract)

[9] 張宋超,薛新宇,秦維彩,等. N-3型農用無人直升機航空施藥飄移模擬與試驗[J]. 農業工程學報,2015,31(3):87-93.

Zhang Songchao, Xue Xinyu, Qin Weicai, et al. Simulation and experimental verification of aerial spraying drift on N-3 unmanned spraying helicopter[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(3): 87-93. (in Chinese with English abstract)

[10] 姚偉祥,蘭玉彬,王娟,等. AS350B3e直升機航空噴施霧滴飄移分布特性[J]. 農業工程學報,2017,33(22):75-83.

Yao Weixiang, Lan Yubin, Wang Juan, et al. Droplet drift law of aerial spraying of AS350B3e helicopter[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(22): 75-83. (in Chinese with English abstract)

[11] 王昌陵,宋堅利,何雄奎,等. 植保無人機飛行參數對施藥霧滴沉積分布特性的影響[J]. 農業工程學報,2017,33(23):109-116.

Wang Changling, Song Jianli, He Xiongkui, et al. Effect of flight parameters on distribution characteristics of pesticide spraying droplets deposition of plant–protection unmanned aerial vehicle[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(23): 109-116. (in Chinese with English abstract)

[12] Kang T G, Lee C S, Choi D K, et al. Development of aerial application system attachable to unmanned helicopter: basic spraying characteristics for aerial application system[J]. Journal of Biosystems Engineering, 2010; 35(4): 215-223.

[13] 陳盛德,蘭玉彬,李繼宇,等. 植保無人機航空噴施作業有效噴幅的評定與試驗[J]. 農業工程學報,2017,33(7):82-90.

Chen Shengde, Lan Yubin, Li Jiyu, et al. Evaluation and test of effective spraying width of aerial spraying on plant protection UAV[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(7): 82-90. (in Chinese with English abstract)

[14] 楊知倫,葛魯振,祁力鈞,等. 植保無人機旋翼下洗氣流對噴幅的影響研究[J]. 農業機械學報,2018,49(1):116-122.

Yang Zhilun, Ge Luzhen, Qi Lijun, et al. Influence of UAV rotor down-wash airflow on spray width[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(1): 116-122. (in Chinese with English abstract)

[15] Bruno S F, Freitas H, Gomes P H, et al. An adaptive approach for UAV-based pesticide spraying in dynamic environments[J]. Computers and Electronics in Agriculture, 2017, 138: 210-223.

[16] Bruno S. F, Pessin G, Filho G P R, et al. Fine-tuning of UAV control rules for spraying pesticides on crop fields: An approach for dynamic environments[J]. International Journal on Artificial Intelligence Tools, 2016, 25(1): 1-19.

[17] Ramasamy M, Lee T E, Leishman J G. Flow field of a rotating-wing micro air vehicle[J]. Journal of Aircraft, 2012, 44: 1236-1244.

[18] Wall B G, Richard H. Analysis methodology for 3C-PIV data of rotary wing vortices[J]. Experiments in Fluids, 2006, 40(5): 798-812.

[19] Yeo D, Shrestha E, Paley D A, et al. An empirical model of rotorcraft UAV downwash for disturbance localization and avoidance[C]//Kissimmee, Florida: AIAA Atmospheric Flight Mechanics Conference, 2015.

[20] 汪沛,胡煉,周志艷,等. 無人油動力直升機用于水稻制種輔助授粉的田間風場測量[J]. 農業工程學報,2013,29(3):54-61.

Wang Pei, Hu Lian, Zhou Zhiyan, et al. Wind field measurement for supplementary pollination in hybrid rice breeding using unmanned gasoline engine single-rotor helicopter[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(3): 54-61. (in Chinese with English abstract)

[21] 胡煉,周志艷,羅錫文,等. 無人直升機風場無線傳感器網絡測量系統設計與試驗[J]. 農業機械學報,2014,45(5):221-226.

Hu Lian, Zhou Zhiyan, Luo Xiwen, et al. Development and experiment of a wireless wind speed sensor network measurement system for unmanned helicopter[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(5): 221-226. (in Chinese with English abstract)

[22] 李繼宇,周志艷,蘭玉彬,等. 旋翼式無人機授粉作業冠層風場分布規律[J]. 農業工程學報,2015,31(3):77-86.

Li Jiyu, Zhou Zhiyan, Lan Yubin, et al. Distribution of. canopy wind field produced by rotor unmanned aerial vehicle pollination operation[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(3): 77-86. (in Chinese with English abstract)

[23] 陳盛德,蘭玉彬,Bradley K F,等. 多旋翼無人機旋翼下方風場對航空噴施霧滴沉積的影響[J]. 農業機械學報,2017,48(8):105-113.

Chen Shengde, Lan Yubin, Bradley K F, et al. Effect of wind field below rotor on distribution of aerial spraying droplet deposition by using multi–rotor UAV[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(8): 105-113. (in Chinese with English abstract)

[24] Yang Fengbo, Xue Xinyu, Zhang Ling, et al. Numerical simulation and experimental verification on downwash air flow of six-rotor agricultural unmanned aerial vehicle in hover[J]. International Journal of Agricultural and Biological Engineering, 2017, 10(4): 41-53.

[25] 楊風波,薛新宇,蔡晨,等. 多旋翼植保無人飛機懸停下洗氣流對霧滴運動規律的影響[J]. 農業工程學報,2018,34(2):64-73.

Yang Fengbo, Xue Xinyu, Cai Chen, et al. Effect of downwash airflow in hover on droplet motion law for multi-rotor unmanned plant protection machine[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(2): 64-73. (in Chinese with English abstract)

[26] 文晟,韓杰,蘭玉彬,等. 單旋翼植保無人機翼尖渦流對霧滴飄移的影響[J]. 農業機械學報,2018,49(8):127-137,160.

Wen Sheng, Han Jie, Lan Yubin, et al. Influence of wing tip vortex on drift of single rotor plant protection unmanned aerial vehicle[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(8): 127-137, 160. (in Chinese with English abstract)

[27] Endalew A M, Hertog M, Delele M A, et al. CFD modelling and wind tunnel validation of airflow through plant canopies using 3D canopy architecture[J]. International Journal of Heat and Fluid Flow, 2009, 30(2): 356-368.

[28] Endalew A M, Debaer C, Rutten N, et al. A new integrated CFD modelling approach towards air–assisted orchard spraying-Part I. Model development and effect of wind speed and direction on sprayer airflow[J]. Computers and Electronics in Agriculture, 2010, 71(2): 128-136.

[29] Endalew A M, Debaer C, Rutten N, et al. A new integrated CFD modelling approach towards air–assisted orchard spraying-Part II: Validation for different sprayer types[J]. Computers and Electronics in Agriculture, 2010, 71(2): 137-147.

[30] Endalew A M, Debaer C, Rutten N, et al. Modelling pesticide flow and deposition from air–assisted orchard spraying in orchards: A new integrated CFD approach[J]. Agricultural and Forest Meteorology, 2010, 150(10): 1383-1392.

[31] Hong S W, Zhao Lingying, Zhu Heping. CFD simulation of airflow inside tree canopies discharged from air-assisted sprayers[J]. Computers and Electronics in Agriculture, 2018, 149: 121-132.

[32] Hong S W, Zhao Lingying, Zhu Heping. CFD simulation of pesticide spray from air-assisted sprayers in an apple orchard: Tree deposition and off-target losses[J]. Atmospheric Environment, 2018, 175: 109-119.

[33] Duga A T, Dekeyser D, Ruysen K, et al. Numerical analysis of the effects of wind and sprayer type on spray distribution in different orchard training systems[J]. Boundary Layer Meteorology, 2015, 157(3): 1-19.

[34] Duga A T, Delele M A, Ruysen K, et al. Development and validation of a 3D CFD model of drift and its application to air-assisted orchard sprayers[J]. Biosystems Engineering, 2017, 154(2): 62-75.

[35] Salcedo R, Vallet A, Granell R, et al. Eulerian-Lagrangian model of the behaviour of droplets produced by an air–assisted sprayer in a citrus orchard[J]. Biosystems Engineering, 2017, 154(2): 76-91.

[36] Salcedo R, Granell R, Palau G, et al. Design and validation of a 2D CFD model of the airflow produced by an airblast sprayer during pesticide treatments of citrus[J]. Computers and Electronics in Agriculture, 2015, 116: 150-161.

[37] Hong S W, Zhao Lingying, Zhu Heping. SAAS, a computer program for estimating pesticide spray efficiency and drift of air–assisted pesticide applications[J]. Computers and Electronics in Agriculture, 2018, 155: 58-68.

[38] 李鵬飛,許敏儀,王飛飛. 精通CFD工程仿真與案例實戰:FLUENT GAMBIT ICEM CFD Tecplot[M]. 北京:人民郵電出版社,2011:186-192.

[39] 肖中云. 旋翼流場數值模擬方法研究[D]. 綿陽:中國空氣動力研究與發展中心,2007.

Xiao Zhongyun. Investigation of Computational Modeling Techniques for Rotor Flowfields[D]. Mianyang: China Aerodynamics Research and Development Center, 2007. (in Chinese with English abstract)

[40] 王福軍. 計算流體動力學分析:CFD軟件原理與應用[M]. 北京:清華大學出版社,2004.

[41] 程湞湞,祁力鈞,吳亞壘,等. 矮化密植果園搖擺變量噴霧機參數響應面法優化[J]. 農業機械學報,2017,48(增刊):22-29.

Cheng Zhenzhen, Qi Lijun, Wu Yalei, et al. Parameter optimization on swing variable sprayer of orchard based on RSM[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(Supp.): 22-29. (in Chinese with English abstract)

Distribution characteristics of rotor downwash airflow field under spraying on orchard using unmanned aerial vehicle

Zhang Hao, Qi Lijun※, Wu Yalei, Cheng Zhenzhen, Liu Wanwan, Elizabeth Musiu, Xiao Yu, Yang Zepeng

(C,,100083,)

when the plant protection unmanned aerial vehicle (UAV) is used to spray pesticides on orchard, the distribution of rotor downwash airflow filed has significant influence on the spatial movement of the droplet and the adhesion and penetration of the droplet inside the canopy. Based on computational fluid dynamics (CFD) method, combined with RNG-turbulence model, porous model and sliding mesh technology, the rotor downwash airflow field of a six–rotor plant protection UAV in hover when spraying on orchard was simulated. The simulation was done in the constructed virtual orchard. The characteristics of the airflow field were analyzed in different hovering heights of the UAV, fruit growth stages and natural wind speeds. Verification experiments were carried out through measuring the downwash airflow velocity at marked points. The results showed that: 1) it no longer met the conditions for spraying of the plant protection UAV in hover when the natural wind speed was greater than 3 m/s due to the downwash airflow under the rotor submerged in the natural wind of the environment. 2) Natural wind destroyed the central symmetry of downwash airflow of the rotor, and airflow diffusion appeared along the downwind direction. With the increase of natural wind speed and hovering height, the backward lift distance increased. When the hovering height was 3 m, and the natural wind speed was 1 and 2 m/s, the trailing distance of the rotor under the airflow reached 1 and 2 m respectively; when the natural wind speed was 3 m/s, the trailing distance had been more than 2 m. Under the condition of natural wind speed of 2 m/s, and the hovering height was 3 and 3.5 m, the trailing distances of the rotor under the airflow were not much different, both were both 2 m, but the former was in contact with the target canopy layer, and the latter was not in contact; When the hovering height was 4 m, the trailing distance of airflow had exceeded 2 m. 3) Compared with the state of no natural wind, the velocity distribution at the nozzle was mainly affected by natural wind, but the effect of fruit growth stages was not significant. In addition, the vertical-direction airflow played a leading role in the target transport of the droplets. The spray head should be installed near 0.2 m directly below the rotor to make the droplets to have large–direction velocity. 4) After the hovering position of the UAV was adjusted in the upwind direction, the average velocity of the upper, middle and lower airflows in the canopy increased from 1.36, 0.80 m/s, and 0.81 to 3.04 m/s, 2.37 and 1.63 m/s, respectively. The coefficient of variation of velocity distribution in the upper and lower layers decreased from 74.26% and 35.80% to 45.39% and 22.70%, respectively, and the middle layer increased slightly, which was beneficial to achieve target spray. The experimental results showed that there was a good consistency between the experimental and simulated values of downwash airflow velocity at marked points. In conclusion, this paper should provide further reference for the development of the adaptive control technology of plant protection UAV hovering target spraying in a dynamic environment.

unmanned aerial vehicle; pesticide; models; downwash airflow; fruit tree canopy; porous model

張 豪,祁力鈞,吳亞壘,程湞湞,劉婠婠,Elizabeth Musiu,肖 雨,楊澤鵬. 無人機果樹施藥旋翼下洗氣流場分布特征研究[J]. 農業工程學報,2019,35(18):44-54.doi:10.11975/j.issn.1002-6819.2019.18.006 http://www.tcsae.org

Zhang Hao, Qi Lijun, Wu Yalei, Cheng Zhenzhen, Liu Wanwan, Elizabeth Musiu, Xiao Yu, Yang Zepeng. Distribution characteristics of rotor downwash airflow field under spraying on orchard using unmanned aerial vehicle[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(18): 44-54. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2019.18.006 http://www.tcsae.org

2019-03-19

2019-08-09

科技部國家重點研發計劃項目“現代果園智能化精細生產管理技術裝備研發”(2017YFD0701400);科技部國家重點研發計劃項目“地面與航空高工效施藥技術及智能化裝備(2016YFD0200700)”。

張 豪,博士生,主要從事植保機械研究。Email: zhanghao08@cau.edu.cn

祁力鈞,博士,教授,主要從事植保機械研究。Email:qilijun@cau.edu.cn

10.11975/j.issn.1002-6819.2019.18.006

S49; S25

A

1002-6819(2019)-18-0044-11

猜你喜歡
風速
邯鄲市近46年風向風速特征分析
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
基于時間相關性的風速威布爾分布優化方法
陜西黃土高原地區日極大風速的統計推算方法
陜西氣象(2020年2期)2020-06-08 00:54:38
基于GARCH的短時風速預測方法
快速評估風電場50年一遇最大風速的算法
風能(2016年11期)2016-03-04 05:24:00
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
GE在中國發布2.3-116低風速智能風機
考慮風速分布與日非平穩性的風速數據預處理方法研究
主站蜘蛛池模板: 国产美女久久久久不卡| 成人福利在线观看| 亚洲国产清纯| 亚洲中文字幕无码爆乳| 不卡无码h在线观看| 人禽伦免费交视频网页播放| 国产JIZzJIzz视频全部免费| 91激情视频| 国产午夜福利在线小视频| 国产a v无码专区亚洲av| 亚洲视频在线青青| a级毛片免费在线观看| 久久国产免费观看| 日韩久草视频| 2021无码专区人妻系列日韩| 欧美日韩免费在线视频| 免费av一区二区三区在线| 国内熟女少妇一线天| av午夜福利一片免费看| 少妇露出福利视频| 免费人欧美成又黄又爽的视频| 中文字幕 欧美日韩| 狼友视频一区二区三区| 色135综合网| 欧美成人午夜视频免看| 91区国产福利在线观看午夜 | 丁香五月亚洲综合在线| 日韩av无码DVD| 在线观看亚洲精品福利片| 高清久久精品亚洲日韩Av| 日本在线视频免费| 成人在线综合| 91在线视频福利| 在线观看免费人成视频色快速| 日韩精品成人在线| 亚洲综合第一区| 色哟哟国产精品一区二区| 伊人久久精品无码麻豆精品| 日韩第九页| 亚洲第一区欧美国产综合| 国产在线专区| 欧美成人日韩| 久久伊伊香蕉综合精品| 国产XXXX做受性欧美88| 99热这里只有精品2| 欧美中文字幕在线视频| 亚洲精品中文字幕无乱码| 免费国产一级 片内射老| 米奇精品一区二区三区| 狠狠综合久久| 日本91在线| 午夜在线不卡| 久久精品人人做人人爽电影蜜月| 无码专区第一页| 精品无码一区二区在线观看| 992tv国产人成在线观看| 少妇露出福利视频| 国产成人高清精品免费5388| 91视频国产高清| 永久免费无码日韩视频| 1024你懂的国产精品| 精品一区二区久久久久网站| 亚洲日本在线免费观看| 精品无码日韩国产不卡av| 五月激激激综合网色播免费| 日韩不卡高清视频| 美女毛片在线| 97青草最新免费精品视频| 香蕉综合在线视频91| 国产无码网站在线观看| 欧美a级在线| 亚洲欧美在线综合图区| 熟妇人妻无乱码中文字幕真矢织江 | 一级毛片免费高清视频| 青青青国产视频手机| 日本午夜在线视频| 亚洲日韩精品伊甸| 欧美天堂在线| 九色在线观看视频| 国产成人夜色91| 国产三级成人| 97人人做人人爽香蕉精品|