Phenotype analysis method of tree under strong wind stress

Abstract

The invention discloses a phenotype analysis method for trees under strong wind stress, which comprises the following steps: a phenotype platform comprises a fan and trees, and the trees are arranged in front of the fan; the laser scanner collects point cloud data of the tree at a certain wind speed; modeling a canopy of the target tree by taking leaf clusters in the point cloud data as units to obtain a crown graphic model, and calculating phenotypic characteristic parameters of each leaf cluster; and calculating wind characteristic parameters in the crown through the k-epsilon turbulence model, the constructed crown graphic model and a numerical solution. According to the invention, the limitation that a conventional wind tunnel experiment is high in cost and small in size and cannot meet the normal scale of an experiment tree is overcome, so that a vivid open wind environment can be reproduced; a method based on a K-epsilon turbulence model special for the forest is proposed to solve wind characteristic parameters in a forest canopy; a digital rubber tree model is reconstructed through a point cloud driving technology, and changes of phenotypic characteristics of an experimental tree under a certain wind speed condition are visually displayed.

Description

Technical Field

[0001] The invention relates to the field of tree phenotyping, in particular to a method for phenotyping trees under strong wind stress. Background Art

[0002] The rubber tree (Hevea brasiliensis) is of great economic importance due to its high yield of natural rubber, which is an important commodity in numerous industries. In Asia, rubber trees are usually planted in coastal areas where typhoons occur every year, bringing strong winds, heavy rains and storm surges that can cause devastating effects. Since rubber trees are not inherently resistant to strong winds, wind forces them to undergo a series of profound changes in phenotypic characteristics and physiological processes, resulting in consequences beyond mere cosmetic damage. Some studies have shown that when the average wind speed is below 2 m / s, it stimulates transpiration in rubber trees, enhances the absorption of inorganic nutrients by the roots, speeds up metabolic processes, and helps to dry the bark faster, thereby reducing diseases on the tapping surface. On the contrary, if the average wind speed exceeds 2 m / s, excessive foliar transpiration of water occurs, which may damage tender leaves and affect the growth of rubber trees and the flow of latex. In the absence of adequate windbreaks, average wind speeds between 8 and 17 m / s can adversely affect the growth rate and wood quality of rubber trees, leading to problems such as bent branches, asymmetrical crowns, and skin damage in young trees. In addition, winds of a certain intensity or speed exceeding 17 m / s can cause leaf tearing and falling, broken branches or trunks, and even occasionally uprooting of trees.

[0003] Tree phenotyping and environmen-tation studies involve assessing the responses of trees to various climatic environments and different cultivation conditions to provide guidance for tree risk prevention and selection of stress-resistant varieties, including those in strong wind fields. In fact, strong winds are often accompanied by gloomy skies, thunder and heavy rain, which makes it particularly challenging to measure tree phenotypic variation during heavy rains. Therefore, many researchers have designed alternative strategies, such as industrial wind tunnels, electronic sensors mounted on trees, and computer simulation algorithms to provide arrangements for rapid exploration of complex wind-tree interaction mechanisms.

[0004] One challenge of using wind tunnels is that they are usually small in size due to their complex operation and high maintenance costs. Since it is difficult to put a mature whole tree into a wind tunnel, a common approach is to use a miniature replica of the tree's size, made of various materials, such as multi-scale models based on polymers and metal materials or 3D printed fractal models according to the rules of L-systems. These models are used to infer the possible wind-induced tree structure reconfiguration and canopy-induced wind energy transfer. When real trees are used, these are usually small, young, or pruned to fit the size of the wind tunnel. These methods have been used to evaluate the safety of green roof systems, the ventilation performance of urban residential areas, and the wind resistance function of forest windbreaks. In addition, electronic sensors are often used to monitor the movement of trees under different wind conditions, for example, strain gauges are attached to the trunk to measure stress distribution and local deformation; wireless tree sensing systems and three-axis accelerometers are used to detect the tilt angle of the trunk; and two-axis geomagnetic compasses and cameras equipped with 3D magnetic trackers are used to monitor the swaying behavior of the canopy.

[0005] Computer simulations offer an alternative approach to assess the aerodynamic parameters of trees in wind conditions. The simulated trees can be treated as physical models or porous media and subjected to wind dynamics. Typically, physical models include a mass-spring system to simulate the movement of leaves, a flexible cantilever beam or a simple pendulum model for conifers, rigid and flexible branching patterns to assess the mechanical stability of the tree, and various anchoring conditions for the root-soil system. In addition, close-range laser scanning technology is used as an auxiliary tool to characterize the crown structure, using a simplified approach that considers the crown as a geometric primitive of a porous medium and is coupled with computational fluid dynamics to predict the evolution of the three-dimensional wind field within the crown.

[0006] Following the above experiments, rigid and flexible trees that exhibit phenotypic variation under wind stress were also documented. For example, the crowns of trees reorient and reorganize their canopies and shed needles, leaves, and branches to reduce drag. Higher porosity leads to greater reorganization of canopy elements, minimizing failures through active adjustments. At the same time, as the leaf area index increases, the area of ​​the windward surface decreases, while the streamline effect reduces wind drag. As wind speed increases beyond the elastic limit of the tree, clustered plastic deformations of various crowns align with the prevailing wind direction. In addition, short wind pulses persistently change the morphology of roots, such as an increase in root length and surface area, while the density of root tissue decreases. On the other hand, some studies have shown that trees reduce the impact of wind on the trunk and roots through structural movements (such as shock absorption or vibration), helping them to disperse energy more efficiently. Studies have reported that as wind loads increase, local torsional vibrations of the crown are weakened, while bending vibrations are more pronounced throughout the tree. There is also some literature that there is a strong correlation between wind loads imposed by coherent structures in turbulent flows and the displacement of spruce trunks. Trees vibrate under wind forces until the roots or trunk fail, and the failure mechanisms are similar between species, but the specific details may vary.

[0007] Although we have made significant progress in the study of tree crown behavior under wind, existing methods still have many limitations. First, due to the limitation of wind tunnel size, only small artificial or pruned trees with a height of less than 2 meters can be used, which cannot truly reflect the scale of natural trees. At the same time, the airflow behavior caused by the closed structure of the wind tunnel is different from that in the open natural environment. Secondly, there is currently a lack of intelligent, automated and refined descriptions of tree shape changes, as well as methods to reveal the gas-tree coupled fluid equations and their aerodynamic parameter extraction mechanisms under different wind stress conditions. Finally, current technology mainly relies on cameras to record the fluctuations of trees disturbed by wind, but it is still challenging to provide a more vivid and comprehensive depiction of wind-driven tree behavior. Summary of the invention

[0008] The technical problem to be solved by the present invention is to provide a phenotypic analysis method for trees under strong wind stress in view of the deficiencies of the above-mentioned prior art. The phenotypic analysis method for trees under strong wind stress overcomes the high cost and small size limitations of traditional wind tunnel experiments that cannot meet the normal scale of experimental trees. An on-site phenotyping platform is established, which is equipped with large industrial fans and medium-sized trees, and can truly reflect the scale of natural trees to reproduce a realistic open wind environment; a laser scanner is used to record the changes in the crown under different continuous wind forces; a tree graphic model is designed to automatically decompose the crown into multiple leaf clusters, and the phenotypic characteristics of each leaf cluster are calculated based on the scanning points; a method based on a forest-specific K-epsilon turbulence model is proposed, and the constructed tree graphic model is combined with a numerical solution to solve the aerodynamic parameters (i.e., wind characteristic parameters) in the forest canopy; finally, a digital rubber tree model is reconstructed through point cloud driven technology to intuitively display the changes in the phenotypic characteristics of the experimental tree under certain wind speed conditions.

[0009] In order to achieve the above technical objectives, the technical solution adopted by the present invention is:

[0010] A method for phenotyping trees under strong wind stress specifically comprises the following steps:

[0011] Step 1: Build a field phenotyping platform, which includes a fan and target trees. The target trees are arranged in front of the fan. Start the fan to maintain a certain wind speed.

[0012] Step 2: The laser scanner collects point cloud data of the target tree under a certain wind speed;

[0013] Step 3: Model the canopy of the target tree based on the leaf clusters in the point cloud data, obtain the canopy graphic model, and calculate the phenotypic characteristic parameters of each leaf cluster;

[0014] Step 4: Establish a k-ε turbulence model for forest environment;

[0015] Step 5, calculate the wind characteristic parameters in the tree canopy through the k-ε turbulence model, the constructed tree canopy graphic model and the numerical solution;

[0016] Step 6: Establish a digital twin model based on the wind characteristic parameters within the crown and the phenotypic characteristic parameters of each leaf cluster to describe the wind-induced response of the tree under a certain wind speed.

[0017] As a further improved technical solution of the present invention, in step 1, the fan includes an industrial forced ventilation fan, and the industrial forced ventilation fan is connected to the lifting truss through a lifting device, and the lifting device is used to adjust the height of the industrial forced ventilation fan. There are 3 target trees in total, and the 3 target trees are arranged in a straight line in front of the fan.

[0018] As a further improved technical solution of the present invention, the step 3 is specifically:

[0019] The crown is divided into leaf clusters according to the affiliation of the first-level branches, and the point cloud data corresponding to each first-level branch is marked with different colors. The crown layer is then modeled based on leaf clusters to obtain a crown graphic model, in which each leaf cluster corresponds to a first-level branch. The crown graphic model contains phenotypic characteristic parameters extracted from the point cloud, including the volume and leaf area density LAD of each leaf cluster.

[0020] As a further improved technical solution of the present invention, the step 4 is specifically as follows:

[0021] A k-ε turbulence model for forest environment is established, and a unified simplified mathematical expression of the k-ε turbulence model is obtained:

[0022]

[0023] in Represents the air flow velocity vector in three-dimensional space, that is U i represents the air flow velocity along the x direction; U j represents the air flow velocity along the y direction; U k represents the air flow velocity along the z direction; φ represents the relevant aerodynamic parameters, including U i , U j , U k , turbulent kinetic energy κ and turbulent dissipation ε; q φ represents the energy source term; B φ Indicates an additional term;

