采用雙向流固耦合方法構(gòu)建輔助氣流作用下棉花葉片變形模型

劉興華，苑 進(jìn)，李 揚，索金政，劉雪美

采用雙向流固耦合方法構(gòu)建輔助氣流作用下棉花葉片變形模型

劉興華1,2，苑 進(jìn)1,3，李 揚1,3，索金政1，劉雪美1,2※

（1. 山東農(nóng)業(yè)大學(xué)機(jī)械與電子工程學(xué)院，泰安 271018；2. 山東省農(nóng)業(yè)裝備智能化工程實驗室，泰安 271018； 3. 山東省園藝機(jī)械與裝備重點實驗室，泰安 271018）

氣流輔助噴霧中，輔助氣流能夠使冠層葉片變形、拓寬霧滴輸運通道，有助于提高霧滴在作物冠層內(nèi)沉積均勻性。為進(jìn)一步研究輔助氣流作用下葉片變形機(jī)理，以棉花葉片為研究對象，首先，在測定葉片力學(xué)參數(shù)的基礎(chǔ)上，建立了棉花葉片雙向流固耦合模型，并利用高速攝像技術(shù)驗證了模型的有效性；進(jìn)一步，借助耦合模型，以輔助氣流風(fēng)速、葉傾角、葉柄及葉面彈性模量為試驗因素，以葉柄和葉面撓度值為評價指標(biāo)，開展正交試驗分析，得出試驗因素影響葉片變形的顯著性順序依次為：氣流速度、彈性模量、葉傾角；最后，基于非線性最小二乘法構(gòu)建棉花葉片變形的參數(shù)辨識模型，所建的葉柄撓曲線方程的平均絕對百分比誤差（Mean Absolute Percentage Error, MAPE）為5.13%，葉面主葉脈撓曲線方程的MAPE為10.43%，證明了葉片變形數(shù)學(xué)模型的有效性。依據(jù)建立的葉片變形模型，量化分析了不同風(fēng)速與初始葉傾角參數(shù)組合下葉片迎風(fēng)面積的動態(tài)變化過程，揭示了速度不恰當(dāng)?shù)妮o助氣流作用下冠層郁閉度增加的原因，明確了初始葉傾角為20°和30°時，棉花葉片的輔助氣流末速度應(yīng)分別不小于2.4和 2.9 m/s。本研究為理解輔助氣流噴霧中作物冠層郁閉度的動態(tài)變化、合理選取施藥工作參數(shù)提供參考。

輔助氣流；模型；葉片變形；流固耦合；參數(shù)辨識

0 引 言

在生長發(fā)育中后期，棉花、玉米、大豆等作物具有枝葉茂密、葉面積指數(shù)大、冠層郁閉度高等特點，導(dǎo)致植保噴霧作業(yè)面臨霧滴輸運通道閉塞、冠層霧滴沉積不均勻難題。針對上述問題，目前多采用氣流輔助噴霧[1-2]和吊桿噴霧[3-4]兩種方式突破枝葉屏蔽，拓寬霧滴輸運通道，以在冠層內(nèi)上中下層實現(xiàn)更好的沉積均勻性。氣流輔助噴霧方式提高藥液霧滴在冠層內(nèi)的沉積均勻性是通過以下兩個方面來實現(xiàn)：一是對作物枝葉施加風(fēng)力載荷，引起葉片變形，拓寬霧滴輸運通道；二是提高藥液霧滴運動速度，加強霧滴在冠層內(nèi)的運移穿透能力[5]。因此，針對拓寬霧滴輸運通道，分析輔助氣流引起的葉片變形規(guī)律是理解作物冠層內(nèi)部霧滴運移沉積規(guī)律和優(yōu)化調(diào)控氣流輔助噴霧作業(yè)參數(shù)的重要前提。

針對高郁閉度作物氣流輔助噴霧，國內(nèi)外學(xué)者在作物冠層與輔助氣流場、霧滴沉積分布的相互影響方面開展了大量研究。魏新華等[4]針對棉花冠層施藥時枝葉遮擋問題，提出冠內(nèi)冠上組合風(fēng)送式噴桿結(jié)構(gòu)，有效提高了冠層霧滴沉積均勻性；李龍龍等[6]根據(jù)LiDAR掃描獲取的果樹形態(tài)特征，采用仿形變量風(fēng)送式噴霧方式，提高了藥液利用率和沉積均勻性。上述機(jī)型設(shè)計與試驗結(jié)果表明，作物冠層特征影響霧滴沉積分布。許童羽等[7]分析無人機(jī)霧滴在水稻冠層沉積分布，發(fā)現(xiàn)霧滴沉積分布與下洗氣流關(guān)聯(lián)度高，并受到水稻冠層枝葉的影響；孫誠達(dá)等[8]研究發(fā)現(xiàn)風(fēng)送噴霧模式下霧滴穿透比例隨葉密度、取樣深度的增加而減少；Duga等[9]開展對蘋果、梨樹的田間噴霧試驗，研究發(fā)現(xiàn)冠層孔隙率和體積顯著影響輔助氣流噴霧過程中霧滴的靶標(biāo)沉積狀況，并且霧滴沉積與輔助氣流存在強關(guān)聯(lián)性?？紤]到輔助氣流對霧滴沉積分布的重要影響，Endalew等[10-11]和Hong等[12]分別建立果樹冠層枝葉的等效多孔介質(zhì)模型，探究了冠層枝葉分布和密度對氣流場分布和風(fēng)速衰減的重要影響；Li等[13]分析了無人機(jī)旋翼風(fēng)場在水稻冠層內(nèi)部的分布，發(fā)現(xiàn)冠層高度越低，風(fēng)場等效面積衰減率越大。上述研究側(cè)重于作物枝葉對冠層風(fēng)場分布和風(fēng)速衰減的影響，并證實了作物枝葉對風(fēng)送霧滴沉積分布的影響，證明了開展輔助氣流與冠層枝葉交互研究的必要性。但是，目前輔助氣流作用下作物葉片變形的研究較少。劉希鳳等[14]研究了風(fēng)載激勵下王棕、散尾葵等6種葉片的最大變形量，揭示了氣流與葉面夾角對變形影響顯著。目前，尚未發(fā)現(xiàn)針對棉花葉片的風(fēng)力載荷變形研究。

