Dynamic stratified porosity computation from canopy interaction simulation between airflow and leaves

The main goal of wind-driven spraying is to use assisted airflow to disrupt the structure of branches and leaves and broaden the air delivery channel, so as to achieve uniform droplet deposition in the middle and lower parts of the canopy. Due to the complex branch and leaf structure inside the canopy, there is currently no effective method to express the dynamic changes of canopy porosity and the law of airflow attenuation under assisted airflow. In this study, based on the two-way fluid-structure interaction numerical simulation method, the relating between the assisted airflow and the structural parameters of the cotton canopy is analyzed, and a new method for predicting and simulating the dynamic porosity of the canopy is proposed. Firstly, a two-way fluid-structure interaction model based on Lattice Boltzmann (LB) solver and Finite Element (FE) solver is developed to simulate the deformation motion of cotton leaves and the spatial distribution of airflow field, and the correctness of the numerical simulation is verified based on indoor measurement data. Secondly, the post-processing method of Computational Fluid Dynamics (CFD) is used to obtain images of leaves at different canopy positions under assisted airflow, and the porosity changes are calculated and analyzed by image processing. The research results show that under different initial wind speeds (5 m·s-1, 10 m·s-1, 15 m·s-1), the maximum normalized mean absolute error (NMAE) between the simulated values and the measured values is 13.99%, 20.72% and 16.08%, respectively. The coefficient of determination (R2) for linear fitting between simulated values and measured values is 0.9221. These validation results indicate the effectiveness of the numerical simulation method. The validated CFD model is applied to predict leaf deformation and porosity changes within the canopy under various wind loads and times. The application results have well revealed the interaction between crop leaves and airflow, and will be beneficial to make a better understanding of the effect of assisted airflow on droplet deposition.