[0024] The momentum equation, turbulent kinetic energy κ equation and turbulent dissipation ε equation are derived based on the unified simplified mathematical expressions;

[0025] The momentum equation is:

[0026]

[0027] Where, μ represents the air viscosity coefficient, μ t represents the turbulent kinematic viscosity, ρ represents the air density, C d represents the drag coefficient of the forest; A represents the leaf area density of the leaf cluster; κ represents the turbulent kinetic energy, P represents the fluid pressure in the calculation space; |U| represents the absolute value of the total velocity of the three velocity components;

[0028]

[0029]

[0030]

[0031]

[0032]

[0033]

[0034] and All of them conform to the Einstein summation convention;

[0035]

[0036] Among them, U1 is equivalent to the air flow velocity along the x direction, U2 is equivalent to the air flow velocity along the y direction; U3 is equivalent to the air flow velocity along the z direction;

[0037] The turbulent kinetic energy κ equation is:

[0038]

[0039] Among them, β d =4, indicating the amount of energy loss when interacting with obstacles; σ κ =1.0, indicating the turbulent Prandtl number; ε indicates the turbulent dissipation;

[0040]

[0041]

[0042] It represents the deformation rate or equivalent strain rate and conforms to the Einstein summation convention;

[0043]

[0044] The equation for turbulent dissipation ε is:

[0045]

[0046] Among them, σ ε =1.31, c ε1 =1.44, c ε2 =1.92,α p =1.7,α d =1.7,β p =1;

[0047]

[0048]

[0049] As a further improved technical solution of the present invention, the step 5 is specifically as follows: meshing the air flow area containing the tree crown graphic model, using a numerical solution and combining formula (4), formula (7) and formula (11) to calculate the U in each grid unit in the tree crown i , U j , U k , turbulent kinetic energy κ, turbulent dissipation ε, and the pressure difference between two grid cells.

[0050] As a further improved technical solution of the present invention, the step 5 is specifically as follows:

[0051] The unified simplified mathematical expression, that is, the integral expression of formula (2) is as follows:

[0052]

[0053] Where V represents the volume of each grid cell;

[0054] Convert formula (12) into the following form:

[0055]

[0056] represents the outward pointing unit normal vector relative to the surface of each grid cell; S represents the sum of the areas of the six faces on each grid cell;

[0057] After discretizing formula (13), the flow field is converted into a grid structure, which is specifically expressed as follows:

[0058]

[0059] where δ x represents the distance between the centers of adjacent grid cells along the x direction, δ y represents the distance between the centers of adjacent grid cells along the y direction, δ z represents the distance between the centers of adjacent grid cells along the z direction, represents the aerodynamic parameters of the right side surface through the grid cell (α, β, λ); represents the aerodynamic parameters of the left side of the grid cell (α, β, λ); represents the aerodynamic parameters of the upper side of the grid cell (α, β, λ); represents the aerodynamic parameters of the lower side through the grid cell (α, β, λ); represents the aerodynamic parameters passing through the outer surface of the grid unit (α, β, λ); represents the aerodynamic parameters passing through the inner surface of the grid unit (α, β, λ); φ α，β，λ represents the aerodynamic parameters passing through the grid unit (α, β, λ); U1，α，β，λ represents the air flow velocity through the grid unit (α, β, λ) along the x direction, U 2，α，β，λ represents the air flow velocity through the grid unit (α, β, λ) along the y direction; U 3，α，β，λ represents the air flow velocity through the grid unit (α, β, λ) along the z direction; s right represents the area of ​​the right side of the grid unit (α, β, λ), s left represents the area of ​​the left side of the grid cell (α, β, λ), s upper represents the area of ​​the upper side of the grid unit (α, β, λ); s below represents the area of ​​the lower side of the grid unit (α, β, λ); s out represents the area of ​​the outer surface of the grid unit (α, β, λ); s inner represents the area of ​​the inner side of the grid unit (α, β, λ);

[0060] Convert formula (14) into the following form:

[0061]

[0062] where ζ is the adjacent coefficient derived from equation (14), Represents the corresponding variable set in each grid unit, and the corresponding variables include U i , U j , U k , turbulent kinetic energy κ and turbulent dissipation ε;

[0063] Decompose ζ into the sum of the diagonal matrix D, the lower triangular matrix B, and the upper triangular matrix J, that is, ζ = D + B + J. The mathematical derivation is as follows:

[0064]

[0065] The weighted Jacobi iteration method is used, where the weight is ω, and formula (16) is iteratively propagated as follows:

[0066]

[0067] Where ω = 2 / 3, according to the relationship ζ-D = +B + J, equation (17) is expressed as follows:

[0068]

[0069] The converged solution of formula (18) can be determined as the spectral radius of the iteration matrix is ​​less than 1, that is, ζ(I-ωD -1ζ)<1; Finally, the inlet wind speed of the air flow area after meshing is used as the initial value, and a series of gradually improved approximate solutions are generated through iterative calculation to extract the U in each grid cell in the tree crown. i , U j , U k , turbulent kinetic energy κ, turbulent dissipation ε, and the pressure difference between two grid cells.

[0070] The beneficial effects of the present invention are:

[0071] In the present invention, a phenotypic analysis method for trees under strong wind stress is proposed, combining point cloud analysis with digital twin technology. The present invention overcomes some limitations of existing methods and provides a practical solution to meet this challenge. The key points include: First, the custom-designed industrial fan system is able to reproduce variable airflow conditions, such as gusts and intermittent winds, thereby accurately simulating the pressure exerted on trees by real typhoons. Secondly, the crown modeling is performed based on the leaf cluster characteristics in the point cloud, and a fine-grained description is achieved based on the phenotypic characteristics to capture the dynamic morphological changes of the crown at different wind speeds. In addition, by synergistically considering aerodynamics and mechanical principles, we use the forest-specific k-ε turbulence equation, combined with the constructed tree model and numerical solution to calculate the wind characteristic parameters within the crown, and verify some of the calculation results through field measurements. Finally, a digital twin model was established, which can combine the wind characteristic parameters obtained in the tree crown, the phenotypic characteristic parameters of each leaf cluster and the obtained plant mechanical characteristics to depict the actions of the target trees under various wind forces, thereby providing comprehensive and intuitive information for understanding their biophysical processes under wind-induced effects. Compared with the existing technology that relies on cameras to record the fluctuations of wind-disturbed trees, the digital twin model can provide a more vivid and comprehensive description of wind-driven tree behavior. The method of the present invention realizes intelligent, automated and refined description of tree shape changes, and also provides a method for revealing the gas-tree coupling fluid equation and its aerodynamic parameter extraction mechanism under different wind stress conditions.

[0072] The method of the present invention establishes a field phenotyping platform equipped with a large customized industrial fan system to reproduce a realistic open wind environment. The platform is combined with transplanted medium-sized rubber trees equipped with sensors to detect their swaying behavior in different wind environments, which is low in overall cost, effectively overcoming the high cost and small size limitations of traditional wind tunnel experiments that cannot meet the normal scale of experimental trees. The method of the present invention uses a terrestrial laser scanner to record the changes in the tree crown under different continuous wind forces. The method of the present invention designs a computer graphics algorithm that can automatically decompose the tree crown into multiple leaf clusters and calculate the phenotypic characteristics of each leaf cluster based on the scan points. In addition, the method of the present invention proposes a method based on a forest-specific K-epsilon turbulence model, combining the constructed tree graphics model with a numerical solution to solve the aerodynamic parameters within the forest canopy. Finally, the digital rubber tree model is reconstructed through point cloud driven technology to intuitively display the changes in the phenotypic characteristics of the experimental trees under different wind speed conditions, providing a more vivid and comprehensive depiction of wind-driven tree behavior. This provides a benchmark implementation plan for simulating tree behavior under high wind stress conditions in the future and designing adaptive management strategies to reduce the impact of wind damage.

[0073] In summary, our method effectively quantifies the unstable behavior of rubber trees by elucidating the complex interaction between wind and tree shape. Although this method is limited to small-scale experiments on specific tree species, its detailed depiction of tree organs and in-depth exploration of the technical details of wind-induced damage are noteworthy. This method and its in-depth analysis constitute the important highlights of this invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 Schematic diagram of the field phenotyping platform and data collection methods for demonstrating tree responses to wind.

[0075] Figure 1 (a) is a schematic diagram of two large fans installed on a steel frame structure.

[0076] Figure 1 (b) shows a schematic diagram of the RIEGL laser scanner set up next to the experimental tree to obtain the corresponding point cloud data during wind processing.

[0077] Figure 1 (c) is a schematic diagram of the inertial motion sensor and mini anemometer fixed to the upper, middle and lower parts of the tree trunk, respectively, to record the aerodynamic parameters during the experiment in real time.

[0078] Figure 1 (d) shows the overall layout of the fan experiment site, where three rubber trees are exposed to the industrial fan system. At the same time, three laptops are connected to sensors tied to branches through wired connections to record the swinging behavior of the tree canopy under the influence of wind.

[0079] Figure 1 (e) is a schematic diagram of point cloud data of three static trees collected using laser scanning technology.

[0080] Figure 1 (f) in the figure is a schematic diagram of dividing the tree crown into leaf clusters according to the affiliation of the first-level branches, and marking the point cloud data corresponding to each first-level branch with different colors.

[0081] Figure 2 A programming flow chart showing how to use numerical methods to divide the experimental domain into smaller units to solve the governing equations of fluid dynamics.

[0082] Figure 2 (a)-(c) in the figure are schematic diagrams showing the flow field grid generation pattern obtained by laser scanning point mapping at the leaf cluster scale for tree models with different morphologies and structures.

[0083] Figure 2 (d) in Figure 2 A partial enlarged view of (a).

[0084] Figure 2 (e) is a schematic diagram showing the local neighboring grid used when solving the aerodynamic equations within the unit.

[0085] Figure 3 Schematic diagram of experimental trees subjected to continuous strong winds at different wind speeds (0m / s, 15m / s, 30m / s and 45m / s).