本文以棉花葉片為研究對象，建立風(fēng)力載荷下棉花葉片彎曲變形的雙向流固耦合模型，以葉片彈性模量、葉傾角、輔助氣流流速三因素開展正交試驗，量化分析輔助氣流下棉花葉片的變形行為，建立棉花葉柄、葉面彎曲變形的數(shù)學(xué)模型，以葉片迎風(fēng)面積變化率為指標(biāo)，分析了不同參數(shù)組合下葉片迎風(fēng)面積的變化規(guī)律。

1 棉花葉片仿真模型構(gòu)建與驗證

1.1 棉花葉片參數(shù)測定

1.1.1 葉片尺寸

棉花葉片取自山東省泰安市山東農(nóng)業(yè)大學(xué)農(nóng)學(xué)實驗站（117.168177°E，36.165072°N），棉花品種為魯棉研37號。由于花鈴期枝葉茂密、冠層郁閉度高，是氣流輔助噴霧的主要作業(yè)時期，因此在花鈴期前期、中期、后期3個批次采集試樣。為考慮不同部位葉片的葉齡差異，將棉花冠層沿垂直方向分為上、中、下3層，在每層隨機(jī)采集3片葉面和3根葉柄試樣，兩種試樣各有27份。

為后續(xù)仿真建模和力學(xué)參數(shù)測定的需要，用直尺測量葉面的長度與寬度、葉柄長度，用游標(biāo)卡尺（150 mm，0.01 mm）測量葉柄直徑（始端、中部、末端測量值取平均），尺寸數(shù)據(jù)見表1。同時，為避免試樣失水影響參數(shù)測定準(zhǔn)確性，對葉面和葉柄試樣編號并裝入自封袋內(nèi)，迅速轉(zhuǎn)移至試驗室內(nèi)進(jìn)行力學(xué)參數(shù)測定。

表1 不同棉花冠層的葉面和葉柄試樣尺寸參數(shù)

1.1.2 葉片彈性模量

借助物性分析儀（TA-XT.Plus 英國Stable Micro Systems公司），采用三點彎曲法測定棉花葉片葉面和葉柄部分的彈性模量[15]。將待測試樣放置在支撐立壁上，探頭垂直向下運動，探頭下壓過程中采集試樣中點撓度值和載荷增量值。物性分析儀參數(shù)設(shè)定：探頭：P5（曲率半徑5 mm）；測試模式：壓縮；目標(biāo)模式：位移；數(shù)據(jù)采集頻率：400 Hz；支撐立壁跨度：20 mm；測前速度和測中速度設(shè)為1 mm/s，測后速度設(shè)為10 mm/s。在葉柄始端和中部截取2段長度為25 mm的試樣，借助物性分析儀獲取其應(yīng)力-位移曲線，利用公式（1）求得葉柄彈性模量p；沿葉面主葉脈方向，在葉面根部、中部、尾部截取3塊長度為25 mm、寬度為10 mm的試樣，借助物性分析儀獲取其應(yīng)力-位移曲線，利用公式（2）求得葉面彈性模量l[16]，測量過程如圖1所示。

式中Δ為載荷增量，N；Δ為中心撓度增量，mm；為水平支架跨距，mm；為慣性矩，mm4；為泊松比，葉面泊松比0.32，葉柄泊松比0.34[17]；

1.1.3 葉片含水率

試樣含水率測定采用105 ℃烘干法，參照GB/T1931-2009標(biāo)準(zhǔn)的規(guī)定開展試驗。將測定完力學(xué)參數(shù)的新鮮棉花試樣快速轉(zhuǎn)移至精密電子天平（FA1004 寧波市鄞州華豐電子儀器廠）內(nèi)進(jìn)行稱量，測量精度為0.000 1 g，然后將試樣轉(zhuǎn)至真空干燥箱（DZF-6050上海博訊實業(yè)有限公司）內(nèi)，在（105±2）℃溫度下烘干8 h，取出放干燥器中冷卻后立即稱量，重復(fù)烘干3 h，冷卻后再次稱量，連續(xù)兩次干質(zhì)量之差不大于0.001 g，停止烘干，獲得試樣恒質(zhì)量，計算試樣含水率。

1.1.4 參數(shù)測定結(jié)果與分析

統(tǒng)計葉面和葉柄材料的彈性模量，剔除異常數(shù)據(jù)樣本，獲得葉面彈性模量中位數(shù)為46.5 MPa，95%置信區(qū)間為[28.5, 64.5] MPa；葉柄彈性模量中位數(shù)為244.6 MPa，95%置信區(qū)間為[215.5, 273.69] MPa。葉柄含水率中位數(shù)為87.2%，95%置信區(qū)間為[82.5%, 91.9%]，葉面含水率中位數(shù)為80.7%，95%置信區(qū)間為[72.3%, 89.1%]。根據(jù)獲取的試驗數(shù)據(jù)，對25份葉面材料的彈性模量相關(guān)性分析表明，Pearson相關(guān)系數(shù)為=0.171 (>0.05)，相關(guān)性不顯著；對24份材料的彈性模量與含水率相關(guān)性分析表明，Pearson相關(guān)系數(shù)為=?0.174 (>0.05)，相關(guān)性不顯著。其原因可能在于，棉花葉面、葉柄含水率普遍較高而纖維含量均比較低；試樣彈性模量主要受纖維含量的影響，而受含水率的影響較小。因此，后續(xù)正交試驗不將含水率納入試驗因素。