Introduction
In the process of plant protection application, a large amount of pesticide misapplication will reduce the effectiveness of pesticide application and increase environmental pollution (Gil and Sinfort, 2005;Damalas and Eleftherohorinos, 2011;Tudi et al., 2021).The main goal of precision pesticide application is to achieve uniform coverage and deposition of pesticides in the target crop canopy (Khot et al., 2012;Li et al., 2021;Grella et al., 2022).The canopy characteristics of target plants directly affect the application mode and the droplet deposition effect (Duga et al., 2015;Xun et al., 2022).A complete understanding of the canopy characteristics of the target plant is important to evaluate the airflow velocity and turbulence levels within the canopy.
Canopy parameters are not only important indicators of growth and yield, but also important factors affecting pesticide interception and deposition (Olesen et al., 2003;Liu et al., 2020).In the late stage of crop growth, the stems and leaves of the plant population cover each other.Air-assisted sprayers can effectively deliver pesticides within dense canopies (Reichard et al., 1979;Derksen et al., 2008).The disturbance of the assisted airflow on the canopy branches and leaves can change the porosity of the canopy, thus widening the transmission channel of the pest control agent, which helps to achieve the droplet deposition at the lower canopy (Müller et al., 2018;Qiu et al., 2022).Scholars have conducted a large number of spray deposition field tests and wind tunnel tests on different target crops with different air-assisted sprayers.Cross et al. (2001aCross et al. ( , 2001b) ) conducted a series of field experiments on apple trees of different sizes to study the complex interaction between air-volume flow rate, spray-liquid flow rate, spray quality (droplet size distribution) and crop characteristics.Chen et al. (2013a;2013b) measured the spray retention and nontarget deposition at three crown growth stages (i.e.leaf stage, half leaf stage and full leaf stage), and found that increasing canopy density significantly reduced the amount of drift from the target.Duga et al. (2015) analyzed the spray deposition profiles in different pome fruit trees and concluded that tree characteristics such as total leaf cover, leaf wall porosity and tree volume strongly influenced total on-target deposition.These empirical spray studies indicate that spray deposition is caused by the complex interaction between the canopy and the air in canopy.
Crop spraying is a complex process involving the interaction of many parameters, such as pesticide dose and spray volume, sprayliquid distribution, droplet spectrum, air volume, sprayer speed, meteorological conditions and crop characteristics.In the past few decades, modeling approaches, especially computational fluid dynamics (CFD) models, have been effectively used to understand and characterize the crop spraying process (Badules et al., 2018;Zhang et al., 2018).Scholars adopted the averaging procedure to model airflow within a plant canopy without considering the flow details of individual elements (Wilson and Shaw, 1977;Raupach and Shaw, 1982).The canopy is considered as a porous medium to study the transport of airflow and droplets in canopy.The properties of porous medium are determined by the structural parameters of the canopy (porosity, leaf density) (Da Silva et al., 2006;Cui et al., 2022;Zhang et al., 2022).Endalew et al. (2009), 2010 simulated the effect of large branches and airflow by adding a resistance term in the canopy.The leaves in the canopy have a large effect, especially on the crop canopy.Dorr et al. (2008) developed a turbulence probability model, which can combine the motion model of fog droplets with the three-dimensional structure of plants, and simulate the drift of pesticide droplets around different plant structures.These studies further illustrate the influence of canopy structure on airflow distribution and droplet deposition.
Optical porosity is an important indicator for quantitatively measuring canopy structure parameters.It is defined as the ratio of leaf gap area on a projection plane to the contour area of the canopy leaves (Loeffler et al., 1992;Zhu et al., 2003).To accurately describe and calculate porosity, researchers have conducted extensive studies using hemispherical photography (Qu et al., 2016), laser point clouds (Escolà et al., 2017), hyperspectral or thermal infrared techniques (Neinavaz et al., 2016) and physical and mathematical models (Giusti and Marsili-Libelli, 2006).However, these sensorbased calculation methods are costly and can only calculate the static porosity of canopy.In addition, the spatial distribution of leaves in canopy cannot be fully captured due to high canopy leaf density and heavy shading, especially in the late stages of crop growth.In fact, crops cultivated in the field have traceable patterns in growth and spatial distribution of leaves (Guo et al., 2009).According to these morphological characteristics (Liu et al. (2020(Liu et al. ( , 2021a) ) specially studied the canopy porosity during spraying, proposed a 3D model to calculate the changes of canopy porosity, and realized the rapid prediction of crop canopy porosity.However, there is a strong interaction between leaves and airflow in the actual process of air-assisted spraying.The assisted airflow can not only move and deform the leaves, thus widening the droplet transport channel, but also improve the droplet transfer speed.It can effectively improve droplet deposition in the plant canopy and reduce drift (Panneton and Piche, 2005;Tang et al., 2021).Existing porosity calculation methods are difficult to express this physical process.
At present, the interaction between airflow and canopy is not clear.It is reflected in the following two aspects: one is how the airflow changes the deformation of the leaves, and the other is how the leaves affect the distribution of the airflow field.However, current research mainly calculates porosity based on the static structural characteristics of the canopy, and porosity is constantly changing in the actual spraying process.There is a lack of research on the quantitative description of porosity under assisted airflow, especially the dynamic changes in porosity.In our previous research, we studied a numerical simulation method for the fluidstructure interaction of leaves and spray airflow (Cui et al., 2023).In this study, we attempt to introduce the entire plant structure into numerical simulation and to calculate the real-time dynamic porosity.The relating between assisted airflow and canopy structural parameters is analyzed, and a new method for predicting and simulating canopy dynamic porosity is proposed.

3D virtual cotton plant model and artificial cotton
In this paper, in order to build a 3D virtual plant model with reasonable simplification of the canopy, the main stems, fruiting branches, petioles, and leaves are mainly considered.According to previous studies (Liu et al., 2021a;Cui et al., 2022) and cotton field measurement data (Figure 1), the structural parameters of cotton, such as plant height, petiole length and diameter, and leaf shape are determined.The main stem nodes of cotton are divided into 20 segments.The 5-20 nodes are set as fruiting branch growth positions.The multiaxial branching of fruiting branches is simplified to straight branches, and the distance between nodes on the one fruiting branch is the same.The ratio of ovate, 3-lobed and 5-lobed leaves to the total number of leaves in a single cotton plant is 20%, 55% and 25%, respectively.Based on the collected phenotype data, the windward areas of ovate, 3-lobed and 5-lobed leaves are determined to be 4206.25 mm 2 , 7859.10 mm 2 and 9823.00 mm 2 , respectively.
The 3D phenotypic plant models are constructed in SolidWorks (2019, Dassault Systemes, FR) software.The cotton plant has a height of 130cm and is stratified into three layers (i.e., 0-50 cm along the main stem height direction as the lower layer, 50-90 cm as the middle layer, and 90-130 cm as the upper layer).Based on the morphological characteristics and growth pattern of cotton plants, a 3D virtual plant is constructed as shown in Figure 2A.
An artificial plant is built according to the morphological characteristics of the 3D virtual plant.The material of artificial cotton leaves is selected from cloth material based on the previous experiments of leaf deformation and spray retention (Liu et al., 2021a;Liu et al., 2021b).The main stem is made of PVC plastic pipe.Wire of suitable stiffness is used for the leaf stem and fruiting branch materials.The fruiting branches and leaves are placed on the pre-drilled holes in the main stem.The single artificial cotton plant is shown in Figure 2B.

Numerical approach 2.2.1 Lattice Boltzmann model
Lattice Boltzmann (LB) model is suitable for solving many complex scientific problems.In particular, it does not require tracing the interface between different phases when dealing with multiphase and multi-component flows.In addition, it has been proven to have reliable accuracy in dealing with problems on microscopic and macroscopic scales (Ran and Xu, 2009;Men et al., 2017;Cui et al., 2023).
The important basis of the LB model is the theory of molecular motion.And LB model has the following assumptions.
a.The velocity component of each moving molecule is calculated without considering the influence of adjacent molecules.
b.The collision between two molecules is only considered.
c.The trajectory of each moving molecule is calculated without considering environmental factors.
Based on the above assumptions, the equation of molecular motion calculation function f is obtained.The independent variables of the equation are spatial velocity position vector, molecular velocity vector and time.The assumption on the molecular collision term can be simplified to a single relaxation time Bhatnagar Gross Krook (BGK) collision operator (Aidun and Clausen, 2010;Zhang et al., 2020).This simplification can effectively reduce the computational power.The simplified equation is the Boltzmann-BGK equation, expressed as (Aidun and Clausen, 2010): Where, r is the position vector; t is the time, s; f a is a discrete velocity distribution function; K a is discrete particle velocity vector; d t is the time step; t is the dimensionless relaxation time, t = t 0 =d t ; f eq a is a local equilibrium distribution function; F a is the external force term.
Through the discrete process, particles may move and collide, meaning that particles can move from one node to another in adjacent time steps while colliding with other adjacent particles.In addition, the LB method can calculate the macroscopic characteristics of the fluid through the statistical analysis of particles in the computational domain, and establish the relating between microscopic particles and macroscopic phenomena.The model D n Q m can be expressed as n dimensions and m discrete velocities.In this study, the octree lattice structure of D 3 Q 27 is used, as shown in Figure 3.

Turbulence model
Large Eddy Simulation (LES) is used to model the turbulence distribution.This approach introduces an additional viscosity, called turbulent eddy viscosity to model the sub-grid turbulence (Weickert et al., 2010).The LES scheme adopts a wall-adapting local eddy viscosity model, which provides a consistent local eddyviscosity and near-wall behavior (Ducros et al., 1998;Men et al., 2017;Zhang et al., 2020).The specific formulas are as follows: Where, U t , B w D, D, and j ab are turbulent eddy viscosity, filter scale, unit grid scale, and Kronecker symbol, respectively.Q ab and Q d ab are the resolving scale strain rate tensors; B w (0.325) is a constant.ɡ ab ɡ ba and ɡ rr are the components of the strain rate tensor obtained from the second-order moment via the LB model.In the above equations, the subscripts a, b and r denote directions in space.The µ and x are velocity at a given distance from the wall and local flow direction tangential to the wall.

Boundary conditions and computational domain 2.2.3.1 Solid domain
In the fluid-structure coupling collaborative simulation, the setting size of stem and leaf parameters and the shape quality of mesh subdivision will have a significant impact on leaf deformation.The data interaction between the solid domain and the fluid domain should be completed by setting appropriate boundary conditions.It is necessary to establish the coupling environment required for collaborative simulation to ensure the accuracy of the simulation.The elements are divided in Finite Element (FE) solver, and the stem and leaf parameters and boundary conditions are set.Through twoway fluid-structure coupling, the change trend and deformation amount of leaves can be calculated, and then the dynamic porosity of leaves after deformation can be obtained.
The accurate and reasonable segmentation of FE mesh is the basis of FE solver analysis.High-quality mesh can not only ensure reasonable analysis results, but also shorten the simulation time.Therefore, selecting a suitable element segmentation method is particularly important.The numerical integration method of FE solver adopts Gaussian numerical integration, and the integration points of different element shapes are different.The C3D8R element of linear hexahedron reduction integral is used for the element subdivision of stems, leaves and petioles, which is used as a threedimensional eight-node linear solid element, namely hexahedron element.Each node has six degrees of freedom and can bend in any direction.The combined element subdivision of cotton plant is shown in Figure 4.The number of elements is 16966, and the number of nodes is 41981 and the element size is 9mm.

A B
3D virtual cotton plant (A) and artificial cotton plant (B).
In the fluid-structure coupling analysis of assisted airflow and leaves, it is necessary to set the calculation attributes of the FE solver for the solid domain of stems and leaves.
(1) Material properties The bending of stems and leaves affected by the assisted airflow is an elastic deformation phenomenon, so the elastic material parameters are used to simulate the stems and leaves.According to previous studies (Ma and Li, 2014;Liu et al., 2021a;Cui et al., 2023), the elastic modulus of cotton leaves is set to 46.5Mpa, the Poisson's ratio of leaves is 0.32, and the density of leaves is 700 kg•m -3 .
(2) Analysis step settings In FE solver, the dynamic display analysis is set up.Since the spraying device drives the air curtain and nozzle to move in the actual spraying process.The initial value is verified according to the maximum time step of the fluid domain pilot site.The solid domain analysis step is set to 1s, and the field output adjusts 1s to 200 uniform time intervals, that is, there are 200 imaging effects within 1s.
(3) Boundary condition According to the distribution and connection of stems and leaves, the location of stem is defined as a completely fixed constraint in the boundary conditions setting, namely 0 degrees of freedom.In the simulation process, it only bears the influence of airflow and no external force, so it is only necessary to set the corresponding contact attributes.

Fluid domain
The fluid domain model is the air domain model.The plant should be placed in the air domain.The air domain model is built in 3D modeling software.The size is 1500 mm long, 1500 mm wide and 1500 mm high.This model is saved in a file format that can interface with the LB solver, and is named as 'air.stl'.The imported air fluid domain uses meshless modeling.The coordinate position of the fluid domain is adjusted to ensure that the branches and leaves are all in the flow field analysis domain.The location distribution of solid domain and fluid domain models is shown in Figure 5.The gravity acceleration applied to the fluid is set to -9.81m•s -2 .The default material 'Material 1' is used in the fluid domain, and its relevant parameters are set as follows: the molecular weight of air is 28.996 g•mol -1 , the density of air is 1.225 kg•m -3 , and the operating temperature is 289.35K (16.2°C).The gas flowing at low speed is a Newtonian fluid, so the dynamic viscosity is set to 1.7894e-05 Pa •s.In air-assisted spraying, the downward airflow can open the upper branches and leaves, which plays an important role in increasing the amount of droplets deposition in the lower and middle layers in dense canopy (Zhang et al., 2022;Zhu et al., 2022).In this study, we mainly consider the effect of downward airflow on porosity changes.

Fluid-structure coupling module
In FE solver, the contact surface with the fluid domain is set as the fluid-structure co-simulation boundary, as shown in Figure 6.After completing the above steps, set the calculation file in the Job module, name it Job-1, and export the.inp file.The recognition subroutine of fluid-structure coupling interface is set up in the.inp file.The flowchart of the fluid-structure coupling approach is shown in Figure 6.

Validation of numerical simulation
In order to make the experimental conditions controllable and avoid interference of uncontrollable factors in the natural environment, in this paper, we construct an artificial cotton plant based on the 3D virtual cotton plant (Dekeyser et al., 2014;Cui et al., 2022).To validate the numerical simulation, the differences between simulated and experimental values at the same position in the canopy are compared.The positions of the sampling points are  In simulation, the airflow velocity at the sampling point is output through the detection line.The normalized mean absolute error (NMAE) between the measured and the simulated values of the upper, middle, and lower layers of the canopy were compared.The total difference was analyzed by fitting a linear equation.

Dynamic changes in canopy porosity based on image processing
Optical porosity is defined as the ratio of canopy leaf void area to contour plane area on a projection plane (Loeffler et al., 1992;Zhu et al., 2003).It is related to plant density, structure and environmental conditions, and is a structural parameter closely related to flow and resistance characteristics near plants.In pesticide spraying, leaves are the main organ of droplet deposition.In addition to the size and position of the leaf, the leaf will bend and deform due to airflow disturbance, and the optical porosity will change accordingly.
Based on the 3D model of the plant and the fluid-structure coupling process, we propose a method to calculate the target canopy porosity by layering and zoning based on image processing.This method can obtain canopy dynamics images during the interaction between canopy and airflow through CFD post-processing, and then use image processing to calculate porosity at any canopy height.The specific process is as follows: (1) Image acquisition.Using the layer plane as the reference plane, images of different canopy positions of threedimensional cotton plant targets along the reference plane are captured and stored in LB solver (Figures 8A-C).

FIGURE 6
Flow chart of the fluid-structure coupling approach.

A B FIGURE 7
Distribution of sampling points (A) and indoor airflow measurement (B).
(3) Determination of canopy outer edge.Determine the canopy projection along the outermost leaf edge of the canopy and calculate the projection area (S i ).
(4) Calculation of windward area.The pixel value occupied by the plant projection is counted, and according to the image and 3D plant leaf scale, the leaf windward area (A i ) at different canopy positions of the plant is calculated.
(5) Calculation of stratified porosity.The porosity of each layer can be calculated using the following formula: 3 Results and discussions

Simulation validation results
In this study, in order to verify the reliability of the numerical simulation results, several sampling points are set up in the upper, middle, and lower layers of the artificial cotton planting target area.The initial airflow velocity values are 5 m•s -1 , 10 m•s -1 and 15 m•s -1 , respectively.The results are shown in Figure 9A.Under the initial air velocity of 5 m•s -1 , the normalized mean absolute error (NMAE) between the simulated value and the measured value of the upper, middle and lower parts of the canopy are 7.54%, 12.51% and 13.99%, respectively.Under the speed of 10 m•s -1 , the NMAE are

A B FIGURE 9
The comparative analysis between the measured values and the simulated values.(A) is the NMAE between simulated and measured values; (B) is the linear fitting analysis between simulated and measured values.
9.36%, 13.14% and 20.72%, respectively.Under the speed of 15 m•s -1 , the NMAE are 8.69%, 10.54% and 16.08%, respectively.The results show that the simulated airflow velocity distribution can reflect the attenuation of airflow velocity in the canopy.To further verify the accuracy of the simulation results, linear fitting analysis is performed on all simulated values (S) and measured values (M), as shown in Figure 9B.The expression of the linear equation is M=0.9849S+0.0246,and the coefficient of determination (R 2 ) is 0.9221.The fitting results indicate that the fluid-structure coupling method can effectively reflect the dynamic changes of canopy porosity.

Leaf deformation and porosity changes in the canopy at different airflow velocities
Figure 10 shows the leaf deformation and porosity changes within the canopy under the influence of different airflow velocities at 0.1s.Under 0 m•s -1 , 5 m•s -1 , 10 m•s -1 and 15 m•s -1 , the upper, middle and lower canopy leaves have different degrees of deformation, and the airflow velocity and vortex distribution changes are different.
In the absence of wind (0 m•s -1 ), we export images at different heights of the canopy and calculated the porosity of the upper, middle, and lower layers as 49.95%, 59.85% and 49.16%,  Cui et al. 10.3389/fpls.2023.1238360Frontiers in Plant Science frontiersin.orgrespectively (Figure 10A).When the airflow velocity reaches 5 m•s -1 , we find that the leaves at different heights showed slight twisting or bending deformation (Figure 10B).Compared with the velocity of 0 m•s -1 , the porosity at the same canopy position is 52.27%, 52.23% and 49.74%, respectively.Figure 10C shows the changes in leaves within the canopy at the same time when the airflow speed increases to 10 m•s -1 .The porosity at different canopy heights is 66.83%, 64.16% and 38.71%, respectively.As the flow rate increases, the plant shrinks its shape, rolls up its leaves and bends downstream, resulting in a decrease in its cross-sectional area and an increase in its fluidization degree.Figure 10D shows the changes in leaves within the canopy at the same time when the airflow speed increases to 15 m•s -1 , with the porosity of each layer being 68.07%, 73% and 62.35%, respectively.As the wind load increases, the leaves sway more and become more random.Compared with the initial state of 5 m•s -1 , the change ratio has increased by 36.28%,21.97% and 26.83%, respectively.From the perspective of leaf deformation, when subjected to increased airflow load, the distribution of canopy branches and leaves will be reconfigured to reduce wind force and absorb momentum in the airflow, which is a reconstruction phenomenon (Ennos et al., 2000;Vollsinger et al., 2005;Kane et al., 2008;Miri et al., 2018).This behavior is often observed in plants with flat and thin leaves (Gillies et al., 2002).At higher wind speeds, plant resistance decreases and leaf deformation increases.However, in actual spraying operations, the distribution of vortices around the canopy changes greatly, which increases the drift of droplets in the atmosphere.We derive the windward area and porosity at different canopy relative heights (CRHs) (10%, 20%,…, 100%) to further describe the changes in leaf deformation and porosity along the depth of the canopy.It is assumed that the initial porosity is defined as the porosity at different positions of the canopy under the condition of no wind (0 m•s -1 ).The dynamic variation of porosity is defined as the change in porosity relative to the initial porosity at different relative heights of the canopy under wind conditions (5 m•s -1 , 10 m•s -1 , 15 m•s -1 ). Figure 11A shows the changes in windward leaf area at different relative heights of the canopy under different initial velocities.From Figure 11A, it can be seen that the changes in leaf area at the bottom (CRH=10%) and top (CRH=90%) of the canopy are relatively small.While the changes in windward area in the middle of the canopy (CRH=20% -80%) are relatively large.Under low wind speed (5 m•s -1 ), the maximum change in the windward area of leaves at different CRHs is 20%.As the airflow speed increases, the changes in the windward area of leaves within the canopy under 10 m•s -1 and 15 m•s -1 are significantly larger than those under low wind speed conditions.These results indicate that low wind speeds cause slight vibration of the leaves, while high wind speeds cause significant bending or deformation of the leaves.Figure 11B shows the dynamic changes in porosity at each relative height of the canopy.Under the three wind speed conditions of low (5 m•s -1 ), medium (10 m•s - 1 ), and high (15 m•s -1 ), the relative heights of the canopy with the largest changes in porosity are 40%, 70%, and 40%, respectively, and the maximum changes in porosity are 11%, 26%, and 16%, respectively.The maximum change in porosity occurs in the middle of the canopy because the interaction between airflow and leaves becomes stronger in areas with denser leaves.Comparing Figures 11A, B, we find that there is consistency between the position with the largest change in leaf windward area and the position with the largest change in porosity.The specific modification form is as follows: 'the position with the largest change in porosity.Figures 11C-D and show the changes in velocity and vorticity at different relative heights of the canopy.As the depth of the canopy increases, the airflow velocity continues to decay.At CRH=80%, a significant change in vorticity occurred, which is related to the obstruction of airflow by the upper branches and leaves of the canopy.

Dynamic changes in porosity at different times
Figure 12 shows the dynamic changes in porosity at different times (0-1s) in the upper, middle and lower layers of the canopy at a wind speed of 5 m•s -1 (light breeze).At 0.1s, the upper leaves of the canopy begin to deform first due to the influence of the upper airflow, and the upper layer undergoes changes earlier than the middle and lower layers.From 0.2s to 0.5s, the leaves within the canopy continuously bend downward and deform.From the perspective of the curvature amplitude of most leaves in the canopy, the curvature degree of the canopy leaves is the highest at 0.6s.After 0.7s, the leaves gradually rebound.Due to the obstruction of airflow by the upper and middle branches and leaves within the canopy, the intensity of airflow decreases greatly when it reaches the lower part.The airflow ultimately causes only slight deformation of the leaves in the lowest layer, indicating that under the action of a small airflow (light breeze), the deformation amplitude of the canopy leaves is small and they sway slightly up and down.
Figure 13 shows the deformation of internal branches and leaves in the canopy under airflow disturbance at different times when the assisted airflow velocity is 10 m•s -1 (moderate wind).At 0.1s, the airflow has just reached the top of the canopy.At 0.2s to 0.4s, the canopy leaves begin to gradually bend downward.Compared to the initial wind speed of 5 m•s -1 , the leaf reaches its maximum deformation at 0.5s.At this moment, the maximum porosity of the entire canopy indicates that the moment when the wind speed increases and the maximum leaf deformation occurs is relatively late.After 0.6s, the deformation amplitude of most leaves in the canopy shows a slight rebound.A small number of leaves have a large rebound amplitude and even show a reversal trend, especially the leaves at the bottom of the canopy have obvious deformation (0.9s to 1.0s).This may be caused by the vortices formed under the canopy.Due to the obstruction of the upper leaves, the assisted airflow under the canopy gradually decreases, and the intensity of the vortices on the back of the lower leaves is greater than that of the assisted airflow on the front.This leaf inversion is very important for the uniform deposition in the middle and lower parts of the canopy during spraying.
Figure 14 shows the deformation of internal branches and leaves in the canopy under airflow disturbance at different times when the assisted airflow velocity is 15 m•s -1 (high wind speed).At 0.2s, the airflow has reached the lower part of the canopy and caused deformation of the lowest leaves.As the airflow gradually enters the interior of the canopy, the leaf deformation amplitude reaches its maximum value at 0.4s, which is earlier than the maximum moment of leaf deformation under medium (10 m•s -1 ) and low (5 m•s -1 ) wind speeds, indicating that the leaves are susceptible to deformation under high wind speeds.After 0.5s,  Cui et al. 10.3389/fpls.2023.1238360Frontiers in Plant Science frontiersin.orgthe leaves in the canopy become disorderd, floating up and down, and some of the leaves undergo reversals, especially for 0.9s and 1.0s.From the perspective of leaf deformation, the degree of leaf deformation will not increase as the airflow velocity continues to increase.However, a higher velocity can lead to an increase in airflow disturbance outside the canopy, which leads to the possibility of spray droplets drift.

Conclusions
In this study, a new method based on CFD simulation and image processing is proposed to calculate the dynamic changes in porosity caused by leaf deformation under assisted airflow.A twoway fluid-structure interaction model is developed based on LB solver and FE solver.The model achieves a quantitative and intuitive analysis of the deformation of leaves in the canopy under the action of assisted airflow.The fluid-solid interaction model is validated by indoor experiments.The results of CFD postprocessing are analyzed using image processing algorithms, and the stratified dynamic porosity in the canopy is calculated at different velocities and different times.This study provides an idea to clarify the dynamic changes of porosity in canopy during air-assisted spraying and to analyze the mechanism of increasing droplet deposition in canopy by assisted airflow.
Compared with previous static porosity studies, this study can better reveal the dynamic interaction phenomenon between crop leaves and airflow.However, the developed CFD model still needs further improvement.In order to accurately predict the dynamic porosity of the canopy, 3D leaf modeling should consider more details.In addition, the effect of different directional assisted airflow on the canopy structure should be considered.

FIGURE 3
FIGURE 3Schematic diagram of the LB model of D 3 Q 27 .

FIGURE 4
FIGURE 4Schematic diagram of elements distribution in the structural explicit FE solver.

FIGURE 5
FIGURE 5Schematic diagram of the distribution of cotton plant and air domain locations.
FIGURE 8 Calculation process of canopy stratified porosity.(A) The front view of a stratified cotton plant.(B) The front view of the upper layer.(C) The top view of the upper layer.(D) The image processing result of the top view.

FIGURE 12
FIGURE 12Changes in canopy leaf morphology at different times under 5 m•s -1 .

FIGURE 13
FIGURE 13Changes in canopy leaf morphology at different times under 10 m•s -1 .

FIGURE 14
FIGURE 14Changes in canopy leaf morphology at different times under 15 m•s -1 .