[0086] Figure 3 (a1) is a schematic diagram showing the experimental photos taken when the trees were subjected to wind mechanical loads under the condition that the fan blade rotation speed was 0 m / s.

[0087] Figure 3 (b1) is a schematic diagram showing experimental photos taken when the trees were subjected to wind mechanical loads under the condition of a fan blade rotation speed of 15m / s.

[0088] Figure 3 (c1) is a schematic diagram showing experimental photos taken when the trees were subjected to wind mechanical loads under the condition of a fan blade rotation speed of 30m / s.

[0089] Figure 3 (d1) is a schematic diagram showing the experimental photos taken when the trees were subjected to wind mechanical loads under the condition of the fan blade rotating speed of 45m / s.

[0090] Figure 3(a2) is a schematic diagram of point cloud data showing a relatively stable compression state reflected by the trees at a certain moment during the wind action when the fan blades rotate at a speed of 0 m / s.

[0091] Figure 3 (b2) is a schematic diagram of point cloud data showing a relatively stable compression state reflected by the trees at a certain moment during the wind action when the fan blades rotate at a speed of 15m / s.

[0092] Figure 3 (c2) is a schematic diagram of point cloud data showing a relatively stable compression state reflected by the trees at a certain moment during the wind action when the fan blades rotate at a speed of 30m / s.

[0093] Figure 3 (d2) is a schematic diagram of point cloud data showing a relatively stable compression state reflected by the trees at a certain moment during the wind action when the fan blades rotate at a speed of 45m / s.

[0094] Figure 3 (a3) is based on Figure 3 The point cloud data in (a2) shows a schematic diagram of the segmentation point cloud of each leaf cluster that constitutes the tree crown.

[0095] Figure 3 (b3) is based on Figure 3 The point cloud data in (b2) shows a schematic diagram of the segmentation point cloud of each leaf cluster that constitutes the tree crown.

[0096] Figure 3 (c3) is based on Figure 3 The point cloud data (c2) in the figure shows a schematic diagram of the segmentation point cloud of each leaf cluster that makes up the tree crown.

[0097] Figure 3 (d3) is based on Figure 3 The point cloud data (d2) in the figure shows a schematic diagram of the segmentation point cloud of each leaf cluster that makes up the tree crown.

[0098] Figure 4 Overview of inertial motion sensors fixed at different heights on the tree trunk to detect trunk vibrations.

[0099] Figure 4 (a)-(b) in the figure show the actual experimental scene diagrams of three target trees.

[0100] Figure 4 (c)-(e) in the figure are vibration detection curves of tree trunks at different heights along the main wind direction when the fan blade rotation speed is 15m / s.

[0101] Figure 4(f)-(h) in the figure are vibration detection curves of tree trunks at different heights along the main wind direction when the fan blade rotation speed is 30m / s.

[0102] Figure 4 (i)-(k) in the figure are vibration detection curves of tree trunks at different heights along the main wind direction when the fan blade rotation speed is 45m / s.

[0103] Figure 5 To show the cross-section of the interaction between wind and trees, the influence of the shape of the virtual tree model in the simulated wind field on the aerodynamic parameters is revealed along the wind direction.

[0104] Figure 5 (a1) in the figure shows the tree model diagram constructed based on the fan blade rotation speed of 15m / s.

[0105] Figure 5 (b1) in the figure shows the tree model diagram constructed based on the fan blade rotation speed of 30m / s.

[0106] Figure 5 (c1) in the figure shows the tree model diagram constructed based on the fan blade rotation speed of 45m / s.

[0107] Figure 5 (a2) in the figure shows the instantaneous wind speed distribution between tree crowns when the fan blades rotate at a speed of 15 m / s.

[0108] Figure 5 (b2) in the figure is a diagram showing the instantaneous wind speed distribution between tree crowns when the fan blades rotate at a speed of 30 m / s.

[0109] Figure 5 (c2) in the figure is a diagram showing the instantaneous wind speed distribution between tree crowns when the fan blades rotate at a speed of 45 m / s.

[0110] Figure 5 (a3) shows the pressure gradient diagram caused by the deformation of the leaf cluster under the action of asymmetric stress when the fan blade rotates at a speed of 15m / s.

[0111] Figure 5 (b3) shows the pressure gradient diagram caused by the deformation of the leaf cluster under the action of asymmetric stress when the fan blade rotates at a speed of 30m / s.

[0112] Figure 5 (c3) in the figure shows the pressure gradient diagram caused by the deformation of the leaf cluster under the action of asymmetric stress when the fan blade rotates at a speed of 45m / s.

[0113] Figure 5(a4) shows the turbulent kinetic energy (TKE) between tree canopies when the fan blades rotate at a speed of 15 m / s.

[0114] Figure 5 (b4) shows the turbulent kinetic energy (TKE) between tree canopies when the fan blades rotate at a speed of 30 m / s.

[0115] Figure 5 (c4) shows the turbulent kinetic energy (TKE) between tree canopies when the fan blades rotate at a speed of 45 m / s.

[0116] Figure 6 A highly realistic 3D tree model image generated for different wind excitation conditions.

[0117] Figure 6 (a) shows the experimental photo taken when the fan blade rotates at a speed of 15 m / s, the virtual tree digital twin model diagram, and a partial enlarged view.

[0118] Figure 6 (b) shows the experimental photos taken when the fan blades rotate at a speed of 30 m / s, the virtual tree digital twin model diagram, and a partial enlarged view.

[0119] Figure 6 (c) shows the experimental photos taken when the fan blades rotate at a speed of 45 m / s, the virtual tree digital twin model diagram, and a partial enlarged view. DETAILED DESCRIPTION

[0120] The specific embodiments of the present invention are further described below according to the accompanying drawings:

[0121] 1. Abstract:

[0122] Rubber trees in coastal habitats are often subjected to high-intensity wind stress. In order to study and predict the phenotypic characteristics of rubber trees under strong winds (the phenotypic characteristics under strong winds include wind characteristic parameters within the crown and phenotypic characteristic parameters of each leaf cluster), this study developed an integrated approach that collaborates with algorithms and hardware. The framework includes the following parts: (1) A custom-designed industrial fan is used to reproduce variable airflow fields under wind speeds of 15, 30, and 45 m / s, and a ground laser scanner and multiple motion sensors are used to obtain point cloud and vibration data; (2) A crown graphics model is developed based on point cloud data. The model constructs an approximate crown morphology through leaf clusters and extracts the phenotypic characteristics of rubber trees under the action of aerodynamic drag; (3) By combining the constructed tree model with gridded wind field fluid partitioning, a k-ε turbulence model specifically for forest environments is developed to calculate the wind characteristic parameters of the forest crown; (4) A digital twin model is established, which integrates detailed crown morphology with the gridded wind field fluid partitioning. The wind characteristic parameters within the model tree, the phenotypic characteristic parameters of each leaf cluster and the plant mechanical properties (plant mechanical properties, as described in Section 3.4 below) were integrated into the physical-virtual integration by defining mechanical properties for the photosynthetic and non-photosynthetic components of the model tree. For example, elastic moduli were assigned to green wood (5.4-8.9 GPa) and branches (2.2-4.6 GPa); in addition, the masses of fresh leaves (0.57-0.81 g) and senescent leaves (1.1-1.7 g) were determined; the elastic-gravity number of rubber leaves (0.11-0.58) was calculated, and the petiole bending stiffness value reflecting its flexibility was established (2.6-12.3 N·m 2 ); these values ​​define a biomechanical model that includes flexible pinnate leaves, elastically deformed petioles, and branches that bend under buckling stress) to describe the wind-induced response of the target trees under different wind conditions. The results show that trees with extended crowns are more likely to produce larger swing amplitudes than trees with compact crowns, but the crown volume of trees directly exposed to the wind shrinks significantly more than that of trees in sheltered areas. Inside the crown, when the inlet wind speed doubles, the crown compression effect is enhanced (the crown volume decreases by approximately 17%), the pressure gradient increases by 2.1 times, and the turbulent kinetic energy increases by approximately 60%. In addition, the adaptive behavior of the experimental trees under wind disturbance is intuitively restored using digital twin technology, providing a complete behavioral dataset for subsequent data-driven decision-making. In summary, this study proposes a comprehensive methodology that can analyze the phenotypic characteristics of trees under wind damage environments, and has potential for application in phenotypic analysis or environmental phenotyping studies for wind resistance assessment of rubber trees.

[0123] This study aims to develop a system for comprehensive evaluation of rubber tree motion under wind, combining actual measured data with model development to achieve complete simulation. First, a field phenotyping platform equipped with a large custom industrial fan system was established to reproduce realistic open wind environments. The platform was combined with transplanted medium-sized rubber trees equipped with sensors to detect their swaying behavior in different wind environments, effectively overcoming the high cost and small size limitations of traditional wind tunnel experiments that cannot meet the normal scale of experimental trees. Secondly, a terrestrial laser scanner was used to record the changes in the tree crown under different continuous wind forces. A computer graphics algorithm was designed to automatically decompose the tree crown into multiple leaf clusters and calculate the phenotypic characteristics of each leaf cluster based on the scanned points. In addition, a method based on a forest-specific K-epsilon turbulence model was proposed to combine the constructed tree graphics model with a numerical solution to solve the aerodynamic parameters within the canopy (i.e., wind characteristic parameters within the canopy). Finally, a digital rubber tree model was reconstructed using point cloud-driven technology to intuitively display the changes in the phenotypic characteristics of the experimental tree under different wind speed conditions. This provides a baseline implementation plan for simulating tree behavior under future high wind stress conditions and designing adaptive management strategies to reduce wind damage impacts.

[0124] 2. Research methods:

[0125] 2.1. Research area and data collection:

[0126] The experimental base was established in Hainan Province (3.30°-20.07°N, 108.15°-120.05°E). The region is located in East Asia, adjacent to the western coast of the Pacific Ocean, and has a tropical monsoon climate with an average annual temperature between 23℃ and 29℃. Hainan is also known for its abundant precipitation and abundant sunshine, making it an ideal area for rubber tree cultivation. As an important source of latex and timber, rubber trees are widely planted in this region, with a planted area of ​​330,800 tons as of 2020. From May to September, the weather in Hainan is hot and rainy, accompanied by heavy rainfall and thunderstorms. From June to mid-November, typhoon activity increases, bringing strong winds and heavy rains. Under normal circumstances, wind speeds can reach 15m / s, while during typhoons, wind speeds can reach as high as 30 to 60m / s. Wind damage is an important cause of death or loss of rubber tree yield. The potential risk of economic losses prompted us to gain a deeper understanding of the specific mechanisms of wind damage and its mitigation strategies for its impact on trees.

[0127] The field phenotyping platform includes an experimental device, which consists of two forced ventilation fans (industrial fans), a steel lifting truss and several steel wire ropes. In order to simulate the typhoon environment as realistically as possible, we selected two forced ventilation fans (see Figure 1), because of its ability to manufacture large fans. Each forced ventilation fan is model HTF-I-(A)-11A, with a diameter of 1.5m, a weight of 0.5t, a maximum power of 15kW, a maximum blade speed of 1450r / min, and an air flow rate of 50,000m 3 / h. The heights of the two forced ventilation fans were adjusted individually to ensure that the fan outlets were aligned with the experimental trees, and a cubic steel lifting truss with a length and width of 1.5m and a height of 6m was used. Four sets of steel brackets were fixed on four pillars, and the other ends were secured to the ground by steel nails. A pulley was installed on the top of the steel frame, supplemented by an iron chain, so that the height of the induced ventilation fan could be manually adjusted. In addition, an intelligent power switch timer was programmed to automatically open or close every 40 seconds, thereby generating gusty or continuous strong airflows at fixed time intervals to simulate the scenario of typhoons periodically attacking trees. The experiments were conducted in a closed studio to avoid the influence of external environmental airflow on the results.

[0128] Figure 1 Experimental setup and data collection methods to demonstrate the response of trees to wind forces. Figure 1 (a) shows two large fans installed on a steel frame structure. Figure 1 (b) shows that the RIEGL laser scanner is set up next to the experimental tree to obtain the corresponding point cloud data during wind processing. Figure 1 (c) shows that inertial motion sensors and mini anemometers are fixed to the upper, middle and lower parts of the tree trunks, respectively, to record the aerodynamic parameters in real time during the experiment. (d) The overall layout of the fan experiment site, with three rubber trees exposed to the industrial fan system. At the same time, three laptops are connected to sensors tied to branches through wired connections to record the swaying behavior of the tree canopy under the influence of wind. Figure 1 (e) in the figure shows the point cloud data of three stationary trees collected using laser scanning technology. Figure 1 (f) in the figure shows that the crown is segmented into leaf clusters according to the affiliation of the first-level branches, and the point cloud data corresponding to each first-level branch is marked with different colors.

[0129] Three trees were transplanted to the experimental site and firmly fixed to the ground after filling the soil. Figure 1As shown in (d), from right to left, they are named tree 1, tree 2 and tree 3, and are arranged in a straight line in front of the fan. Tree 1 is located closest to the fan, while tree 3 is farthest from the fan. The distance between the three trees is 2 meters, and the distance between tree 1 and the industrial fan system is also 2 meters. The height of these trees is about 6 meters, which is equivalent to the height of the steel lifting truss. Throughout the experiment, we deliberately adjusted the height of the fan to ensure that the artificial wind covers the full vertical range of each tree canopy. At the same time, by analyzing the typhoon track and wind speed data of the past 50 years, we found that once a typhoon enters a terrestrial forest area, its wind speed is usually reduced to less than 50 meters / second due to terrain, vegetation, mountains and other surface friction resistance. Therefore, we set and maintained three target wind speeds for the experiment by precisely adjusting the rotation speed of the fan blades: 15m / s, 30m / s and 45m / s.

[0130] Recording vibration data and point cloud data from tree trunks is essential for studying the behavior of trees under wind loads, a task accomplished by inertial sensors and laser scanners, respectively. The inertial motion sensors used are WT901 WIFI. These sensors are capable of recording three-dimensional acceleration, angular acceleration, and angle changes in real time. Sensors were installed at the bottom (1m), middle (3m), and top (5m) of three tree trunks to capture wood motion data in real time. The time series data was resampled to 2Hz, which is sufficient to capture trunk displacement and vibration without excessive resolution. In addition, all inertial motion sensors are susceptible to measurement pseudo-offsets, which result in bias errors. We filter this bias using a high-pass Butterworth filter to achieve noise suppression.

[0131] LiDAR point cloud data were collected using a RIEGL VZ-6000 terrestrial laser scanner. The scanner has a laser emission frequency of 300kHz and a measurement accuracy of 2mm. Two sets of medium-resolution (0.1m point spacing, 0.057 degree angular resolution) scanners were placed symmetrically on opposite sides of the three trees at a distance of 4m from the middle crown, pointing to the middle tree. The laser scanner usually completed the scan of the experimental trees in about 20 seconds. The wind fan was controlled by an intelligent power switch timer to apply continuous wind resistance so that the tree structure remained in a suppressed state. Under constant wind speed, point cloud data corresponding to the suppressed state of the tree shape were recorded to further characterize the phenotypic characteristics of the trees and their fine morphological modeling. By fusing the point cloud data acquired during the experiment, a computerized approximate model (by computer graphics to depict the changes of trees at the leaf cluster scale, that is, the compression state of trees at the leaf cluster scale under different wind conditions) was constructed to simulate the compression state of trees under different wind speeds. The detailed modeling steps will be provided in the subsequent subsections.

[0132] 2.2. Tree modeling at leaf cluster scale:

[0133] Point cloud data obtained from terrestrial laser scanning needs to be preprocessed for use in computer simulations. The subsequent air-solid interaction model is very time-consuming to simulate the movement of individual leaves and branches, so a balance must be struck between resolution and computational simplicity. The model models the canopy in leaf clusters, each of which corresponds to a primary branch. The model contains parameters extracted from the point cloud to specify the location, morphology, volume, and leaf area density of the leaf clusters.

[0134] First, a leaf classification-based algorithm was used to extract the trunk of each tree in the control and treatment conditions, and a small amount of manual correction was performed. Subsequently, a spatial colonization algorithm was used to generate a 3D trunk model based on the connectivity from the root to the end of each branch, and the first-level branch points were used as the center points of the leaf clusters. Finally, the k-means clustering algorithm was applied to determine the affiliation between leaf points and branch points to determine the association between corresponding leaf points in each cluster of leaves.

[0135] The volume of the leaf cluster is obtained by applying the alpha shape algorithm to the leaf point cloud of each first-level branch (see Figure 1 (f) in Figure 3). The computational scheme uses a global polynomial level set parameterization to smooth the closed surface of each leaf cluster, generating a convex hull appearance and a closed surface that fits the outline of the leaf cluster from the set of sampled points.

[0136] The leaf area density (LAD) of each leaf cluster is a key parameter that measures the total leaf area per unit volume and is intended to describe the characteristics of each tree crown. Higher LAD values ​​indicate greater wind resistance, associated with increased windward surface and higher wind resistance area, resulting in increased air flow resistance, changes in streamlined profiles, and increased turbulent kinetic energy within the canopy. In contrast, lower LAD values ​​promote smoother air passages and stable flow behavior.

[0137] The symbol "A" represents LAD, which is defined as:

[0138]

[0139] Among them, N in Represents the total number of intercepts recorded as scan points in the convex hull representing the corresponding leaf cluster, that is, the number of scan points generated in the convex hull. N tot Represents the total number of laser beams that reach the convex hull, which can be derived from the scanning configuration (including the angular resolution of the laser beam and the distance between the scanner and the target convex hull). It represents the average path length of the laser beam actually passing through the convex hull. tot and The quantification of can be achieved by ray tracing algorithm. G(θ) is the leaf projection function, which quantifies the projection coefficient per unit leaf area on the plane perpendicular to the beam direction θ. Here, G(θ) is assigned a mean value of 0.5 to represent the random (spherical) orientation of leaves within the crown.

[0140] Each leaf cluster is considered as a pseudo-turbulent medium containing the calculated leaf distribution density (LAD) and the fitted leaf cluster morphology, and is included in the CFD model to achieve a more realistic simulation and improve the prediction accuracy of the air flow pattern and microclimate conditions in the vegetation area. Before the wind treatment, the leaf cluster parameters of the three experimental trees and their components are listed in Table 1.

[0141] Table 1. Quantitative description of individual crowns and leaf clusters of the three stationary experimental trees before wind treatment. Specifically, tree 1 had the shortest distance to the fan, while tree 3 was the farthest from the fan.

[0142]

[0143]

[0144] Note: LAD is leaf area density.

[0145] 2.3 Aerodynamic model of forest environment:

[0146] Under strong wind conditions, the complex structure of the tree canopy with its uneven distribution of leaves, branches and trunks significantly disturbs the ordered laminar flow through it. This disturbance leads to the generation of turbulence, which manifests itself as irregular fluctuations in wind speed and direction, characterized by vortices, swirls or other irregular motion patterns. In general, the flow of air through trees can usually be regarded as an incompressible flow state because the wind speed is less than 100m / s and the Mach number is less than 0.3. Based on this, we propose a modified k-epsilon model for forest environments, which introduces an additional term in the framework of the governing equations. The k-epsilon model is often used to simulate flows at higher Reynolds numbers, where turbulence is the main feature and the molecular viscosity effect is relatively negligible. The unified simplified mathematical description of our modified k-epsilon model is as follows:

[0147]

[0148] in, Represents the air flow velocity vector in three-dimensional space, that is Specifically, U i represents the air flow velocity along the x direction, U j represents the air flow velocity along the y direction, U krepresents the air flow velocity along the z direction, i, j, k represent the components in three directions, that is, along the x, y, and z directions respectively. φ represents the relevant aerodynamic parameters, and q φ It represents the energy source term, and its specific definition is shown in Table 2. φ is an additional term used to describe the reverse energy interaction between each leaf cluster and contains information derived from the leaf area density and geometry. In addition, the K-epsilon equation uses an ensemble average solution method to calculate the average parameters and consider transient effects. Therefore, the K-epsilon equation does not involve the time factor at all in its expression and omits the first term on the left side of formula (2), that is, We then generated the following table format, where each project variation shows the aerodynamic model we designed.