1.2 葉片變形流固耦合模型構(gòu)建

流固耦合現(xiàn)象是由流體域與固體域相互影響而產(chǎn)生的，其分析求解包括流場的非定常計算和結(jié)構(gòu)體的瞬態(tài)動力學(xué)分析兩部分構(gòu)成。本文關(guān)注的輔助氣流作用下棉花葉片的變形行為分析是典型的流固耦合現(xiàn)象。在本文的分析過程中，流體計算采用Fluent 模塊，固體計算采用 Transient Structural模塊，利用ANSYS的system coupling模塊實現(xiàn)流場和結(jié)構(gòu)體的數(shù)據(jù)交互，并采用順序求解方式，依次求解流體域與固體域。

1.2.1 控制方程

考慮輔助氣流作用下棉花葉片的變形過程，輔助氣流的流動是三維不可壓縮流體的低速流動，其雷諾平均N-S方程為[18-19]

輔助氣流載荷作用下，棉花葉片的瞬態(tài)動力學(xué)分析式為[20]

1.2.2 計算模型

棉花葉片外形特征：闊卵形，長寬近相等，多見3淺裂，少有5裂[21]。建立3裂棉花葉片三維模型，如圖2所示。依據(jù)氣流輔助噴霧主要作業(yè)期（花鈴期）的魯棉研37號棉花的掌狀3裂葉片外形和采集的葉片試樣的統(tǒng)計均值，確定模型尺寸參數(shù)：葉柄長110 mm，葉面寬120 mm，葉面長110 mm。由于本文主要關(guān)注棉花葉片在輔助氣流速度、葉傾角、葉片力學(xué)參數(shù)等因素影響下葉片變形及迎風(fēng)面積變化規(guī)律等內(nèi)容，同時器官生長的同伸關(guān)系使得葉面葉柄尺寸比例基本保持不變，因此分析過程不對葉片尺寸參數(shù)進(jìn)行調(diào)整。流體計算域邊界：4個側(cè)立面距離葉片邊緣均為150 mm，設(shè)置為自由出口；上平面距離葉片200 mm，設(shè)置為速度入口；下平面距離葉片300 mm，設(shè)置為自由出口。

棉花葉片由葉柄和葉面兩部分構(gòu)成，采用分塊網(wǎng)格劃分方式生成六面體網(wǎng)格；流體域采用ANSYS Meshing的CFD網(wǎng)格劃分方式，網(wǎng)格增長率為1.15。流體域求解過程中，使用動網(wǎng)格技術(shù)模擬棉花葉片在風(fēng)力載荷下的動態(tài)變形過程，網(wǎng)格光順采用Diffusion方式，葉片上下迎風(fēng)面設(shè)定為耦合交界面，壓力速度耦合采用Coupled算法；固體域求解過程中，求解器控制開啟Large Deflection功能；在System Coupling模塊中設(shè)置固體域和流體域的時間步長為0.001 s，最小迭代步數(shù)為1步，最大迭代步數(shù)為5步，仿真計算時間為0.5 s。在輔助氣流下棉花葉片變形高速圖像試驗中，當(dāng)風(fēng)速處于適宜棉花施藥農(nóng)藝要求的范圍時（1.5～3 m/s）[22]，葉片尖端的變形量小，而葉片沿葉柄和葉面主葉脈彎曲變形明顯，并且葉柄約葉面結(jié)合部的變形幅度最大。因此，沿葉片中軸線，選擇葉柄中點、葉柄末端、葉面形心、葉尖4個監(jiān)測點，用于記錄分析葉片不同位置的變形情況。

1.3 耦合模型有效性驗證

采用高速攝像機(jī)（Phantom VEO410L，Vision Research），搭配微距鏡頭（Nikon AF-S VR，105 mm f/2.8G），使用數(shù)碼燈（LED-200B，5600 K，20 000 Lm）對棉花葉片進(jìn)行補光，保持?jǐn)z像機(jī)光軸垂直于參照平面和運動平面，拍攝和記錄棉花葉片的動態(tài)變形過程。借助TEMA圖像后處理軟件分析葉片不同部位的變形量，對流固耦合模型的有效性進(jìn)行驗證。

氣流輔助噴霧施藥過程中，冠層內(nèi)的棉花葉片以不同的迎風(fēng)角度承受風(fēng)力載荷，發(fā)生彈性變形。因此，以氣流速度、葉傾角（即葉片腹面法線與天頂軸的夾角）、葉柄彈性模量p和葉面彈性模量l為變量，以葉柄末端、葉柄中點、葉面形心、葉尖撓度為指標(biāo)，對耦合模型進(jìn)行驗證。依據(jù)末速度原則[22]，風(fēng)送噴霧機(jī)的風(fēng)送氣流末速度為2 m/s[4]，考慮不同部位葉片的風(fēng)速差異，風(fēng)速范圍定為1.5～3.0 m/s。依據(jù)棉花葉傾角分布概率密度函數(shù)[23]，盛花期葉片的葉傾角處于0°～30°的概率為83.3%，其余角度所占概率之和為16.7%，本文選取的葉傾角范圍為0°～30°。由于器官發(fā)育的同伸關(guān)系，葉柄與葉面的彈性模量隨葉齡具有相同的變化規(guī)律[5]。因此，葉柄和葉面彈性模量視作1個參數(shù)組合。

試驗研究中發(fā)現(xiàn)，由于輔助氣流吹動時葉片快速失水、萎蔫影響變形準(zhǔn)確性，輔助氣流作用下葉片變形試驗的開展難度大。因此，為充分驗證耦合模型在風(fēng)速取值域內(nèi)、彈性模量分布范圍內(nèi)的有效性，同時減小驗證試驗工作量，測定和遴選了4種不同力學(xué)參數(shù)的葉片，與不同的風(fēng)速、葉傾角組合，開展驗證試驗，具體參數(shù)組合見表2。

