Simulation for Fitting the Bending Shape of Fruit Branches of Lycium barbarum Based on the Finite Element Method

Abstract

The accurate modeling of wolfberry plant morphology is the basis for theoretical and simulation analyses of the wolfberry picking process. The curved shape of the fruit branches makes it challenging to model Lycium barbarum (wolfberry) plants. This paper establishes a three-dimensional model of the branches under no gravity through field measurements, and then assesses the morphology of the branches under gravity load, fruit load, and branch load using finite element simulation. An orthogonal rotation combination experiment determined the relationship between branch morphology, length, growth angle, and growth mode parameters. The p-values of the prediction model were 0.0001, 0.0067, and 0.0203, respectively. Finally, the bending shape of the actual branches was verified against the branches generated by the prediction model. The experimental results show that the prediction model accurately models the fruit-bearing branches of Lycium barbarum. This paper introduces a method to quickly predict the bending shape of fruit-bearing branches of Lycium barbarum, providing a theoretical basis for rapid modeling of the L. barbarum plant and a simulation analysis for its harvesting.

1. Introduction

The interaction mechanism between the Lycium barbarum plant and the picker is the basis for the research and development of the L. barbarum picker. The interaction between the Chinese wolfberry plant and the harvester is at the center of the research and development of all kinds of harvesting equipment. To determine the key operating parameters of L. barbarum fruit abscission, it is necessary to analyze the kinematics of the L. barbarum plant and the picking parts. Due to the irregular structure of plants, it is difficult to obtain analytical solutions by theoretical analysis. Numerical simulations analyze the motion and stress responses of different parts of plants under picking and have been widely used in the study of fruit abscission [1,2]. Therefore, an accurate three-dimensional plant model is a requirement of fruit abscission simulation analysis [3].

Since the 1960s, researchers have developed various methods to create three-dimensional models of vegetables, field crops, and trees [4,5]. These mainly include automata, L-systems and particle systems, and image, data, and sketch reconstructions [6,7,8,9,10,11,12,13,14]. In 1968, the American biologist Aristid proposed the L-system [6,7]. In 1971, the Japanese scholar Honda created a computer simulation of the tree branch structure for the first time. In 1983, Reeves put forward the concept of a particle system [8,9,10,11]. In 1988, deReffye and the French Agricultural International Cooperation Research and Development Center (CIRAD) developed reference axis technology, known also as the automaton model [13,14]. The L-system, particle system, and automata create models according to the growth rules of plants. They are generally used to study the growth and development of plants, including their branching model, canopy distribution, and tendency. Image reconstruction is used to reconstruct the three-dimensional model of the plant by drawing on a single picture, multiple pictures, or videos with certain constraints [15]. Data reconstruction obtains the characteristic size or point cloud information of plants employing ground lidar, remote sensing, and three-dimensional scanning, and then divides different plants and organs using data segmentation to restore the three-dimensional shape of plants [16]. Image and data reconstruction create three-dimensional models of plants by using information obtained from natural plants. Both methods efficiently obtain a large amount of plant morphological information, but they are restricted by plant occlusion and data noise. Sketch reconstruction involves generating the corresponding three-dimensional plant model by drawing a sketch. One can flexibly modify the geometric size and topology of the plant on the mobile device to obtain realistic plant effects. This method can be used for the rapid implementation of various virtual scenes, such as games and gardens. L-system and automata are generally used to construct fruit harvesting and wind load response simulations and data reconstructions in combination with size measurements and constraints. In the fruit harvesting simulation process, different branch shapes and mechanical properties have significant differences in response to vibration, impact, and other picking actions. The bending shape of branches is one of the difficulties in plant modeling [17,18,19,20,21,22].

At present, there are three main methods used to construct the bending shape of branches. The first is to use the parameter curve to fit the measured data. This method obtains the approximate bending shape of branches, but the measurement workload is large, and the data are not universal [23]. The second method uses small vector segments to generate branches and model the bending of branches through the direction change of vector segments. The advantage of this method is that the bending direction of branches is very flexible; however, many parameters need to be given to control the bending of branches, and the direction of nonlinear change means that the branches are not smooth enough [24]. The third method uses the mechanical and finite element methods to determine the bending shape of branches through stress analysis, obtaining a more accurate branch shape [25,26]. The method of vector line segment may distort the shape of the whole branch because the direction vectors of each segment are set artificially. The method of parameter measurement and mechanical model calculation can obtain more realistic branches, but the measurement and model calculation work is large, the drawing period is long, and the measurement and calculation for specific branches lacks universality.

This paper aims to provide a fast and accurate modeling method for the bending shape of L. barbarum fruit branches. Firstly, the prediction model of branch shape is obtained by finite element analysis, and then the relationship between the parameters of the prediction model and branch length, growth angle, and planting mode is determined by a simulation test. Finally, the prediction model is verified with reference to real branches.

2. Materials and Methods

2.1. Classification of Branches

Figure 1 shows the branches of the Lycium barbarum at all levels, which was taken in Zhongning, Ningxia in August 2021. The Lycium barbarum in the picture is the main cultivated variety in Ningxia, the seven-year-aged Ningqi 7. The frame of the Lycium barbarum tree is composed of a trunk, primary branches, and secondary branches, with fruit-hanging branches growing on the primary branches and secondary branches. The trunk grows vertically. The primary and secondary branches are thick and grow along the initial growth direction, and the drooping deformation is small under the load of the lower branches and leaves. The fruit-bearing branches are slender and soft, and experience significant deformation under the load of fruit and leaves. Understanding this deformation is the research object of this paper.