[0149] Table 2. This study uses a unified equation and an improved K-epsilon model to conduct an in-depth analysis of the distribution of wind-related parameters within the tree canopy. The parameters used in equation (2) are as follows:

[0150]

[0151]

[0152] Note: '-' means the corresponding term does not exist. Index variables that appear twice in a monomial, such as i or j, follow the Einstein summation convention, i.e., all indices in the term are summed. i It is an abbreviation, and can also be understood as the speed in a certain direction (such as x, y, z, etc.).

[0153] The continuity equation describes the principle of mass conservation in a pipe: the mass flowing into the pipe minus the mass flowing out of the pipe equals the mass change. Therefore, substituting the items corresponding to the second row of Table 2 into equation (2), the continuity equation (3) can be derived as follows:

[0154]

[0155] The momentum equation is a statement of Newton's second law and involves the sum of the fluid or air forces acting on the tree unit, which causes the change in the rate of change of momentum. As mentioned in Section 2.2, each leaf cluster that constitutes the canopy is regarded as a porous medium with specific properties, such as position, volume, and LAD. Therefore, the drag effect of the canopy in the airflow domain can be considered using the additional source terms defined in Table 2. Subsequently, the terms corresponding to the third row of Table 2 are substituted into equation (2) to derive our momentum equation (4), which is as follows:

[0156]

[0157] where μ = 1.81*10 -5 kg / (m·s) is the air viscosity coefficient, μ t =c1κ 2 / ε is the turbulent kinematic viscosity, ρ=1.204kg / m 3 Represents the air density, C d is the drag coefficient of the forest, and its value varies from 0.15 to 0.36 depending on the density of the leaves, and is linearly related to the LAD of each leaf cluster. A represents the leaf area density of each leaf cluster derived from equation (1), P represents the fluid pressure in the calculation space; κ represents the turbulent kinetic energy; U represents the absolute value of the total velocity of the three velocity components; the absolute value of the total velocity is equal to the square root of the sum of the squares of the three components.

[0158] In equation (4) and The expanded form is as follows:

[0159]

[0160] and The expanded form of is similar to formula (5), and The expanded form of is similar to formula (6);

[0161] and All of them conform to the Einstein summation convention;

[0162]

[0163] Among them, U1 is equivalent to the air flow velocity along the x direction, U2 is equivalent to the air flow velocity along the y direction; U3 is equivalent to the air flow velocity along the z direction;

[0164] The fourth and fifth rows in Table 2 show the K–epsilon (k-ε) turbulence model used in computational fluid dynamics (CFD) to simulate the average flow characteristics under turbulent conditions. The specific formula for the turbulent kinetic energy κ shown in the fourth row of Table 2 is as follows:

[0165]

[0166] In equation (7), the additional source term where β d =4 represents the amount of energy lost when interacting with an obstacle, assumed to represent a branch or leaf in the canopy. κ =1.0 is the turbulent Prandtl number in the k transport equation, and ε represents the turbulent dissipation.

[0167] In formula (7) and The expanded form of is explained as follows:

[0168]

[0169] in It represents the deformation rate or equivalent strain rate and conforms to the Einstein summation convention;

[0170] Using the Einstein summation convention, the derived expression is as follows:

[0171]

[0172] For the turbulent dissipation ε shown in the fifth row of Table 2, the specific formula is as follows:

[0173]

[0174] In formula (11), and The expanded form of is the same as that in formula (7) and Therefore, the corresponding expansion form in mathematics will not be repeated. ε1 =1.44 and c ε2 =1.92 is the coefficient of the corresponding item. p =1 indicates the proportion of the blade shroud resistance that converts the average flow kinetic energy into the energy generated by the blade shroud components (i.e. the flow of the fluid around the blade shroud). ε =1.31 is the turbulent Prandtl number in the ε transport equation, α p and α d They represent the adjustment system of this model respectively and are set to 1.7 here.

[0175] 2.4. Gridding of wind-forest field:

[0176] Analytical solutions to the aerodynamic model equations of air flow through complex tree canopy structures are difficult to obtain due to the presence of nonlinear dynamics and complex high-order derivatives. Therefore, we use numerical methods to approximate solutions that tend to the exact solution in the fluid-structure interaction system. In physical simulations, such as finite element analysis or computational fluid dynamics, meshing is a standard method. The meshing of the experimental site in this study aims to discretize the air flow area containing the fan and tree components to determine the solid-gas coupling interaction between adjacent grids.

[0177] Commonly used flow field mesh generation software include ANSYS ICEM CFD and GAMBIT. However, these software tools are mainly targeted at rigid and regular geometries, such as cars, buildings, and pipelines. For trees with 3D geometric complexity, mesh generation often faces problems such as mesh distortion and incomplete coverage caused by specific details of the tree structure. Therefore, after the initial mesh data is generated by existing software, it is necessary to perform additional mesh optimization steps. To this end, we developed a self-written MATLAB code for mesh optimization in the forest field, with the goal of achieving mesh smoothing, refinement, and topology improvement.

[0178] The grid size is set to an average side length of 10 cm. Figure 2 As shown, a comprehensive grid is constructed in the simulated flow field area (the simulated flow field area includes the tree model we constructed). Figure 2 (a)-(c) in which the pink vertical area represents the wind inlet (also known as the wind inlet), indicates the location of the fan outlet, and assumes a constant average wind speed along the x-axis as the specified fan output speed. Figure 2 (d) shows Figure 2 A local enlarged view of (a), where the continuously distributed geometric grid is clearly visible. Figure 2 (e) in Figure 1 presents the mesh set used for the numerical solution, including the central mesh and its six neighboring meshes, where α, β, and γ represent the spatial arrangement of these meshes. An adaptive meshing strategy is used to take into account shape and structural features. Large leaf clusters and open spaces are subdivided into larger mesh cells to reduce the computational burden, while smaller mesh cells are used at the boundary junctions between leaf clusters and branches to capture fine-grained details more accurately. In the airflow scenario, the trunk part is set to be solid and no air circulation is allowed, while the leaf cluster part is given the properties derived from equation (1).

[0179] Figure 2 A programming flow chart showing how to use numerical methods to divide the experimental domain into smaller cells to solve the governing equations of fluid dynamics. Figure 2 (a), (b) and (c) in the figure respectively show the flow field grid generation mode obtained by laser scanning point mapping at the leaf cluster scale for tree models with different shapes and structures (15m / s, 30m / s, 45m / s). The pink on the left wall in each flow field represents the wind speed inlet (i.e. the airflow inlet), and the opposite wall is the wind speed outlet (i.e. the airflow outlet), and the wind speed inlet is the air inlet generated by the fan. Figure 2 (d) in Figure 2 A magnified partial view of (a) showing the geometric mesh of the tree crown and the surrounding void area. Figure 2 (e) in the figure shows a schematic diagram of the local neighboring grid used when solving the aerodynamic equations within the unit.

[0180] 2.5. Numerical solution method of the aerodynamic model in this study:

[0181] After subdividing the fluid domain into multiple cells, we further quantified the distribution of aerodynamic-related parameters in the experimental area, such as the airflow velocity, pressure difference (the pressure P is set to an initial value, and then it changes at each iteration. The interpolation of different grids gives P, and the pressure difference is the difference between the P calculated for each grid and another grid), turbulent kinetic energy, and turbulent dissipation for each grid cell in the canopy. Due to the nonlinear characteristics and complex high-order derivatives in aerodynamic models, it is usually difficult to obtain analytical solutions. Therefore, numerical analysis, a mathematical method for solving complex equations, is used to meet these challenges. Numerical methods transform complex equations into discrete numerical problems and use computer programs to obtain approximate solutions. Nonlinear terms are effectively expanded through continuous grid geometry to achieve the generation of stable solutions to turbulence models.

[0182] First, define the inlet and outlet boundaries of the airflow for each flow field (i.e., the inlet and outlet are boundary conditions), and then iteratively solve the difference equation to calculate the discrete solution of the κ-ε equation group, i.e., the numerical solution of formula (4), (7), and (11). The equations of formula (4), (7), and (11) can be converted in sequence according to formula (12)-(18) recorded in Section 2.5, that is, the corresponding numerical values ​​can be obtained. The numerical solutions of formula (4), (7), and (11) include the initial airflow velocity and initial pressure P, and solve for κ and ε. The initial value of the grid in the pink area is set to the inlet wind speed, and then at each iteration step, the wind propagation and the corresponding k and epsilon values ​​can be calculated.

[0183] By first specifying the inlet and outlet boundary conditions for each grid cell, air flows in or out, and then iteratively solving the difference equations to obtain the discrete solution of the turbulence equations. In the entire spatial range, all parameters of the four equations shown in Table 2 need to be solved to complete the integration of the fluid domain. Formula (2), the general form of the quaternary equation, is expressed as follows:

[0184]

[0185] Where V represents the volume of each grid cell mentioned in Section 2.4, and Represents the outward pointing unit normal vector relative to the surface of each grid cell. According to the divergence theorem, and the RANS equation as the average motion equation of fluid flow without considering the time-related terms, formula (12) can be converted into the following form:

[0186]