試驗前，在葉片夾持工作臺前設(shè)置四象限標(biāo)記點，用于建立參照平面；在葉柄中點、葉柄末端、葉面形心及葉尖設(shè)置四象限標(biāo)記點作為運動追蹤點，用于記錄葉片變形量，試驗布置如圖3a所示。試驗中，對于每一葉片，重復(fù)拍攝3段高速動態(tài)圖像，用于后續(xù)的變形數(shù)據(jù)提取。試驗后，利用TEMA軟件進(jìn)行圖像后處理。由于輔助氣流作用下棉花葉片處于動態(tài)變形過程，因此采用多重平行面方法在參照平面進(jìn)行尺寸標(biāo)定，根據(jù)參照平面與葉片運動平面的距離進(jìn)行縱深補償，利用Tracking功能獲得葉柄和葉面標(biāo)記點的變形數(shù)據(jù)，取3次試驗的均值。

表2 用于模型驗證的3種試驗參數(shù)組合

根據(jù)表2試驗參數(shù)，設(shè)置葉片雙向流固耦合模型的仿真參數(shù)，在葉柄末端和葉面形心設(shè)置監(jiān)測點，獲取4種參數(shù)組合下的葉柄和葉面變形數(shù)據(jù)，如圖3b所示。對比耦合仿真與高速攝像結(jié)果，如表3所示，4種參數(shù)組合下相對誤差最大值為9.85%；同時，分析耦合仿真中棉花葉片的動態(tài)變形特點，發(fā)現(xiàn)葉柄和葉面結(jié)合部發(fā)生較大幅度應(yīng)變，進(jìn)而導(dǎo)致葉面沿風(fēng)速方向彎轉(zhuǎn)，與實際試驗中棉花葉片的變形特點一致。上述結(jié)果表明本文所建立的棉花葉片雙向流固耦合模型能夠有效仿真輔助氣流作用下棉花葉片的變形過程，可以為變形規(guī)律分析與數(shù)學(xué)模型構(gòu)建提供支持。

表3 基于耦合模型的葉片撓度值仿真誤差統(tǒng)計

2 棉花葉片彎曲變形數(shù)學(xué)模型構(gòu)建

利用搭建的棉花葉片流固耦合模型，考慮輔助氣流風(fēng)速、葉傾角、葉柄及葉面力學(xué)參數(shù)等因素，對輔助氣流作用下棉花葉片變形開展正交試驗分析，依據(jù)材料力學(xué)原理并結(jié)合正交試驗因素顯著性排序，確定葉片變形數(shù)學(xué)模型框架，使用葉片變形仿真數(shù)據(jù)，通過最小二乘參數(shù)辨識方法識別模型參數(shù)，建立輔助氣流作用下棉花葉片變形數(shù)學(xué)模型。

2.1 正交試驗

2.1.1 試驗方案

在正交試驗設(shè)計中，以氣流速度、葉傾角、葉柄與葉面彈性模量（p，l）為試驗因素，各因素在參數(shù)取值域內(nèi)均勻選取4個水平，確定正交試驗水平及方案如表4。選用L16(45)正交試驗表，以葉柄中點撓度、葉柄末端撓度、葉面形心撓度和葉尖撓度為試驗指標(biāo)，開展正交試驗。

表4 正交試驗因素水平表

2.1.2 試驗結(jié)果與分析

正交試驗結(jié)果如表5。

表5 正交試驗方案及結(jié)果

對葉柄中點、末端和葉面形心、葉尖的變形量開展極差分析。由表6知，4個監(jiān)測點的撓度變化過程中，風(fēng)速、葉傾角和葉片力學(xué)參數(shù)（p，l）3個因素對4個撓度試驗指標(biāo)的影響從大到小均為：、（p，l）、。排序結(jié)果表明輔助氣流噴霧作業(yè)中風(fēng)速對葉片變形影響最大。

表6 葉片變形影響因素極差分析

2.2 葉片變形數(shù)學(xué)模型

2.2.1 葉片受力變形基本方程

棉花葉片由葉柄和葉面兩部分構(gòu)成，由于兩部分的生物力學(xué)參數(shù)差異較大，同一個撓曲線模型、同一組模型參數(shù)不易取得理想的建模精度?？紤]到葉面的變形和葉面積變化是在葉柄彎曲變形基礎(chǔ)上由葉面的彎曲變形引起的，本文對耦合模型提供的葉面變形數(shù)據(jù)分成葉柄變形分量和葉面變形分量兩部分，先后建立了葉柄和葉面主葉脈的撓曲線方程，從而獲得輔助氣流作用下葉片變形的數(shù)學(xué)模型。葉柄呈長軸狀，末端承受葉面帶來的輔助氣流風(fēng)力載荷，可視為懸臂梁結(jié)構(gòu)[14]。構(gòu)成葉面的葉肉和葉脈兩部分的剛?cè)狁詈媳戎递^高，主葉脈是主要的機(jī)械支撐結(jié)構(gòu)；氣流作用下葉片彎曲過程，可視為葉面均布載荷下主葉脈懸臂梁的彎曲變形[24-25]。因此，本文對兩部分均以材料力學(xué)中懸臂梁受力變形的撓曲線方程為框架，利用正交仿真獲取的多監(jiān)測點變形量數(shù)據(jù)，確定最優(yōu)辨識參數(shù)，建立輔助氣流作用下棉花葉片變形數(shù)學(xué)模型。

懸臂梁在集中載荷作用下發(fā)生彈性小變形，其撓曲線方程為[26]