L. barbarum fruit branches can be divided into independent branches without embranchment and compound branches with embranchment. These two types of branches grow on the vertical and horizontal planes, respectively. They can therefore be subdivided into four types of branches, named type I, type II, type III, and type IV. As shown in Figure 2 (by the red lines), type I and type III branches have no embranchments, type II branches have embranchments, and their superior branches grow along the vertical plane; type IV branches have embranchments, and their superior branches grow in the horizontal direction (approximately).

2.2. Parameters Measurement

When a large number of Lycium barbarum matured in June, three Lycium barbarum trees with well-developed crowns were selected, and the length, base diameter, top diameter, and growth angle of fruiting branches were measured. June is still the peak period for the growth of fruit branches, and fruit branches can grow about 10 mm every day. In order to ensure the reliability of branch size distribution, the data of three Lycium barbarum plants were measured in 2 days. The number of fruiting branches, the number of nodes on the fruiting branches, and the number of leaves and fruits growing on the nodes were counted, while the quality of ripe fruit, green fruit, flowers, leaves, and fruit stalks was also measured.

2.2.1. Length of Branches

The fruiting branches of the Lycium barbarum can be divided into current fruiting branches and biennial fruiting branches. The fruiting branches of the same year develop and grow in spring and autumn, and the length of the fruiting branches which germinated this spring was close to 400 mm when the fruit matured in the middle of June, and the length of the fruiting branches which germinated in August was about 300 mm when the autumn fruit was harvested in the middle of September. The length of some biennial branches with a large fruiting capacity can reach 800 mm, while the length of other biennial fruiting branches can be kept between 300 mm–500 mm by pruning.

2.2.2. Base Circle and Top Circle of Branches

The base diameter of the branch refers to the diameter at the beginning of the branch, and the top diameter refers to the diameter at the end of the branch. With the increase in the age of the Lycium barbarum, the annual growth of the branch diameter decreases year by year. Annual growth is further affected by precipitation, fertilization, and other factors [27]. There was little difference between the top diameter and the base diameter of the current-year-bearing branches that were measured for this study, but there was a large difference between the base diameter and the top diameter of the biennial branches due to repeated pruning and growth.

2.3. Fruits and Leaves on Fruiting Branches

Every 5 cm on the fruiting branch, there is a node where leaves grow. Ripe fruit, unripe fruit, and flowers will grow at the nodes during the full production stage. The quality of fruits and leaves growing on different lengths of fruit branches is different, which makes the bending shape of fruit branches very different [28]. We randomly selected 100 branch nodes, and counted the growth numbers of leaves, ripe fruits, unripe fruits, and flowers on each node.

2.4. Simulation Models

2.4.1. Geometric Model and Mesh Generation

The schematic diagram of the wolfberry tree without gravity is shown in Figure 3a. In a separate geometric model, part of the trunk and hanging fruit branches were cut off to reduce the amount of calculation and facilitate simulation analysis. A geometric model of a fruiting branch without gravity was created in Solidworks 2013 (Dassault Systems Simulia Corp., Concord, MA, USA). Four parameters were used to create the three-dimensional branches model: branch length, base circle diameter, circle diameter, and growth angle. The median length of the fruit-bearing branch was 550 mm, and the diameter of the base and top circles were determined using Formula (1). The branch growth angle was taken as 60°, and the three-dimensional model of the branch under no gravity was obtained, as shown in Figure 3b. The file was imported into Ansys Workbench 18.1 (ANSYS, Inc., Pittsburgh, PA, USA), which discretizes the 3D model through meshing, and creates the finite element model of branches. Because the branches have cylindrical structures, to better adapt to the arc of the cylinder, the tetrahedral mesh was used to divide the branches. In order to give consideration to calculation accuracy and efficiency, the mesh size of the branch was set at 5 mm, and the mesh size was set at 1 mm at the junction of branches and supports and at the junction of branch and fruits, as shown in Figure 3c. The model contains 16,458 elements and 33,366 nodes.

2.4.2. Simulation Parameters

The parameters used in the finite element model of the branches include material density, elastic modulus, and Poisson’s ratio. The density of the material was measured using the drainage method (Standard: ASTMD854-10). In the measurement, the hanging fruit branches of the Lycium barbarum were cut into small segments that were 20 mm long, and five segments were selected from the beginning to the end of the branches. Toluene solution was used as the measuring liquid in order to reduce volume measurement error [1]. We placed 50 mL of toluene solution into a 100 mL measuring cylinder (type: 1601; measurement range: 0–100 mL; precision: 1 mL; manufactured by Tianjin Glass Instrument Factory, Tianjin, China) and dipped the sample into the toluene solution to obtain the volume. An electronic balance (type: Professional digital jewelry scale 8028-series; measurement range: 0–100 g; precision: 0.001 g; manufactured by Shenzhen Diheng Electronics Co., Ltd., Shenzhen, China) was used to measure the change in mass before and after the branch was put in, and the mass difference was calculated to obtain the mass of the branch. The ratio of mass to volume was calculated to obtain the branch density. The elastic modulus of fruit-bearing branches was determined by the stress–strain method [1]. The formula for calculating the elastic modulus of the branch is as follows:

$$E=\frac{Fl}{A{l}_0}$$

(1)