[0187] Equation (13) represents the surface integral of the vector field of the closed surface S on each grid cell, which is equal to the internal area V of the closed surface S (e.g. Figure 2 The closed surface S represents the sum of the areas of the six faces on each grid cell.

[0188] After discretizing equation (13), the flow field is converted into a grid structure, and the aerodynamic model is also converted into a discrete form, such as Figure 2 The specific expression is as follows:

[0189]

[0190] where δ x represents the distance between the centers of adjacent grid cells along the x direction, δ y represents the distance between the centers of adjacent grid cells along the y direction, δ z represents the distance between the centers of adjacent grid cells along the z direction, represents the aerodynamic parameters of the right side surface through the grid cell (α, β, λ); represents the aerodynamic parameters of the left side of the grid cell (α, β, λ); represents the aerodynamic parameters of the upper side of the grid cell (α, β, λ); represents the aerodynamic parameters of the lower side through the grid cell (α, β, λ); represents the aerodynamic parameters passing through the outer surface of the grid unit (α, β, λ); represents the aerodynamic parameters passing through the inner surface of the grid unit (α, β, λ); φ α,β,λ represents the aerodynamic parameters passing through the grid unit (α, β, λ); U 1,α,β,λ represents the air flow velocity through the grid unit (α, β, λ) along the x direction, U 2,α,β,λ represents the air flow velocity through the grid unit (α, β, λ) along the y direction; U 3,α,β,λ represents the air flow velocity through the grid unit (α, β, λ) along the z direction; s right represents the area of ​​the right side of the grid unit (α, β, λ), s left represents the area of ​​the left side of the grid cell (α, β, λ), s upper represents the area of ​​the upper side of the grid unit (α, β, λ); s below represents the area of ​​the lower side of the grid unit (α, β, λ); s out represents the area of ​​the outer surface of the grid unit (α, β, λ); s inner represents the area of ​​the inner side of the grid unit (α, β, λ). -q in formula (14) φ ·VB φV is the value corresponding to the grid unit (α, β, λ), where B φ It is obtained from A, where A is the leaf area density of the leaf cluster within the grid unit (α, β, λ).

[0191] Since the entire flow field space covers all grid cells, formula (14) can be converted into the following form:

[0192]

[0193] in represents the corresponding set of variables in each grid cell in Table 2, including the speed U i , U j , U k , turbulent kinetic energy κ and turbulent dissipation ε. ζ and Then are the adjacent coefficients derived from equation (14): α = 1, 2, 3 ... l, β = 1, 2, 3 ... n, γ = 1, 2, 3 ... m.

[0194] Subsequently, the Jacobi iteration method, an iterative algorithm for solving complex systems represented by a system of linear equations, is used to reduce computational complexity. ζ can be decomposed into the sum of the diagonal matrix D, the lower triangular matrix B, and the upper triangular matrix J, i.e., ζ = D + B + J. Then, the mathematical derivation is as follows:

[0195]

[0196] The weighted Jacobi iteration method is used, where the weight is ω, and formula (16) is iteratively propagated as follows:

[0197]

[0198] in, represents the set at time t+1. Represents the set at time t. Since there is wind speed initially, the changes in aerodynamic parameters at the corresponding time step can be obtained through repeated iterations.

[0199] Among them, ω = 2 / 3 is usually the selected parameter. According to the relationship ζ-D = +B + J, equation (17) can be expressed as follows:

[0200]

[0201] The converged solution of formula (18) can be determined as the spectral radius of the iteration matrix is ​​less than 1, that is, ζ(I-ωD -1 ζ)<1. Finally, the inlet wind speed boundary condition described in Section 2.4 is used as the initial value, and a series of gradually improved approximate solutions are generated through an iterative calculation process to extract the aerodynamic parameters related to each grid cell in the canopy and its surrounding space.

[0202] Fluent software provides a variety of turbulence model calculations, including the Spallanzani-Alemaras, k-ω, and k-ε model families. However, there is currently no forest-specific fluid model on the Fluent platform that contains the additional term B shown in Equation 2 and Table 2. φ Therefore, after elaborating the above fluid dynamics model, we conducted secondary development by integrating the custom tree model with the user-defined function (UDF) module in Fluent and writing the additional term B of the forest-specific k-ε model in Equation 2 using C language. φ Definition of .

[0203] 2.6. Digital twin structure of experimental wind-shaken tree:

[0204] Digital twin technology can accurately construct a digital copy that represents the research object (such as a tree) at a specific moment or in its life cycle, and make corresponding adjustments based on the changes in phenotypic characteristics caused by the natural environment, cultivation methods or genetic mutations and the aerodynamic parameters extracted in chronological order. To this end, we used a self-developed leaf segmentation algorithm and leaf extraction method to collect point cloud data through a ground laser scanner. In addition, a virtual entity is constructed based on the calculated wind correlation coefficient (i.e., aerodynamic-related parameters) to respond to a model corresponding to an organ of real size, thereby intuitively displaying changes in phenotypic characteristics. Its advantage is that it provides richer information transmission than traditional videos and images. In addition, a graphics engine, including related library functions and supporting frameworks, is required for video support for movie scene rendering, plant element collision response, and animation production. Therefore, we selected MeshLab middleware software from the Italian ISTI-CNR Institute and Speedtree software from IDV in South Carolina, USA as a customized flexible integration platform, combined with tree models and element properties extracted from our algorithm, for the description of tree-wind interaction behavior in bio-interaction applications, system cycles of various degrees of leaf clustering, vibration frequencies of wooden structures, and collisions of tree crowns under wind damage, and other mechanical phenomena.

[0205] 3. Experimental results:

[0206] The three trees were subjected to artificially generated wind speeds of 15, 30, and 45 m / s. In the following sections, we documented in detail the effects of wind on the morphological changes of the tree trunks and leaves, the vibration stresses on the wood structure, and the effects of air flow and turbulence as it passes through the leaves.

[0207] 3.1 Morphology and characteristics of stems and leaf clusters under wind stress:

[0208] Figure 3The displacement of the experimental trees is shown in the wind speed range from 0 to 45m / s. Each row presents the photos, point cloud data and initial crown segmentation results of the experimental trees under different wind speed conditions. As the wind speed increases, the trunks tilt downwind, and due to the strong wind, the leaf clusters are more closely gathered together.

[0209] Figure 3 In the figure, the four rows from top to bottom show experimental trees subjected to continuous strong winds at different wind speeds (0m / s, 15m / s, 30m / s and 45m / s). (a1), (b1), (c1) and (d1) show experimental photos taken when trees were subjected to wind mechanical loads under different wind speed conditions. (a2), (b2), (c2) and (d2) show the point cloud data of the relatively stable compression state reflected by the trees at a certain moment during the wind action. (a3), (b3), (c3) and (d3) show the segmented point clouds of each leaf cluster that constitutes the crown, and are annotated with a semi-transparent convex hull with consistent color coding.

[0210] Table 3 shows the characteristics of the leaf clusters as wind speeds increased. Compared to the records of stationary trees (Table 1), Table 3 shows fluctuations in the height of the three target trees. This may be due to the straightening of the more elastic tissue in the upper stem and the pushing of taller branches on the windward side into a vertical direction, resulting in an apparent increase in tree height (e.g., see Figure 3 The canopy width also shows changes, with the overall trend pointing to a compact canopy profile. This trend is reflected in the overall reduction of the canopy width perpendicular to the wind direction (north-south direction), while the lateral dimension is slightly stretched, resulting in an increase in the canopy width in the air flow direction (east-west direction) (see Figure 3 (c1) and (d1) in the figure).

[0211] Due to conflicting processes, the number of detected points for leaf clusters varies with wind speed. At relatively low wind speeds (15m / s), the movement of leaves and branches creates noise, resulting in more noticeable leaf point artifacts in the scan. As wind speed increases, overlap and clustering between adjacent leaves enhances the occlusion effect, while canopy compression results in a reduction in the scan target area. This results in an overall reduction in the number of points per leaf cluster.

[0212] The wind-blown canopy rearranges the internal plant components and simplifies the structure, thus facilitating air flow through the canopy. Plant flexibility and its natural characteristics promote local aggregation and staggered branches in pursuit of larger gaps between canopies to reduce wind loads. Therefore, this spatial concentration in plant material with uniform characteristics and higher consistency leads to a decreasing trend in the number of points per leaf cluster and the standard deviation of volume as wind speed increases, as shown in Table 3. Finally, due to the compact volume of each leaf cluster, its leaf area density increases.

[0213] Table 3. Adaptive changes of morphological characteristics and forest parameters of three tree species to different wind speeds (15m / s, 30m / s and 45m / s) at the scale of the entire crown and leaf cluster.

[0214]

[0215] Note: LAD stands for leaf area density. The above values ​​can be calculated from point cloud data.

[0216] 3.2. Vibration monitoring of experimental trees under different wind excitations:

[0217] Figure 4 (a) shows the different colored strapped inertial motion sensors installed on each tree trunk. Figure 4 (b) shows the extracted point clouds of the trunks of three trees (pink represents leaves, black represents trunks, and green represents leaves on the main trunk). Figure 4 (c)-(k) in the figure show the deviation of the horizontal vibration displacement of the tree trunk at different heights from the original vertical center position in two time periods. Each time period (80 seconds) is divided into two stages. The first stage (1-40 seconds) is the startup stage, when the capacitor gradually accelerates the motor rotation, and the fan blades gradually accelerate, resulting in a gradual increase in wind speed and tree movement. In the second stage (41-80 seconds), the power outage causes the fan blades to continue to slow down until they stop, which in turn causes the tree movement to slow down, accompanied by attenuated vibration to achieve energy consumption.

[0218] Figure 4 Figures (c), (d), and (e) show the obvious and characteristic movement of the tree trunks in two cycles under the condition of wind speed of 15 m / s. Tree 1, which is closest to the fan, has a relatively high leaf area density and an inverted conical crown, forming a compact leaf cluster close to the trunk, which contains many clustered leaves. Figure 4 As shown in (e) in Figure 1, tree 1 exhibits oscillatory behavior in the amplitude range of 2 to 4 cm. At the same time, due to the increased traction resistance caused by soil compaction, the swaying amplitude becomes less obvious as it approaches the lower part of the trunk. Figure 4 This is reflected in the red curves of different depths in (e).

[0219] The swing amplitudes of trees 2 and 3 were 2-6 cm and 2-10 cm respectively, which was mainly due to the pyramidal crown and more leaf clusters ( Figure 4 In addition, the interaction between tree 1 and the wind causes its shape to change continuously, and the wind field also changes accordingly, which increases the swing amplitude of the trees behind it. Figure 4(f)-(k) show similar vibrations of three trees with larger swing amplitudes under higher wind speed loading (the average swing amplitude of the experimental trees increased by 2-6 cm at 30 m / s, 2-10 cm at 45 m / s, and 15 m / s as a reference). It is worth noting that the bifurcated crown of Tree 3 is prone to structural imbalance because of its weak structural stability and is prone to collapse or breakage under strong wind conditions. In contrast, Tree 1 has a simple but rich canopy structure, reflecting its enhanced ability to withstand strong winds.

[0220] Figure 4 Overview of the experimental setup showing inertial motion sensors fixed at different heights on the trunk to detect trunk vibrations. Figure 4 (a) and (b) in Figure 3 show actual experimental scenarios of three target trees. Figure 4 (c)-(e) show the vibration detection curves of tree trunks at different heights along the main wind direction when the wind speed is 15m / s, including two stages of fan power increase and decrease, a total of four stages. Figure 4 (f)-(h) and (i)-(k) in the figure present the corresponding graphs when the wind speed is 30m / s and 45m / s, respectively.

[0221] 3.3. Extract aerodynamic parameters in the flow field using computer simulation:

[0222] In the wind tunnel experiment mentioned in Section 2.5, we calculated the aerodynamic model of the tree until the aerodynamic parameters were converged. To ensure stability, we adopted two strategies: adjusting the mesh and removing the mesh with singularities through singular perturbations to improve the fit of the model. The obtained results (see Figure 5 ) shows how a horizontal wind blowing from right to left affects trees under three different aerodynamic conditions, including wind speed profile, air pressure profile, and turbulent kinetic energy (TKE) profile. Figure 5 (a1), (b1) and (c1) show the tree models constructed according to different wind speeds.

[0223] and Figure 5 Figures (a2), (b2), and (c2) show the instantaneous wind speed distribution between tree crowns as wind speed increases. The results show that in the initial stage (15 m / s), the canopy gaps are relatively small, resulting in significant air blockage and the formation of a large dark blue area behind tree 3 (low air velocity of 0-5 m / s). However, as the leaves are compressed and the wind speed increases, the gap expansion is activated, which triggers the morphological changes of tree 3 and the occurrence of strong wind penetration, while a limited weak airflow area appears on the left. In addition, the handheld anemometer verifies the accuracy of our computer simulation program in calculating wind speed within the canopy.

[0224] Figure 5 Figures (a3), (b3) and (c3) show the pressure gradient caused by the deformation of leaf clusters under different wind speed conditions under asymmetric stress. A high-pressure area (300-600Pa) is formed on the right side of tree 1 due to direct exposure to the fan outlet, while a negative pressure area (-600 to -300Pa) appears on the left side of tree 3 due to the difference in wind speed and air flow dispersion. As the wind speed increases, the areas on the right side of the high-pressure area and the left side of the low-pressure area expand in the process of reducing the wind speed passing through each porous leaf cluster, and the pressure difference (286Pa, 653Pa and 1149Pa at 15m / s, 30m / s and 45m / s, respectively) gradually increases.

[0225] Figure 5 (a4), (b4) and (c4) in the figure show the turbulent kinetic energy (TKE) between tree crowns. Turbulent kinetic energy (TKE) is turbulent kinetic energy, which is an important indicator of the intensity and disorder of air flow in turbulence. Since TKE determines the magnitude and direction of wind force acting on leaves and branches, Figure 5 (a4), (b4) and (c4) in the figure show that the TKE is zero before the wind reaches the canopy, indicating that there is no turbulence in the free airflow. When the wind reaches the canopy edge of Tree 1, the minimum average TKE values ​​of Tree 1 at three different wind speeds are 0.4, 1.1 and 2.3, respectively, indicating that turbulence is generated by the collision and rebound between the tree and the airflow. As the wind enters the subsequent canopies, the average TKE increases (3.1, 6.4 and 9.7 m / s at 15 m / s, 30 m / s and 45 m / s, respectively). 2 / s 2 ), forming a wake that decreases with the crown, while dense leaf clusters divide it into multiple small TKE areas. This phenomenon also shows that the crown converts large-scale flows into small-scale vortices through the fragmentation of plant element groups, thereby inducing turbulence. In addition, higher turbulent kinetic energy values ​​were observed at some abnormal structures, accompanied by greater wind speeds. On the other hand, the turbulent kinetic energy at the corners of some leaf clusters is close to zero, reflecting that there is less airflow disturbance or interference at specific local locations. In addition, the distribution of turbulent kinetic energy and its size also change with the change of inlet wind speed. Higher inlet wind speed intensifies the degree of liquid separation and reattachment between irregularly shaped crowns, that is, it enhances the liquid recovery and mixing effects, as shown on the left side of Tree 3, where vortex phenomena appear.

[0226] Figure 5 A cross section of the interaction between wind and trees is shown, revealing the influence of the shape of the virtual tree model on the aerodynamic parameters in the simulated wind field along the wind direction (first row). The three columns from left to right present the wind sweeping of three types of trees at speeds of 15, 30 and 45 m / s. In addition, the second to fourth rows show the wind speed, air pressure and turbulent kinetic energy around the forest.

[0227] 3.4. Results of digital twin modeling:

[0228] Figure 6 A carefully constructed digital twin was presented that accurately replicates the phenotypic characteristics of the target tree and the changes in wind-induced forces under different aerodynamic parameters. The authors restored the tilt angle and azimuth of the leaves in each tree crown model and achieved physical-virtual integration by defining mechanical properties for the photosynthetic and non-photosynthetic components of the model tree. For example, elastic moduli were assigned to green wood (5.4-8.9 GPa) and branches (2.2-4.6 GPa). In addition, the mass of fresh leaves (0.57-0.81 g) and aging leaves (1.1-1.7 g) was determined. The elastic-gravity number of the rubber tree leaf was calculated (0.11-0.58), and the petiole bending stiffness value reflecting its flexibility was established (2.6-12.3 N·m 2 ). These values ​​define a biomechanical model of flexible pinnate leaves, elastically deformed petioles, and branches bending under buckling stresses, in pursuit of realism and accuracy in pattern recovery. A partial enlargement of our model is delivered in Figure 6 shown.

[0229] Our digital twin model demonstrates the various behaviors of trees in wind conditions with different aerodynamic parameters. At wind speeds of 15 m / s, the trees sway gently, almost gracefully following the wind direction. As wind speeds increase to 30 or 45 m / s, the trees begin to bend, especially in the direction of the strongest winds, with significant leaf and branch swaying and shaking, a phenomenon that stems from the trees' own physical resilience. Bending is a natural mechanism that allows trees to minimize their frontal area, thereby reducing severe damage, while their elastic vibrations help to more evenly distribute the forces exerted on them. The extent of these behaviors is limited by wind speed, duration, structural properties of the tree, and soil type and conditions. Overall, digital twins provide an effective way to recover the complex interactions between trees and wind without expensive and time-consuming physical prototypes.

[0230] Figure 6 Represents highly realistic 3D tree models generated for different wind excitation conditions (from top to bottom, the three lines represent wind speeds of 15m / s, 30m / s, and 45m / s, respectively). The left column shows the captured images of the experimental scene; the middle column shows the virtual tree digital twin model we built, whose organ depiction is consistent with the actual size and is set with a background environment; the right column shows a partial enlarged view of the digital twin tree model, including triangular mesh division and detailed depiction of tree organs under different wind threats.

[0231] 4. Discussion:

[0232] Rubber tree ecosystem functions in coastal areas of China are often affected by natural disturbances, such as hurricanes or frost damage; therefore, it is necessary to phenotype rubber tree responses to dynamic wind stress in the field. To pioneer the phenotyping of rubber trees under wind conditions, this study developed a self-built hardware system and analysis method that integrated industrial fans with laser scanners as a platform to obtain field observation data. Based on these observations, we proposed a hybrid computer graphics and forest-specific fluid simulation model that can display the aerodynamic parameters of rubber trees and characterize their phenotypic characteristics under unstable wind stress conditions. The method operates at the leaf cluster level and simulates wind disturbance events iteratively. Our custom fans simulated intermittent wind speeds under different climate conditions, and comprehensively demonstrated the effects of various wind speeds on experimental trees on the experimental platform. In addition, we used digital twin technology to reproduce short-term high-speed wind events and pauses to recover dynamic damping and tree-to-tree collision scenarios. In addition, our method provides a flexible experimental design and a scalable hardware system for any given input tree list to accommodate wind farms of adjustable width. Realism can be further enhanced by introducing more detailed variations, including variable trunk mass, resistance, crown front area, or skeleton structure.

[0233] Several studies have investigated the structural and morphological changes of trees under the effective wind resistance caused by different wind speeds. Regarding tree growth parameters, long-term exposure to high wind speeds leads to a decrease in stem height and leaf area, an increase in DBH, and wider growth rings on the leeward side, while a decrease in the stem length-to-radius ratio due to wind swaying. Regarding wind-induced tree motion, a literature suggests that conifers have smoother shapes and more flexible branches, which dissipate energy more efficiently than broadleaved trees; in addition, small branches exhibit high-frequency motion at lower wind speeds, which is accompanied by larger stem displacements under high wind speed conditions. For crown deformation, its components undergo elastic reconstruction when wind flows through it. Higher porosity allows the crown to reconfigure more easily, thereby reducing the drag caused by wind and making the crown more streamlined. In our approach, advanced terrestrial laser scanning technology can quickly record the stable state of the crown under different mechanical stress conditions and capture it in the form of point clouds. In addition, we designed computer graphics algorithms to reveal the phenotypic characteristics of trees under different wind speeds. Our Table 3 shows that when the wind speed doubles, the crown volume decreases by about 17%, while the leaf area density increases by about 12%. The experimental results also show that the crown volume decreases with increasing distance from the industrial fan, which is consistent with earlier studies that trees in the edge area, which are subject to direct resistance, are significantly more compressed than those in the sheltered forest interior to minimize the kinetic energy transfer effect. In addition, we found that the case of Tree 3, which is farthest from the industrial fan and is partially blocked by the infiltration airflow of Tree 1, has a larger swing amplitude than Tree 1, which is close to the fan and has a compact vase-shaped crown structure. These findings are consistent with previous studies that trees with conical crowns similar to Tree 1 experience less resistance than trees with cylindrical or bifurcated crowns, and that the upper disturbance motion can be significantly reduced by pruning or reducing the crown volume under all wind speed conditions. Indeed, a complex physical process is formed due to multiple factors, such as the material properties of leaves and branches, the mechanical response of the topological structure, the momentum transfer during the collision, the effect of various vegetation types on wind resistance, and the soil composition and its anchoring capacity. Therefore, this paper focuses on methodological innovation and aims at delivering basic results to promote future scientific research covering various complex and interrelated factors.

[0234] In order to create an accurate digital twin of forest wind damage, it is necessary to obtain high-resolution spatial data of forest structure, meteorological data (including wind direction), and modeling capabilities required to simulate wind-forest interactions. This study is the first attempt to combine graphics engineers with extracted tree phenotypic characteristics to build a basic digital twin to show the unsafe interaction between wind and trees, including detailed depiction of trunk vibration behavior, branch bending movements, leaf cluster aggregation, and morphological plastic deformation, which is of great significance for understanding the adaptability of rubber trees and their resistance to typhoon disturbances. In the future, the development of digital twins of forest wind damage can be applied to a variety of purposes: risk assessment, used to identify areas vulnerable to wind damage and prioritize management measures; mitigation planning, aimed at designing forest management strategies to enhance resilience to wind damage; and post-disaster analysis, assessing the impact of wind disasters on forests and identifying damage patterns to guide restoration efforts and future planning.

[0235] 5. Conclusion:

[0236] In this study, we proposed an innovative wind damage phenotyping platform and phenotyping method for rubber trees, combining point cloud analysis with digital twin technology. Our approach overcomes some limitations of existing methods and provides a practical solution to address this challenge. Key points include: First, our custom-designed industrial fan system is able to reproduce variable airflow conditions, such as gusty and intermittent winds, to accurately simulate the pressure exerted by real typhoons on trees. Second, our platform models the crown based on leaf cluster features in the point cloud and achieves fine-grained description based on phenotypic features to capture the dynamic morphological changes of the crown under different wind speeds. In addition, by considering aerodynamics and mechanical principles in a coordinated manner, we use the forest-specific k-ε turbulence equation, combined with the constructed tree shape model and numerical solution to calculate the wind characteristic parameters within the crown, and verify some of the calculation results through field measurements. Finally, we established a digital twin model that combines the measured tree phenotypic characteristics with the obtained plant mechanical properties to depict the actions of the target tree under various wind forces, thereby providing comprehensive and intuitive information for understanding its biophysical processes under wind-induced effects.

[0237] In summary, our framework effectively quantifies the unstable behavior of rubber trees by elucidating the complex interaction between wind and tree shape. Although this study is limited to small-scale experiments on specific tree species, its detailed characterization of tree organs and in-depth exploration of technical details of wind-induced damage are noteworthy. This approach and its in-depth analysis constitute the important highlights of this paper. To further promote this research, we recommend expanding the scope of experiments to include a variety of tree species, densities, and configurations as an important tool for sustainable forest management and stress-resistant cultivation during future storms.

[0238] The protection scope of the present invention includes but is not limited to the above embodiments. The protection scope of the present invention shall be based on the claims. Any replacement, deformation, and improvement of the technology that can be easily thought of by technicians in this field shall fall within the protection scope of the present invention.

Claims

1. A method for phenotyping trees under strong wind stress, characterized in that: include: Step 1: Build a field phenotyping platform, which includes a fan and target trees. The target trees are arranged in front of the fan. Start the fan to maintain a certain wind speed. Step 2: The laser scanner collects point cloud data of the target tree under a certain wind speed; Step 3: Model the canopy of the target tree based on the leaf clusters in the point cloud data, obtain the canopy graphic model, and calculate the phenotypic characteristic parameters of each leaf cluster; Step 4: Establish a k-ε turbulence model for forest environment; Step 5, calculate the wind characteristic parameters in the tree canopy through the k-ε turbulence model, the constructed tree canopy graphic model and the numerical solution; Step 6: Establish a digital twin model based on the wind characteristic parameters within the crown and the phenotypic characteristic parameters of each leaf cluster to describe the wind-induced response of the tree under a certain wind speed.

2. The method for phenotyping trees under strong wind stress according to claim 1, characterized in that: In step 1, the fan includes an industrial forced ventilation fan, which is connected to a lifting truss via a lifting device, and the lifting device is used to adjust the height of the industrial forced ventilation fan. There are 3 target trees in total, and the 3 target trees are arranged in a straight line in front of the fan.

3. The method for phenotyping trees under strong wind stress according to claim 1, characterized in that: The step 3 is specifically as follows: The crown is divided into leaf clusters according to the affiliation of the first-level branches, and the point cloud data corresponding to each first-level branch is marked with different colors. The crown layer is then modeled based on leaf clusters to obtain a crown graphic model, in which each leaf cluster corresponds to a first-level branch. The crown graphic model contains phenotypic characteristic parameters extracted from the point cloud, including the volume and leaf area density LAD of each leaf cluster.

4. The method for phenotyping trees under strong wind stress according to claim 1, characterized in that: The step 4 is specifically as follows: A k-ε turbulence model for forest environment is established, and a unified simplified mathematical expression of the k-ε turbulence model is obtained: in Represents the air flow velocity vector in three-dimensional space, that is U i represents the air flow velocity along the x direction; U j represents the air flow velocity along the y direction; U k represents the air flow velocity along the z direction; φ represents the relevant aerodynamic parameters, including U i , U j , U k , turbulent kinetic energy κ and turbulent dissipation ε; q φ represents the energy source term; B φ Indicates an additional term; The momentum equation, turbulent kinetic energy κ equation and turbulent dissipation ε equation are derived based on the unified simplified mathematical expressions; The momentum equation is: Where, μ represents the air viscosity coefficient, μ t represents the turbulent kinematic viscosity, ρ represents the air density, C d represents the drag coefficient of the forest; A represents the leaf area density of the leaf cluster; κ represents the turbulent kinetic energy, P represents the fluid pressure in the calculation space; |U| represents the absolute value of the total velocity of the three velocity components; and All of them conform to the Einstein summation convention; Among them, U1 is equivalent to the air flow velocity along the x direction, U2 is equivalent to the air flow velocity along the y direction; U3 is equivalent to the air flow velocity along the z direction; The turbulent kinetic energy κ equation is: Among them, β d =4, indicating the amount of energy loss when interacting with obstacles; σ κ =1.0, indicating the turbulent Prandtl number; ε indicates the turbulent dissipation; It represents the deformation rate or equivalent strain rate and conforms to the Einstein summation convention; The equation for turbulent dissipation ε is: Among them, p ε =1.31, c ε1 =1.44, c ε2 =1.92,a p =1.7,a d =1.7, β p =1; 5. The method for phenotyping trees under strong wind stress according to claim 1, characterized in that: The step 5 is specifically as follows: The air flow area containing the tree crown graphical model is meshed, and the U in each grid cell in the tree crown is calculated by numerical solution combined with formula (4), formula (7) and formula (11). i , U j , U k , turbulent kinetic energy κ, turbulent dissipation ε, and the pressure difference between two grid cells.

6. The method for phenotyping trees under strong wind stress according to claim 5, characterized in that: The step 5 is specifically as follows: The unified simplified mathematical expression, that is, the integral expression of formula (2) is as follows: Where V represents the volume of each grid cell; Convert formula (12) into the following form: represents the outward pointing unit normal vector relative to the surface of each grid cell; S represents the sum of the areas of the six faces on each grid cell; After discretizing formula (13), the flow field is converted into a grid structure, which is specifically expressed as follows: where δ x represents the distance between the centers of adjacent grid cells along the x direction, δ y represents the distance between the centers of adjacent grid cells along the y direction, δ z represents the distance between the centers of adjacent grid cells along the z direction, represents the aerodynamic parameters of the right side surface through the grid cell (α, β, λ); represents the aerodynamic parameters of the left side of the grid cell (α, β, λ); represents the aerodynamic parameters of the upper side of the grid cell (α, β, λ); represents the aerodynamic parameters of the lower side through the grid cell (α, β, λ); represents the aerodynamic parameters passing through the outer surface of the grid cell (α, β, λ); represents the aerodynamic parameters passing through the inner surface of the grid unit (α, β, λ); φ α,β,λ represents the aerodynamic parameters passing through the grid unit (α, β, λ); U 1,α,β,λ represents the air flow velocity through the grid unit (α, β, λ) along the x direction, U 2,α,β,λ represents the air flow velocity through the grid unit (α, β, λ) along the y direction; U 3,α,β,λ represents the air flow velocity through the grid unit (α, β, λ) along the z direction; s right represents the area of ​​the right side of the grid unit (α, β, λ), s left represents the area of ​​the left side of the grid cell (α, β, λ), s upper represents the area of ​​the upper side of the grid unit (α, β, λ); s below represents the area of ​​the lower side of the grid unit (α, β, λ); s out represents the area of ​​the outer surface of the grid unit (α, β, λ); s inner represents the area of ​​the inner side of the grid unit (α, β, λ); Convert formula (14) into the following form: where ζ is the adjacent coefficient derived from equation (14), Represents the corresponding variable set in each grid unit, and the corresponding variables include U i , U j , U k , turbulent kinetic energy κ and turbulent dissipation ε; Decompose ζ into the sum of the diagonal matrix D, the lower triangular matrix B, and the upper triangular matrix J, that is, ζ = D + B + J. The mathematical derivation is as follows: The weighted Jacobi iteration method is used, where the weight is ω, and formula (16) is iteratively propagated as follows: Where ω = 2 / 3, according to the relationship ζ-D = +B + J, equation (17) is expressed as follows: The converged solution of formula (18) can be determined as the spectral radius of the iteration matrix is ​​less than 1, that is, ζ(I-ωD -1 ζ)<1; Finally, the inlet wind speed of the air flow area after meshing is used as the initial value, and a series of gradually improved approximate solutions are generated through iterative calculation to extract the U in each grid cell in the tree crown. i , U j , U k , turbulent kinetic energy κ, turbulent dissipation ε, and the pressure difference between two grid cells.

Images