式中()為受力彎矩，N·m；為集中載荷作用力，N；為懸臂梁長度，mm；為撓度值分析點的橫坐標(biāo)，mm；()為懸臂梁處的撓度值，mm；為材料彈性模量，MPa；為懸臂梁慣性矩，mm4。

輔助氣流作用下，棉花葉面是主要的風(fēng)力載荷承受部位；初始瞬間，葉面沒有發(fā)生彈性變形，屬于剛性體，其風(fēng)力載荷為[27]

式中為空氣密度，1.293 kg/m3；為輔助氣流速度，m/s；為棉花葉片面積，mm2；為葉片葉傾角，(°)。

2.2.2 模型參數(shù)辨識算法

棉花葉片風(fēng)載變形建?？蓺w類為已知基的非線性辨識問題，即可以將非線性項寫成待辨識參數(shù)乘以已知基的形式，利用過參數(shù)化方法并結(jié)合線性系統(tǒng)的辨識方法，如最小二乘算法、隨機(jī)梯度算法、極大似然估計算法等來辨識非線性系統(tǒng)的參數(shù)[28]。柔性棉花葉柄的受力變形，可視為在懸臂梁末端施加變化的風(fēng)力集中載荷時懸臂梁的非線性大變形。聯(lián)立式（5）～（7），構(gòu)造撓曲線函數(shù)框架；考慮葉片變形的非線性特點，設(shè)計風(fēng)速二次多項式（包含1、2、3）代替原風(fēng)速平方項，設(shè)計考慮葉片風(fēng)載變形的葉傾角指數(shù)式（包含4）代替原葉傾角常數(shù)項，設(shè)計葉柄長度指數(shù)式（包含5）代替原葉柄長度一次多項式變量[29]。葉柄撓曲線方程為

棉花葉面的受力變形，可視為葉面均布載荷下主葉脈懸臂梁的彎曲變形過程[25]。因此，在考慮風(fēng)速、風(fēng)力載荷動態(tài)變化和非線性大變形等因素的同時，葉面主葉脈撓曲線方程可表示為

基于信賴域的最小二乘法是辨識非線性數(shù)學(xué)模型中待定參數(shù)的有效方法，信賴域迭代法能夠更好地解決局部極值問題[30]。棉花葉柄和葉面撓曲線方程中，參數(shù)向量和的辨識可借助上述方法進(jìn)行，其最小二乘法的目標(biāo)函數(shù)如式（10）所示；利用正交試驗獲取的棉花葉柄和葉面撓度數(shù)據(jù)，其信任域算法迭代計算如式（11）所示。

式中*為根據(jù)試驗數(shù)據(jù)利用最小二乘法獲取的葉柄或葉面撓曲線最優(yōu)辨識參數(shù)向量。

其中

式中()為關(guān)于撓曲線方程()的Jacobian矩陣；為試驗數(shù)據(jù)組數(shù)；為辨識參數(shù)向量長度，即待辨識參數(shù)的數(shù)量；為迭代次數(shù)。

2.2.3模型精度分析與驗證

1）模型精度分析

利用表5中葉柄、葉面的變形數(shù)據(jù)，開展葉柄、葉面主葉脈撓曲線方程的參數(shù)辨識。利用建立的棉花葉片非線性大變形參數(shù)辨識模型，得到葉柄、葉面主葉脈撓曲線方程的最優(yōu)辨識參數(shù)分別為p=[1,2,3,4,5]=[?0.732 4, 5.178 2, ?4.445 8, 5.458 6, 3.056 8]、=[1,2,3,4,5,6,7]=[1.018 1, ?2.891 9, 0.987 3, 1.848 3, 3.592 0, 5.94, ?1.158 2]，葉柄、葉面主葉脈撓曲線預(yù)測值與仿真分析值的平均絕對百分誤差（Mean Absolute Percentage Error，MAPE）分別為5.13%、10.43%，決定系數(shù)2分別為98.99%、89.18%，模型預(yù)測值與仿真值的對比如圖4，說明建立的葉柄、葉面主葉脈撓曲線方程具有較好的精度。

2）模型驗證

利用高速攝像葉片變形試驗的變形量數(shù)據(jù)（見表3），驗證葉柄、葉面主葉脈撓曲線方程的有效性。對比分析可得（見表7），葉柄、葉面主葉脈撓曲線方程的最大預(yù)測誤差分別為8.5%、11.4%，驗證了所建葉片彎曲數(shù)學(xué)模型的有效性。

表7 葉片變形數(shù)學(xué)模型的誤差統(tǒng)計

3 棉花葉片迎風(fēng)面積變化規(guī)律分析

根據(jù)前文建立的棉花葉片變形模型，結(jié)合葉片變形高速圖像和耦合仿真所反映的葉片變形特點，利用葉柄撓曲線方程和葉面主葉脈撓曲線方程，分析輔助氣流、葉傾角等參數(shù)影響下葉片傾角和迎風(fēng)面積的變化規(guī)律。

3.1 葉傾角變化量

根據(jù)圖4葉柄和葉面撓度數(shù)據(jù)，因為葉柄彈性模量遠(yuǎn)高于葉面彈性模量，風(fēng)力載荷下棉花葉片的受力變形主要發(fā)生在葉面部位。同時，根據(jù)雙向流固耦合仿真試驗中葉片等效應(yīng)變及變形量分布，輔助氣流作用下，棉花葉面的彎曲變形主要發(fā)生在葉柄與葉面的連接處，葉片在連接處的彎折下沿風(fēng)速方向傾斜，但葉面仍基本保持平展?fàn)顟B(tài)，葉面彎度較小，如圖5所示。由此，可根據(jù)葉面形心撓度值、葉尖撓度值及初始葉傾角，明確棉花葉片撓度值與葉傾角變化量的幾何關(guān)系，如圖6所示，進(jìn)而可得風(fēng)力載荷作用下葉片葉傾角變化量Δ的計算公式，如式（13）所示。利用葉片變形流固耦合模型，可提取16組葉傾角變化量測量值；根據(jù)公式（13），可得對應(yīng)的16組葉傾角變化量的計算值。兩組數(shù)據(jù)的對比分析如圖7a所示，模型計算值與仿真測量值的趨勢一致，模型計算的MAPE為11.89%，說明葉片葉傾角計算模型的精度良好，可用于不同參數(shù)組合下葉片葉傾角計算。