E is the elastic modulus of the branch in the formula. A is the cross-sectional area of the branch. The diameter of the sample is measured using a Vernier caliper and calculated using the area formula of the circle. l is the length of the branch sample, which is measured using a Vernier caliper. F is the magnitude of the external force acting on the branch, and l₀ is the change in the length of the branch under the force F, which is obtained using the loading test of the universal testing machine. By substituting each piece of data into the formula, the elastic modulus of the hanging fruit branch can be calculated as shown in

Table 1. Simulation parameters.

	Density ρ/kg·m−3	Elastic Modulus E/MPa	Poisson Ratio υ
Fruiting branches	421	585	0.3

. Zhang Zui measured the elastic modulus of branches in the range of 420.96–749.28 MPa by using the similar method in [29]. The average elastic modulus of all branches of Lycium barbarum given in [30] is 665 MPa. This shows that the calculation results of the elastic modulus in this study are reasonable.

The Poisson’s ratio for wolfberry branches has been measured or used by many people. Hu Mingming set the Poisson’s ratio of goji berry branches at 0.3 during his modal analysis of goji berry plants in 2018 [30]. The average Poisson’s ratio of branches determined by Zhang Zui in 2016 was 0.327 [29]. When He Miao studied the dynamics of Lycium barbarum plants in 2018, she set the Poisson’s ratio of Lycium barbarum branches at 0.3 in the simulation calculation [31]. In 2021, Zhao Jian conducted a simulation test of wolfberry fruits falling off by force, and the Poisson’s ratio of hanging fruit branches was 0.3 [1]. Therefore, Poisson’s ratio of branches of the Lycium barbarum was assumed to be 0.3 as shown in Table 1.

2.4.3. Physical Model

L. barbarum branches are subject to gravity and the tensile force of fruit and leaves at the nodes. According to the growth and force directions of branches, branches are mainly subject to shear force and therefore belong to the beam element. The most common three-dimensional beam element models used in ANSYS static analysis are the beam4 constant section and beam187 variable section beams. The diameter of the L. barbarum branches decreases from the base circle to the top circle, meaning they should be categorized as variable section beams. Therefore, the beam187 variable section beam element is selected. The whole branch is regarded as a Euler beam, with six degrees of freedom per solving element, i.e., three linear degrees of freedom and three rotational degrees of freedom. The deflection curve differential equation of Euler beam under load [32] is:

$$\frac{{\mathrm{y}}^{″}}{{\left[1+{({\mathrm{y}}^{\prime })}^2\right]}^{\frac{3}{2}}}=-\frac{M(x)}{EI}$$

(2)

In this formula, EI is the bending stiffness and M(x) is the bending moment at the section. This is a nonlinear equation, and the exact solution is challenging to obtain. Many scholars have simplified this formula for different applications and to improve its practicability [32]. However, the obtained differential equations are still challenging to apply to the direct calculation results to drive the three-dimensional modeling of branches. In this study, the finite element solution was used to obtain the deformation shape of branches.

2.4.4. Constraints and Loads

The root of the superior branch of the fruit-hanging branch is subject to fixed constraints, and the superior branch plays a supporting role with the fruiting branch. The load on the branch includes its gravity, the tension caused by the gravity of the fruit and leaves at each node, and the pressure of the lower branches. The gravity of the branch itself is applied by inertial load, and it can be calculated according to Formula (3):

$$G_1=\rho Vg=\frac{1}{12}\rho g\pi l\left(D^2+d^2+Dd\right)$$

(3)

where G₁ is the branch gravity, ρ is the branch density, and g is the gravitational acceleration. Using the circular platform volume formula and Formula (1), the relationship between the gravity of the branch itself and the branch length can be obtained as follows:

$$G_1=\frac{1}{12}\rho g\pi \left(1.2\ast {10}^{-5}\ast l^3+0.0224l^2+11.2l\right)$$

(4)

The maximum ratio of the number of flowers to the number of Chinese wolfberry fruits could be calculated in the first ten days of June. At this time, the root fruit of the branch is the first to mature, while the fruit that matures at the beginning of July is in a state of flowering, and the ratio of flowers to fruit is close to one. Considering that there are two ripe fruits, four green fruits, six flowers and three leaves on the node, it is calculated that the flower mass accounts for 6.7% of the node mass, and that this ratio decreases rapidly with the increase in the number of mature fruits. Therefore, the flower quality is ignored in the simulation. The number and quality of growing leaves, mature fruits, and immature fruits at each node were averaged to determine the concentrated mass force applied on each node as follows:

$$G_2=n_1{m}_1+n_2{m}_2+\left(n_1+n_2\right)m_3$$

(5)

where G₂ is the equivalent tension at the node, n₁ is the number of ripe fruits at the node, n₂ is the number of unripe fruits at the node, m₁ is the mass of ripe fruits, m₂ is the mass of unripe fruits, and m₃ is the mass of the fruit stalk. If the fruit-hanging branch has lower branches at the node, the force of the branches on the fruiting branch was applied in the same way as the inertial load, and the load was calculated by Formula (3). Figure 4 shows the constraint and load application positions. In the figure, A is the fixed constraint, and the rest represent at each stage the concentrated mass force applied at the respective node.

2.5. Sequential Loading Test

In the experiment, the weight was loaded according to fruit maturity at each stage of the branch to simulate the bending process of the branch in its natural state. Therefore, a total of five loads were performed. First, only the gravity of the branch was loaded, that is, the inertial load marked by K in Figure 4. Second, based on the K load, three concentrated mass forces were applied on nodes B, C, and D, at the root of the branch. E, F, and G loads were then applied. The third load was based on the above loads. The fourth load continued to increase the three concentrated mass forces of H, I, and J. For the fifth load, the branch load was increased at C, E, and G.