3.2 葉片迎風(fēng)面積變化規(guī)律

由風(fēng)力載荷下棉花葉片變形特點可知，沿輔助氣流流向，棉花葉片葉傾角將發(fā)生變化，葉片迎風(fēng)面積是葉片面積在氣流流向法平面上的投影，其將隨葉傾角的變化而改變。為量化葉片迎風(fēng)面積的變化程度，定義葉片迎風(fēng)面積變化率Δ為變形前葉片迎風(fēng)面積與風(fēng)力載荷下變形后迎風(fēng)面積的比值，如公式（14）所示，其值大于1，代表迎風(fēng)面積增大，其值小于1，代表迎風(fēng)面積減小。

分析不同初始葉傾角及輔助氣流流速下葉片迎風(fēng)面積變化趨勢，如圖7b。對于初始葉傾角0°和10°，隨氣流流速增加，葉片迎風(fēng)面積逐步減小，霧滴輸運通道拓寬，有利于霧滴的穿透運移；但是，對于初始葉傾角20°和30°，隨氣流流速增加，葉片迎風(fēng)面積呈現(xiàn)先增加后減小的趨勢。上述數(shù)據(jù)表明，棉花葉片在風(fēng)力載荷作用下發(fā)生彈性變形，但其迎風(fēng)面積并非單調(diào)遞減，也會出現(xiàn)迎風(fēng)葉面積增加的情況。上述分析表明，速度恰當(dāng)?shù)妮o助氣流能夠減小葉片的迎風(fēng)面積，拓寬霧滴輸運通道；反之，則導(dǎo)致迎風(fēng)面積增加，冠層更加郁閉。具體而言，初始葉傾角為20°和30°時，應(yīng)保證棉花葉片的輔助氣流末速度分別不小于2.4和2.9 m/s，以實現(xiàn)減小迎風(fēng)面積、拓寬輸運通道的目標(biāo)。

4 結(jié) 論

本文基于棉花葉片流固耦合模型，建立了葉片變形參數(shù)辨識模型，并結(jié)合葉片變形特點，量化分析了棉花葉片的迎風(fēng)面積變化規(guī)律。

1）建立棉花葉片流固耦合仿真模型，經(jīng)高速攝像試驗驗證，模型最大誤差為9.85%，具有較好的仿真精度。

2）利用棉花葉片流固耦合模型，開展三因素四水平正交試驗分析，利用試驗數(shù)據(jù)分別建立基于非線性最小二乘法的葉柄、葉面主葉脈撓曲線參數(shù)辨識模型，模型平均絕對百分誤差分別為5.13%、10.43%。

3）結(jié)合葉柄、葉面主葉脈撓曲線方程，建立葉片葉傾角計算模型，模型平均絕對百分誤差為11.89%；進(jìn)而通過量化分析不同參數(shù)組合下葉片迎風(fēng)面積的動態(tài)變化過程，明確了初始葉傾角為20°和30°時，應(yīng)保證棉花葉片的輔助氣流末速度分別不小于2.4和2.9 m/s。

[1]田志偉，薛新宇，崔龍飛，等. 植保無人機(jī)晝夜作業(yè)的霧滴沉積特性及棉蚜防效對比[J]. 農(nóng)業(yè)工程學(xué)報，2020，36(5)：69-77.

Tian Zhiwei, Xue Xinyu, Cui Longfei, et al. Comparison of droplet deposition characteristics and cotton aphid control effect of plant protection UAV working during the day and night[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(5): 69-77. (in Chinese with English abstract)

[2]劉雪美，苑進(jìn)，張曉輝，等. 3MQ-600型導(dǎo)流式氣流輔助噴桿彌霧機(jī)研制與試驗[J]. 農(nóng)業(yè)工程學(xué)報，2012，28(10)：8-12.

Liu Xuemei, Yuan Jin, Zhang Xiaohui, et al. Development and experiment on 3MQ-600 type air-assisted boom sprayer with air-deflector[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(10): 8-12. (in Chinese with English abstract)

[3]張弼鋮，喬寧波. 3WP(X)-2000-18吊桿式噴桿噴霧機(jī)的研制[J]. 新疆農(nóng)機(jī)化，2018(2)：21-22，35.

Zhang Bicheng, Qiao Ningbo. Developing of 3WP(X)-2000-18 suspender type boom sprayer[J]. Agricultural Mechanization in Xinjiang, 2018(2): 21-22, 35. (in Chinese with English abstract)

[4]魏新華，邵菁，解祿觀，等. 棉花分行冠內(nèi)冠上組合風(fēng)送式噴桿噴霧機(jī)設(shè)計與試驗[J]. 農(nóng)業(yè)機(jī)械學(xué)報，2016，47(1)：101-107，90.

Wei Xinhua, Shao Jing, Xie Luguan, et al. Design and experiment of air-assisted cotton boom sprayer with separating row and spraying in inside and upper canopy[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(1): 101-107, 90. (in Chinese with English abstract)

[5]燕明德，賈衛(wèi)東，毛罕平，等. 風(fēng)幕式噴桿噴霧霧滴粒徑與速度分布試驗[J]. 農(nóng)業(yè)機(jī)械學(xué)報，2014，45(11)：104-110.

Yan Mingde, Jia Weidong, Mao Hanping, et al. Experimental investigation of droplet diameter and velocity distributions in air-assist boom sprays[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(11): 104-110. (in Chinese with English abstract)

[6]李龍龍，何雄奎，宋堅利，等. 果園仿形變量噴霧與常規(guī)風(fēng)送噴霧性能對比試驗[J]. 農(nóng)業(yè)工程學(xué)報，2017，33(16)：56-63.

Li Longlong, He Xiongkui, Song Jianli, et al. Comparative experiment on profile variable rate spray and conventional air assisted spray in orchards[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(16): 56-63. (in Chinese with English abstract)

[7]許童羽，于豐華，曹英麗，等. 粳稻多旋翼植保無人機(jī)霧滴沉積垂直分布研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報，2017，48(10)：101-107.

Xu Tongyu, Yu Fenghua, Cao Yingli, et al. Vertical distribution of spray droplet deposition of plant protection multi rotor UAV for japonica rice[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(10): 101-107. (in Chinese with English abstract)

[8]孫誠達(dá)，柳長鼎. 風(fēng)送噴霧霧滴冠層穿透模型構(gòu)建及應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報，2019，35(15)：25-32.

Sun Chengda, Liu Changding. Construction and application of droplet canopy penetration model for air-assisted spraying pattern[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(15): 25-32. (in Chinese with English abstract)

[9]Duga A T, Ruysen K, Dekeyser D, et al. Spray deposition profiles in pome fruit trees: Effects of sprayer design, training system and tree canopy characteristics[J]. Crop Protection, 2015, 67: 200-213.

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

[11]Endalew A M, Hertog M, Gebrehiwot M G, et al. Modelling airflow within model plant canopies using an integrated approach[J]. Computers and Electronics in Agriculture, 2009, 66(1): 9-24.

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

[13]Li J, Shi Y, Lan Y, et al. Vertical distribution and vortex structure of rotor wind field under the influence of rice canopy[J]. Computers and Electronics in Agriculture, 2019, 159: 140-146.

[14]劉希鳳. 植物葉片拓?fù)浣Y(jié)構(gòu)與其力學(xué)性能關(guān)系研究[D]. 廣州：華南理工大學(xué)，2010.

Liu Xifeng. Research on Relationship Between Topological Structure and Mechanical Property of Plant Leaf[D]. Guangzhou: South China University of Technology, 2010. (in Chinese with English abstract)

[15]周顯青，張玉榮，褚洪強，等. 糙米機(jī)械破碎力學(xué)特性試驗與分析[J]. 農(nóng)業(yè)工程學(xué)報，2012，28(18)：255-262.

Zhou Xianqing, Zhang Yurong, Chu Hongqiang, et al. Experiment and analysis of mechanical properties of mechanical crushing brown rice[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(18): 255-262. (in Chinese with English abstract)

[16]任鳳章，周根樹，趙文軫，等. 梁三點彎曲法測量薄膜彈性模量[J]. 稀有金屬材料與工程，2004，33(1)：109-112.

Ren Fengzhang, Zhou Genshu, Zhao Wenzhen, et al. Measurement of elastic modulus of films by three-point bend test[J]. Rare Metal Materials and Engineering, 2004, 33(1): 109-112. (in Chinese with English abstract)

[17]馬龍飛，李立州. 植物葉脈自然選擇過程的拓?fù)鋬?yōu)化模擬研究[J]. 應(yīng)用力學(xué)學(xué)報，2014，31(1)：132-136，12.

Ma Longfei, Li Lizhou. The topology optimization simulation study on the process of natural selection of leaf veins[J]. Chinese Journal of Applied Mechanics, 2014, 31(1): 132-136, 12. (in Chinese with English abstract)

[18]Menter F R, Langtry R, Volker S. Transition modeling for general purpose CFD codes[J]. Flow Turbulence and Combustion, 2006, 77(1/2/3/4): 277-303.

[19]Yang Y, Gu M, Chen S Q, et al. New inflow boundary conditions for modeling the neutral equilibrium atmospheric boundary layer in Computational Wind Engineering[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2009, 97(2): 88-95.

[20]洪毅. 高性能復(fù)合材料螺旋槳的結(jié)構(gòu)設(shè)計及水彈性優(yōu)化[D]. 哈爾濱：哈爾濱工業(yè)大學(xué)，2010.

Hong Yi. Structure Design and Hydroelastic Optimization of High Performance Composite Propeller[D]. Harbin: Harbin Institute of Technology, 2010. (in Chinese with English abstract)

[21]倪萬潮. 現(xiàn)代中國棉花生產(chǎn)技術(shù)[M]. 北京：中國農(nóng)業(yè)科學(xué)技術(shù)出版社，2011.

[22]戴奮奮. 風(fēng)送噴霧機(jī)風(fēng)量的選擇與計算[J]. 植物保護(hù)，2008，34(6)：124-127.

Dai Fenfen. Selection and calculation of the blowing rate of air-assisted sprayers[J]. Plant Protection, 2008, 34(6): 124-127. (in Chinese with English abstract)

[23]張軍，王一鳴，趙燕東，等. 基于橢球函數(shù)的棉花葉傾角分布動態(tài)模擬[J]. 農(nóng)業(yè)機(jī)械學(xué)報，2009，40(4)：157-160，183.

Zhang Jun, Wang Yiming, Zhao Yandong, et al. Dynamic simulation of leaf inclination angle distribution based on ellipsoidal function[J]. Transactions of the Chinese Society for Agricultural Machinery, 2009, 40(4): 157-160, 183. (in Chinese with English abstract)

[24]Roth N A, Uhl D, Mosbrugger V, Kerp H. Evolution and function of leaf venation architecture: A review[J]. Annals of Botany, 2001, 87: 553-566.

[25]劉靜靜. 典型植物葉片剛?cè)狁詈狭W(xué)特性及其仿生研究[D]. 長春：吉林大學(xué)，2018.

Liu Jingjing. Research on Rigid-flexible Coupling Mechanical Characteristics and Bionics of the Typical Plant Leaves[D]. Changchun: Jilin University, 2018. (in Chinese with English abstract)

[26]劉鴻文. 高等材料力學(xué)[M]. 北京：高等教育出版社，1985.

[27]Kane B, Smiley E T. Drag coefficients and crown area estimation of red maple[J] Canadian Journal of Forest Research, 2006, 36(8): 1951-1958.

[28] 陳晶. 一類非線性系統(tǒng)的參數(shù)辨識方法研究[D]. 無錫：江南大學(xué)，2013.

Chen Jing. Identification Methods for A Class of Nonlinear Systems[D]. Wuxi: Jiangnan University, 2013. (in Chinese with English abstract)

[29]梁禮明. 優(yōu)化方法導(dǎo)論[M]. 北京：北京理工大學(xué)出版社，2017.

[30]林恬，孟小紅，張致付. 基于約束最小二乘與信賴域的儲層參數(shù)反演方法[J]. 地球物理學(xué)報，2017，60(10)：3969-3983.

Lin Tian, Meng Xiaohong, Zhang Zhifu. The petrophysical parameter inversion method based on constrained least squares and trust region approach[J]. Chinese Journal of Geophysics, 2017, 60(10): 3969-3983. (in Chinese with English abstract)

Deformation modelling of cotton leaf under assisted airflow using bidirectional fluid-structure coupling method

Liu Xinghua1,2,Yuan Jin1,3, Li Yang1,3,Suo Jinzheng1, Liu Xuemei1,2※

(1.’271018,;2.’271018,;3.’271018,)

Assisted airflow can cause the deformation of canopy leaf for the wide transport channel of droplets in the air-assisted spray. The uniformity of droplet deposition can be enhanced in the crop canopy. In this study, a mathematical model of cotton leaf deformation was established under the action of assisted airflow using fluid-structure coupling and parameter identification. Firstly, the petiole and leaf surface samples of cotton were collected to obtain the size parameters. A three-point bending was selected to calculate the elastic modulus of petiole and leaf surface. The moisture content of the petiole and leaf surface was measured during drying at 105 ℃. The results showed that the median elastic modulus of leaf surface was 46.5 MPa, where 95% confidence interval was [28.5, 64.5] MPa, whereas the median elastic modulus of petiole was 244.6 MPa, where 95% confidence interval was [215.5, 273.69] MPa, while the average moisture content of petiole was 87.2%, where 95% confidence interval was [82.5%, 91.9%], and the average moisture content of leaf surface was 80.7%, where 95% confidence interval was [72.3%, 89.1%]. Then, a bidirectional fluid-structure coupling model was established to characterize the deformation process of cotton leaf. The deflection data were obtained in four monitoring points, including the middle and end of the petiole, as well as the center and tip of leaf surface. A high-speed camera was used to carry out the deformation test of cotton leaf with assisted airflow. The deflection test data of monitoring points were captured to verify the coupling model with the maximum simulation error of 9.85%. Furthermore, an orthogonal test was performed on experimental factors, including the assisted airflow speed, leaf inclination angle, elastic moduli of petiole and leaf. A significant order of experimental factors was obtained: wind speed, elastic modulus, leaf inclination angle. Finally, a mathematical model of cotton leaf deformation was constructed using the nonlinear least square. A trust-region iterative was used to obtain the identification parameters. A systematic evaluation was made on the performance of the mathematical model for cotton leaf deformation. The Mean Absolute Percentage Error (MAPE) of petiole deflection was 5.13%, and the MAPE of main vein deflection was 10.43%. It was found that the bending deformation of cotton leaf mainly occurred at the junction of petiole and leaf surface in the assisted airflow, where the leaf surface basically kept flat with the relatively small curvature. A quantitative analysis was made to obtain the initial parameters of blade inclination and dynamic variation in the frontal area of cotton leaf at different airflow speeds. The frontal area of cotton leaf decreased monotonously with the increase of airflow speed, when the inclination angles of cotton leaf were 0° and 10°. At the inclination angles of 20° and 30°, the frontal area of cotton leaf increased first and then decreased as the airflow speed increased. It infers that the canopy density can increase under the inappropriate assisted airflow. This finding can provide a sound reference to understanding the dynamic changes of cotton canopy density, further optimizing spraying parameters in the air-assisted spray.

assisted airflow; model; leaf deformation; fluid-structure coupling; parameter identification

劉興華，苑進(jìn)，李揚，等. 采用雙向流固耦合方法構(gòu)建輔助氣流作用下棉花葉片變形模型[J]. 農(nóng)業(yè)工程學(xué)報，2021，37(5)：68-76.doi：10.11975/j.issn.1002-6819.2021.05.008 http://www.tcsae.org

Liu Xinghua, Yuan Jin, Li Yang, et al. Deformation modelling of cotton leaf under assisted airflow using bidirectional fluid-structure coupling method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(5): 68-76. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2021.05.008 http://www.tcsae.org

2020-10-01

2021-01-13

國家自然科學(xué)基金項目（52075308）；山東省自然科學(xué)基金項目（ZR2019PC024）；山東省農(nóng)業(yè)重大應(yīng)用技術(shù)創(chuàng)新項目（SD2019NJ003）

劉興華，博士生，講師，研究方向為植保裝備與精準(zhǔn)噴霧理論。Email：lxh9357@163.com

劉雪美，博士，教授，博士生導(dǎo)師，研究方向為植保裝備與精準(zhǔn)噴霧理論。Email：lxmywj@126.com

10.11975/j.issn.1002-6819.2021.05.008

TP391.9

A

1002-6819(2021)-05-0068-09