Before solving, ‘path’ was added in WORKBENCH to the branch centerline to obtain the coordinate values of the branch centerline before and after deformation to determine the branch shape curve.

2.6. Orthogonal Simulation Test

Through intuitive understanding, it can be found that the overall shape of fruit-hanging branches is a part of quadratic or higher-order curves. Determining the parameters describing the curve, that is, the coefficients of polynomial items, is the key to accurately describing the shape of the branches. It was assumed that the highest order term of the polynomial describing the shape of fruit-bearing branches is less than cubic; that is, the prediction polynomial coefficient of the shape curve is 3–4, and the expression of the shape curve is assumed as follows:

$$y=a_1{x}^3+a_2{x}^2+a_{3}x+a_4$$

(6)

where, x and y are the horizontal and vertical coordinates of each point on the centerline of the branch under the branch coordinate system, defined in Figure 3a. The branch coordinate system takes the center of the base circle of the branch as the coordinate origin, parallel to the ground and along the branch direction as the positive direction of the x-axis, and downward perpendicular to the ground as the positive direction of the y-axis. The coefficients of the morphological curve are a₁, a₂, a₃, and a₄ (a₁ might be 0). They are all functions that have branch length, growth angle, and implantation position as independent variables.

To investigate the relationship between branch length, growth angle, planting position, and each coefficient, three-factor and three-level orthogonal simulation experiments were carried out. The levels of branch length (A) and growth angle (B) are determined to be 300 mm, 550 mm, and 800 mm, and 40°, 60°, and 80°, respectively, according to the measured size distribution. The fruit-bearing branch attachment position (C) codes are 1, 2, 3, and 4. The Box–Behnken method of orthogonal tests establishes the efficiency and effectiveness of the test.

Table 2. Orthogonal test arrangement.

Test No.	Factor A Length/mm	Factor B Growth Angle/°	Factor C Growth Location
1	300	60	4
2	550	60	2
3	550	80	4
4	550	40	4
5	800	60	4
6	300	40	2
7	550	40	1
8	550	80	1
9	300	60	1
10	800	60	1
11	550	60	2
12	800	80	3
13	550	60	3
14	550	60	3
15	300	80	3
16	550	60	3
17	800	40	2

shows the tests arrangement.

2.7. Comparison and Verification between Branch Shape Prediction Curve and Real Curve

Because the branches are soft, shaking, and sheltered from each other in the full productive stage, it is difficult to ensure the accuracy of the abscissa and ordinate of each point of the branches by directly measuring the real branch shape. Therefore, the image recognition method described in reference [33] was used to determine the coordinates of each point on the branch centerline. The length, growth angle, and growth position of branches were measures using field pictures. Branch edges and non-target areas were manually identified. The true-color image was transformed into a gray image in Matlab R2010a (MathWorks.Inc., Natick, MA, USA), and the gray image was transformed into a binary image using fixed threshold segmentation. The find function was used to derive the coordinates of non-zero pixel values of the picture, that is, the coordinate values of all pixels on the branch. The pixel coordinates of each point were transformed into length coordinates using the image scale, and then the coordinates were transformed from the picture coordinate system into the branch coordinate system to obtain the horizontal and vertical coordinates of 50 equidistant points on the branch centerline.

The length, growth angle, and growth position serial number of the above branches were brought into the morphological prediction curve equations to obtain the prediction curve equation results. Fifty abscissae equidistant in the x direction and were brought into the prediction curve to obtain the coordinate values of 50 points on the prediction curve.

The abscissa and ordinate values of 50 equidistant points on the centerline of the real branch were compared with the abscissa and ordinate values of the 50 points on the shape prediction curve, and the deviation was used to evaluate the accuracy of the prediction curve.

3. Results

3.1. Measurement Results of Fruit Branch Parameters

Figure 5 shows the length distribution of hanging fruit branches, from which it can be seen that 90% of the hanging fruit branches are distributed between 330 and 760 mm. SPSS 22.0 (International Business Machines Corporation, Armonk, New York, NY, USA) was used to conduct a Shapiro–Wilk (S-W) test on the measured results, and the p-value of S-W test was 0.765 (>0.05). It can be seen that the length distribution of branches meets normal distribution.

The measured data show that the top diameter of the fruiting branches is 47.06% to 61.76% of the base diameter. Figure 6 shows the size distribution of the base and top diameters of fruit branches.

The linear relationship between base diameter and branch length is obtained using the least-squares fitting technique. It can be seen from Figure 6 that with the increase in branch length, the fitting line between the beginning diameter and the end diameter is at a certain angle, indicating that the fruit branch of Lycium barbarum has a cone angle which varies with the branch length in a small range. In order to further reduce the model variables, the functions of taper angle and branch length are introduced to represent the end diameter. Branch length and Formula (1) can therefore be used to calculate the base diameter and top diameter of the branches as follows:

$$\{\begin{array}{l}D=0.0032l+2.6151\\ \theta =-0.000005\ast l+0.0082\\ d=\theta l\end{array}$$

(7)

where, D represents the diameter of the base circle of the branch, d represents the diameter of the top circle of the branch, and l represents the length of the fruiting branch, while θ represents the tangent of the taper angle.

Growth angle α represents the included angle between the branch growth direction and the superior branch growth direction, as shown in Figure 2a,c. The growth angle range of fruiting branches is shown in

Table 3. Length, growth angle, and diameter of branches.

Type of Branches	Length/mm		Diameter at the Base Circle/mm		Diameter at the Top Circle/mm		Growth Angle/°	Number of Branches
	Range	Annual Growth	Range	Annual Growth	Range	Annual Growth
Fruiting branches	310–1000	314	1.7–3.4	1.7	0.8–2.1	1.7	40–87	45–97

, and the data were obtained through field measurement.

The quality of 100 randomly selected leaves, ripe fruits, unripe fruits, flowers, and fruit stalks was measured, as shown in

Table 4. The quantity and quality of ripe fruits, leaves, green fruits, flowers, and fruit stalks growing on the nodes.

	Mass/g				Nodes
	Maximum Value	Minimum Value	Average Value	Standard Deviation	Maximum Value	Minimum Value	Average Value	Standard Deviation
Ripe fruit	1.240	0.320	0.812	0.192	6	2	3.21	0.217
Leaves	0.155	0.074	0.121	0.023	4	2	2.57	0.132
Unripe fruit	0.346	0.056	0.157	0.065	4	0	1.95	0.046
Flowers	0.110	0.020	0.055	0.019	6	0	3.41	0.017
Stalks	1.069	0.867	0.978	0.048	7	2	4.17	0.219

.

3.2. Deformation of Branches under Different Loads

The deformation process of fruit-bearing L. barbarum branches under different loads is shown in Figure 7. Figure 7a shows that the branch grows upward along the growth angle, with only a small deformation. Figure 7b–d show the deformation process of the fruit-bearing branches from the beginning of the fruit ripening to the full production stage. In this process, the bending degree of branches gradually increases, and the growth direction changes from oblique upward growth to ground growth. There are two main reasons for this. On the one hand, due to the long force arm at the end of the branch, the load increases, and the torque increases as the fruit at the end of the branch matures. On the other hand, the end of the branch is growing, the branch is thin, and the degree of lignification is low, making the end-bending effect more pronounced. In Figure 7e, the branch load was added. When the length of fruit-bearing branches exceeds 300 mm, the thorns at some nodes will continue to grow and develop into the next level of fruit-bearing branches. These new branches are composed of newly germinating leaves, and green fruits will be born at the node. The branch length in the simulation is 500 mm, and two lower-branch loads are applied. The application of a lower-branch load makes the branch reach the final shape of bending.

Using ANSYS post-processing, the coordinate values of 50 points on the branch centerline before and after deformation are obtained, and the shape curve of the branch before and after deformation is drawn, as shown in Figure 7f. The figure shows that for a 500 mm branch without load, under a load of fruit, leaves, and branches, the branch becomes slenderer, the growth direction changes from oblique upward growth to oblique downward growth, and the vertical height of the end of the branch reaches 441 mm.

Table 5. Result of simulation.

Loads		Branch Gravity	3 Nodes	6 Nodes	9 Nodes	Lower Branch
x/mm	Before deformation	433	433	433	433	433
	After deformation	436	445	491	601	687
	Deformation	3	9	46	110	86
	Cumulative deformation	3	12	58	168	254
	Percentage	1.04%	3.70%	18.26%	43.21%	33.80%
Y/mm	Before deformation	−249	−249	−249	−249	−249
	After deformation	−244	−228	−147	44	193
	Deformation	5	16	81	191	149
	Cumulative deformation	5	21	101	292	441
	Percentage	1.04%	3.70%	18.25%	43.22%	33.79%

provides the contribution rates of different loads to branch deformation. The contribution rates of tip node and branch loads were the largest, with 43.22% and 33.79%, respectively.

Quadratic or cubic curves can describe Figure 7f, the curved shape of fruit-bearing branches, and the fitting equation of the branch shape is as follows:

$$y_1=0.9781x^2+0.0743x-0.0132(R^2=0.9946)$$

(8)

$$y_2=-2.9208x^3+2.8752x^2-0.2509x-0.0021(R^2=0.9999)$$

(9)

Equation (8) is a quadratic fitting curve, and Equation (9) is a cubic fitting curve. The R² values show that both methods can fit the branch shape. To determine the relationship between the parameters of the morphological equation and the fruit-bearing branch length, growth angle, and growth position, it is necessary to analyze the Box–Behnken test further.

3.3. Relationship between Fruiting Branch Parameters and Branch Shape Equation Parameters

Table 6. Simulation results of Box–Behnken test.

Test No.	Cubic Curve					Conic Curve
	a1	a2	a3	a4	R2	b1	b2	b3	R2
1	−3.5772	1.6055	−1.1306	−0.0005	1.0000	0.4525	−1.0346	−0.0019	0.9998
2	−2.635	2.6223	−0.2378	0.0034	0.9998	1.275	−0.0566	−0.0014	0.997
3	−1.3218	0.8902	−0.5767	−0.0001	1.0000	0.342	−0.5168	−0.0014	0.9997
4	−2.7624	2.5061	−0.0279	0.0043	1.0000	1.1941	0.1371	0.0001	0.9987
5	−1.2636	1.9301	−1.1905	−0.0003	1.0000	0.9566	−0.9959	−0.0078	0.998
6	−0.8666	1.3506	−0.6116	0.0002	0.9995	0.5403	−0.4139	−0.0092	0.9625
7	−0.9497	1.5631	−0.5054	−0.0018	0.9994	0.6195	−0.2645	−0.0129	0.9715
8	−0.9497	1.5631	−0.5054	−0.0018	0.9994	0.6195	−0.2645	−0.0129	0.9715
9	−0.9497	1.5631	−0.5054	−0.0018	0.9994	0.6195	−0.2645	−0.0129	0.9715
10	−0.0131	0.033	−0.4192	−0.0002	1.0000	0.0229	−0.4172	−0.0003	1.0000
11	−0.0131	0.0329	−0.4192	−0.0002	1.0000	0.0229	−0.4172	−0.0003	1.0000
12	−0.4578	1.1095	−0.2222	0.0012	0.9999	0.7086	−0.1299	−0.003	0.9986
13	−0.1345	0.258	−0.1391	−0.0002	0.9999	0.1461	−0.1145	−0.0013	0.9950
14	−0.5552	1.2695	−0.9245	−0.0088	0.9978	0.4701	−0.6284	−0.0292	0.9718
15	−0.3598	1.0092	−0.5665	−0.0004	1.0000	0.5062	−0.3862	−0.0122	0.9774
16	−0.5337	1.3275	−0.1867	0.0017	1.0000	0.5502	0.11	−0.0207	0.9972
17	−0.3285	1.0178	−0.113	0.0033	0.9999	0.6042	0.0242	0.0058	0.9991

shows that both morphological curves predict the shape of branches. However, the R² value of the cubic curve is closer to 1, so the prediction effect is better. If the parameters describing the shape curve of fruit-bearing branches can be determined by the parameters of branch length, growth angle, and growth position, rather than through simulation calculation every time, the drawing efficiency of fruit-bearing branches will be greatly improved. The relationship between each factor and each parameter of the cubic curve is analyzed visually and by variance. The analysis results are shown in

Table 7. Results of range and variance analysis of shape parameters.

Factors	A				B				C
Morphological Parameters	a1	a2	a3	a4	a1	a2	a3	a4	a1	a2	a3	a4
Range (R)	2.1298	0.8612	0.0543	0.0028	0.7073	0.3607	0.8388	0.0045	0.617	0.6154	0.4694	0.0047
p-value	0.005	0.3016	0.9808	0.5126	0.6823	0.8279	0.002	0.2183	0.9027	0.8009	0.2147	0.1371

. The range analysis shows that the branch length had the most significant influence on the parameters a₁ and a₂, with a range (R) of 2.1298 and 0.8612. The branch growth angle has the greatest influence on parameter a₃, while the branch growth position has the greatest influence on parameter a₄. Analysis of variance shows that branch length had a significant effect on a₁. The growth angle significantly affects parameter a₃.

Figure 8 shows the changing trend between each variable and the morphological parameters. With the increase in branch length from 300 mm to 800 mm, the shape parameter a₁ first increases and then decreases. The shape parameter a₂ reaches the maximum when the branch length is 300 mm, decreases gradually, and increases gradually when the branch length exceeds 600 mm. The coefficients a₃ and a₄ increased monotonically with the increase in branch length. The shape parameters a₁ and a₂ did not change significantly with the growth angle. The shape parameters a₃ and a₄ gradually increased with the growth angle by 38% and 50%, respectively.

Design-expert 8.0.6 (Stat-Ease, Inc., Minneapolis, MN, USA) directly calculates the regression model with each shape parameter as the response function and each factor as the independent variable. The p-values of the a₁, a₂, and a₃ prediction models are 0.0001, 0.0067, and 0.0203, respectively, indicating that the regression model is statistically significant. The value of parameter a₄ in each process is less than 0.01, so the average value of each process −0.00012 can be used in the final prediction model. The final expression of the prediction model is shown in Formula (10):

$$\{\begin{array}{l}{y}^3{=\mathrm{a}}_1{\mathrm{x}}^3{+\mathrm{a}}_2{\mathrm{x}}^2{+\mathrm{a}}_3\mathrm{x}-0.00012\\ {\mathrm{a}}_1=-16.92676+0.03218A+0.31355B-1.88578C-0.00012AB-{0.00002A}^2-{0.00194B}^2+{0.39538C}^2\\ {\mathrm{a}}_2=8.11901-0.01307A-0.06238B-1.41951C+0.02218BC+0.00001A^2\\ {\mathrm{a}}_3=-2.57623+4.81457A+0.00512B+0.11411C-0.000003A^2\\ {\mathrm{a}}_4=-0.00012\end{array}$$

(10)

3.4. Error Analysis of Branch Shape Prediction Model

Figure 9a–d show the front image and image processing of the branches. Figure 9a is the original picture. After manually removing the non-target area, Figure 9b is obtained. Figure 9b is transformed into a gray image, and threshold segmentation is carried out to obtain the binary in Figure 9c. From the enlarged picture, there are many burrs on the edge of the branch, which will increase the error in the subsequent coordinate picking. A morphological closed operation corrodes the branches to obtain a relatively smooth binary, as shown in Figure 9d. Figure 9e is the scatter diagram of the branches in the pixel coordinate system. To obtain the coordinate value of the branch centerline, it is necessary to take the average value of multiple y values corresponding to the x value of each pixel coordinate. However, the rotation on the left side of the branch overlaps the x coordinates. Therefore, the picture is divided into left and right at pixel abscissa 640, as shown by the red line in Figure 9e. The mean value of multiple x values corresponding to the y value is calculated in the left half of the picture to determine the actual coordinate values of 50 points equidistant from the centerlines of the left and right branches, respectively.

The branch length is 492.4 mm, the growth angle is 74.8°, and the growth position is 3. The branch shape curve parameters a₁, a₂, and a₃, are 0.000032, −0.00369, and −0.28750, respectively. The three parameters are incorporated into Formula (10) to calculate the predicted coordinate values of 50 points on the branch centerline. The scatter diagram of predicted coordinates and real coordinates in the length coordinate system is shown in Figure 9f.

Table 8. Real value, predicted value, and relative error of coordinates of each point on the middle line of fruit-bearing branch.

No.	x/mm	yt/mm	yp/mm	Relative Error	No.	x/mm	yt/mm	yp/mm	Relative Error
1	0.00	0.00	0.00	0.00%	26	149.95	−27.10	−18.27	8.03%
2	6.03	−1.50	−1.86	0.33%	27	155.98	−24.29	−13.27	10.02%
3	12.06	−5.67	−3.94	1.57%	28	162.01	−20.40	−7.45	11.78%
4	18.09	−10.11	−6.22	3.54%	29	168.04	−15.77	−0.76	13.65%
5	23.96	−13.14	−8.56	4.17%	30	174.07	−11.76	6.81	16.89%
6	29.99	−16.63	−11.08	5.05%	31	179.94	−7.61	15.10	20.66%
7	36.02	−19.97	−13.65	5.75%	32	185.97	2.48	24.60	20.12%
8	42.05	−22.84	−16.24	6.00%	33	188.09	5.30	28.18	20.81%
9	48.08	−24.87	−18.80	5.52%	34	189.88	7.40	31.31	21.75%
10	53.95	−26.21	−21.23	4.53%	35	191.51	9.41	34.23	22.58%
11	59.98	−28.75	−23.62	4.67%	36	194.61	12.71	40.01	24.83%
12	66.01	−30.69	−25.86	4.39%	37	196.02	17.28	42.72	23.14%
13	72.20	−32.78	−27.97	4.38%	38	198.02	20.54	46.70	23.79%
14	78.07	−34.02	−29.73	3.90%	39	199.94	24.61	50.63	23.67%
15	83.94	−35.63	−31.23	4.00%	40	202.19	27.55	55.37	25.30%
16	89.97	−36.86	−32.46	4.00%	41	203.99	30.15	59.28	26.50%
17	96.00	−36.80	−33.33	3.16%	42	205.96	35.69	63.69	25.47%
18	102.03	−37.46	−33.79	3.34%	43	208.00	37.98	68.36	27.63%
19	108.06	−38.05	−33.82	3.85%	44	209.20	40.75	71.17	27.67%
20	113.93	−38.59	−33.37	4.75%	45	209.93	44.66	72.92	25.70%
21	120.12	−38.71	−32.36	5.78%	46	210.96	48.57	75.39	24.40%
22	125.99	−38.28	−30.85	6.76%	47	212.02	52.81	77.99	22.90%
23	132.02	−35.91	−28.70	6.56%	48	212.99	54.76	80.39	23.31%
24	138.05	−33.46	−25.89	6.89%	49	214.01	59.00	82.94	21.78%
25	144.08	−30.70	−22.38	7.57%	50	215.31	71.23	86.27	13.68%

shows the actual value y_t, predicted value y_p, and relative error of the coordinates of 50 points on the branch centerline. Figure 9f shows that the prediction curve simulates the shape of branches growing upward first and then bending downward by gravity. Among the 50 points, 16 points had a relative error of less than 5% of the predicted value, 18 points were between 5–20%, and the relative error of 16 points was higher than 20%.

4. Discussion

The fruiting branches begin to sprout and grow in April each year. With branch length and diameter growth, the leaves and fruits on each node of the fruit-bearing branches gradually develop and grow. In this process, the branch load gradually increases, their lignification degree gradually increases, and their shape gradually bends, which is irreversible. Therefore, the branches reached the final bending shape at the full fruit stage with the maximum load. In mid-June, the fruit on the upper end of the fruiting branch develops and ripens first, followed by the fruit in the middle of the branch, and finally by the fruit at the end of the branch. This also explains why the fruit begins to mature from the root of L. barbarum branches. The lignification of fruit branches is a gradual process. Compared with the apical fruits, root fruits are less likely to maintain the stability of the L. barbarum plant system due to the shorter force arms and less load on the branch system when they mature. Although the load increases significantly when the terminal fruit is mature, the plant system is developed enough to maintain its stability with the gradual lignification of fruit-bearing branches. In order to study the damage to broad-leaved trees caused by wind load, the deformation model of trees under the action of gravity was also established in [34]. In that study, the tree data were first obtained by ground lidar, and the tree model was reconstructed. Then the anti-gravity load was applied to the trees to obtain the plant model without gravity. The study shows that the effect of gravity load on the deformation of branches is a cumulative effect. This is similar to the results of this paper.

In order to verify the accuracy of the model, it is necessary to establish a plane coordinate system for the actual branches which is exactly the same as the simulation calculation. However, the actual coordinates of each point on the axis of the actual branch are difficult to obtain by measurement. This is due to the fact that the shape of the actual branch is three-dimensional, and there is a sudden change in the diameter of the branch node. The branches have unsmooth bending to obtain more light. The error of points no. 41, 43, and 44 exceeds 25%, which may be caused by pruning or spatial competition among branches. Catherine Jirasek pointed out in [35] that, under natural conditions, the bending of branches is sagged by their load, and also affected by phototaxis, artificial pruning, and lateral wind, and therefore making calculations is very complex.

The modeling of Lycium barbarum was also studied in [1]. In order to ensure that more ripe fruit and less green fruit is obtained during the harvest process, simulation experiments and field experiments on the shedding sequence of fruit and fruit stalks from branches under tension were carried out. In [1], the fruit-pedicel-stalk-branch system was constructed by reverse engineering. The reverse engineering method constructs the three-dimensional model completely according to the real object, and the model accuracy is therefore high, which provides a reliable geometric model for simulation analysis. However, the reverse engineering method is completely based on physical creation, and the flexibility of the model is poor. Reverse engineering with the help of 3D scanning methods for reconstruction requires the creation of the complex branch structure of the whole plant (including thousands of leaves and dozens of hanging fruit branches), or even different branch types of multiple plants. The modeling process is not efficient enough to meet these needs. In this study, the parameters that determine the branch shape (length, growth angle, and growth mode) can obtain from the actual branches, or can be randomly determined by the research needs within a reasonable range. After these three parameters are determined, the branches in accordance with the growth law can be obtained by the branch shape prediction model. It should be pointed out that the branch morphology predicted in this study is one of the actual growth forms of branches. This is because the branches with the same length, growth angle and growth position will have some differences in growth morphology when affected by water stress, nutrient competition among branches, accumulated temperature during the growth period, and so on. In order to study the power transfer characteristics of Lycium barbarum plants, the three-dimensional model of Lycium barbarum plants was established by measurement, and the decreasing characteristics of power among all levels of branches were determined in [31]. In this study, because of the large number of hanging fruit branches and the difficulties of measuring the bending shape, the three-dimensional model does not include hanging fruit branches in order to facilitate the analysis, and the study does not include the dynamic transfer characteristics from hanging fruit branches to the Lycium barbarum fruit. The actual measurement can truly reflect the geometric shape of specific plants, but the complex bending shape of fruit branches and the large number of fruit branches bring difficulties to the measurement.

Since 1959 [36], researchers have tried to calculate the geometric shape and internal stress of branches under gravity and ecological constraints by means of mechanical analysis. Fournier obtained a semi-analytical solution [37], and Fourcaud developed a numerical finite element model [25,26] to adjust the branch shape. The above methods focus on establishing the mechanical equation and solving the equation to generate the geometric shape of the branch. However, after the branch parameters are given, they still need finite element calculation to obtain the branch shape. In a 2014 study [38], Guillaume reconstructed the biomechanical model of branches considering the effect of progressive maturity on branch morphology. These studies aim to explore the mechanism of branch morphological development and establish the mechanical model of plant morphology. In our study, the progressive, mature finite element calculation is considered to ensure the accuracy of the model prediction, and the prediction equation is determined by orthogonal simulation experiments to ensure the rapidity of modeling. Therefore, we can flexibly and quickly create a three-dimensional model of the hanging fruit branch.

The efficient establishment of a branch bending model is helpful to further study the harvesting mechanism of Chinese wolfberries and improve the development efficiency of Chinese wolfberry harvesters. It can be seen from references [1,31] that due to the limited accuracy and efficiency of branch bending modeling, some studies are carried out only on a single fruit branch, while others directly ignore the role of fruit branch, which reduces the generality and reliability of the research results. By carrying out simulation research on the basis of this study, we can quickly obtain a large number of shapes of fruit-hanging branches in accordance with mechanical law so as to quickly establish the three-dimensional model of fruit-hanging branches. In the analysis of vibration harvesting mechanisms, fruit branches with different lengths and branch angles can be created flexibly, and the finite element method can be used to solve the vibration response of branches. In the simulation analysis of airwave picking, we can create a single fruit branch or multiple fruit branches growing in various ways, so as to study the effect of different air flows. In the simulation research process of medlar pruning and pesticide application, the crown shape or application position of the rich fruit branches can be more intuitively viewed.

The next step of this study is to construct the whole plant model, and explore the dynamic response of leaves and fruits under the action of contact force and air flow. On the other hand, it is necessary to introduce more influencing factors in the construction of the plant model to improve the openness of the prediction model and incorporate different Chinese wolfberry varieties and berries to expand the application of the model according to their own cultivation characteristics.

5. Conclusions

The prediction model of fruit-bearing branch bending of L. barbarum was obtained through the finite element simulation test of branches. The experimental analysis shows that the prediction result of the cubic model is better than the quadratic model. The length of branches has a significant impact on the cubic and quadratic term coefficients of the prediction model, the growth angle of branches has a significant impact on the primary phase coefficient of the prediction model, and the constant term of the prediction model of different fruit-bearing branches has little change. The growth position of branches had no significant effect on the bending shape of branches.

The large deformation and sagging of fruit-bearing branches of L. barbarum is mainly caused by the loads of mature fruits and lower branches. Therefore, the longer the branch, the more nodes on the branch, the more fruits and branches growing on the branch, and the greater the bending deformation of the branch. The bending of fruit-bearing branches is an irreversible process of gradual development, and the branches reach the final shape in the full production stage. The fruit ripening sequence on the fruit-bearing branch is conducive to protecting the stability of the plant structure.

The branch shape prediction model better simulates the branch shape, but the coordinate error of the corresponding points is still large. Adding more factors, such as the phototaxis of branches, artificial pruning, and nutritional competition between branches, into the model to obtain a more accurate prediction model is the next step for this research.