## Abstract

Numerous natural and synthetic systems can be modeled as clusters of interacting cantilever beams. However, a closed-form mathematical model capable of representing the mechanics of multiple interacting cantilever beams undergoing large deflections has yet to be presented. In this work, a pioneering mathematical model of the force–deflection response of multiple, inline, interacting (i.e., contacting) cantilever beams is presented. The math model enables the determination of the force–deflection response of a system of interacting cantilever beams and is predicated upon the “Pseudo Rigid Body Model” concept. The model was validated through data triangulation experiments which included both physical and computational studies. An analysis of the mathematical model indicates it is most accurate with deflections less than 50 deg. In the future, the model may be used in high throughput phenotyping applications for investigating stalk lodging and estimating the flexural rigidity of crop stems. The model can also be used to gain intuition and aid in the design of synthetic systems composed of multiple cantilever beams.

## 1 Introduction

Many structures can be accurately modeled as a single cantilever beam. Consequently, the force–deflection response of cantilever beams has been well studied [1–6]. However, numerous natural and artificial systems are better represented as a group of mutually interacting cantilever beams. Equations for modeling systems of interacting cantilever beams have not been presented previously. For example, consider the case of a wheat field in which each individual wheat stem can be approximated as a cantilever beam. When subjected to external forces, each wheat stem will contact and interact with its neighbors. The same is true for many agricultural cropping systems. Other natural and synthetic systems with similar attributes include grasses, forests, hair, nanotube arrays, brooms, and brushes. A method capable of accounting for the interactions between adjacent cantilever beams in such systems (i.e., adjacent beams contacting one another) is needed.

For example, wind damage (i.e., windthrow) can have major economic impacts on crops and forests [7]. The problem of stalk lodging (crops breaking in windstorms prior to harvest) is particularly troublesome in staple grain crops such as rice, wheat, and corn [8,9]. Numerous modeling efforts have been conducted to generate understanding to alleviate the problems of stalk lodging and windthrow [10,11]. However, these modeling efforts often employ small deflection assumptions and often do not account for plants supporting one another in dense canopies [12]. When plant-to-plant interactions are accounted for less accurate empirical relationships are often utilized as opposed to mechanistic equations [13]. A mechanistic, large deflection model of interacting cantilever beams could be used to improve the accuracy of modeling efforts and generate greater understanding of stalk lodging and windthrow. Increased understanding in this area is needed to develop mitigation strategies to the problems of stalk lodging and windthrow.

Accurate models of large deflection, interacting cantilever beams are also needed to improve crop phenotyping efforts. For example, one method to alleviate the problem of stalk lodging is to selectively breed for stronger stalks. Several electromechanical devices have been developed to measure the bending stiffness and bending strength of plant stems to aid such breeding efforts [14,15]. Unfortunately, selective breeding studies require very large sample sizes and current mechanical measurement devices are relatively low throughput. In other words, they are generally unable to process enough plants per hour to be economically viable. The throughput of these devices could be improved if a methodology of accounting for plant-to-plant interactions were available. Currently, all neighboring plants must be removed from the field to eliminate plant-to-plant interactions when using these devices. The process of removing neighboring plants to eliminate plant-to-plant interactions is time-consuming and costly as it requires destroying a significant portion of the potential plants available for testing. However, these devices could be adapted to test multiple plants in a high throughput manner if mechanistic equations were available that could accurately represent the interactions occurring between the plants being tested. A high throughput device of this nature could significantly advance stalk lodging research.

Lastly, a mechanistic methodology for representing interacting cantilever beams could be used in the design of synthetic systems. It is apparent that the stiffness of a system of cantilever beams is a function of the flexural rigidity of the individual beams, the length of the beams, and the spacing between the beams. However, there is currently no closed-form solution capable of quantifying how the stiffness of the system changes as a function of spacing, length, and stiffness of individual beams.

This paper takes a first step toward addressing these problems by providing a closed-form solution to model the force and deflection response of multiple, inline, interacting cantilever beams undergoing large deflections. Consider Fig. 1 which depicts a row of vertical cantilever beams, of equal length, placed directly inline along the *x*-axis with equal spacing. A rigid body (hereafter referred to as a force bar), oriented parallel to the *z*-axis at a constant height *h*, moves in the *x*-direction applying a displacing force to each beam. As each beam deflects, it will contact the adjacent cantilever beam(s). Each beam will eventually be displaced until the *y*-coordinate of the end of each displaced beam is at height *h* at which point the beam will pass beneath the force bar. The purpose of this paper is to outline a mathematical model to solve the force–deflection response of the system depicted in Fig. 1. In particular, the solution to the simple case of inline, prismatic, rectangular beams that are equally spaced is presented. The solution is compared and validated against physical experiments and a corresponding finite element modeling simulation. More complex cases including non-prismatic beams with circular cross sections that may contain nonlinear material properties require additional considerations (e.g., tribology modeling) and remain the subject of future work.

## 2 Theory

Several approaches to solving the large deflection contact problem depicted in Fig. 1 are plausible. The authors chose to utilize the “pseudo rigid body model” (PRBM) approach [16,17]. The PRBM is an approximation method that provides an accurate and efficient manner to analyze large deflection problems. The PRBM method represents the force–deflection response of flexible members with a combination of rigid bodies, precisely placed pin joint(s) and torsional springs. In this study, the PRBM of a cantilever beam with a point force applied along its length is used to determine the force–deflection response of multiple inline, interacting cantilever beams as shown in Fig. 1. The derivation of the solution follows. Several concepts and terms from the single cantilever PRBM are used in the derivation. For the readers’ convenience, a brief outline of the single cantilever beam PRBM may be found in the Appendix. A complete derivation of the single cantilever beam PRBM can also be found in Refs. [16,17]. For simplicity, the same nomenclature and naming conventions are used in the derivation presented below as are used in Ref. [16].

### 2.1 Multiple Inline Interacting Cantilever Beam Model.

*l*are placed with equal spacing

*s*along the

*x*-axis. Let the force bar, oriented parallel to the z-axis, move across the beams, applying a force at a known, constant height

*h*that deflects all three beams. Weight and friction are considered negligible, and deflection is assumed to be within the linear elastic range. While this system is dynamic in nature, there will be a point at which the frontmost beam is at its maximum deflection (i.e., immediately prior to passing under the force bar and returning to its vertical position). The applied force from the force bar at this point is referred to as

*F*. The PRBM of this specific scenario is depicted in Fig. 2. In addition, free body diagrams of the PRBM of each individual beam are illustrated in Fig. 3. Note that for ease of visualization, the deflected angle of the PRBM will often be described with respect to the angle

_{peak}*β*, rather than Θ, where

*i*indicates a measurement for beam(

*i*).

*F*, the deflected angle of the PRBM (Θ) of each beam must be determined. Starting with the first beam

_{peak}*b*, the length of the fixed base of the PRBM (see Fig. 2), is expressed as

*θ*

_{o1}. Initially, however, the force bar is assumed to apply a horizontal force to yield values of

*n*

_{1}= 0 and

*γ*

_{1}= 0.8517 (see Fig. 2 and Appendix for term definitions) which when input into Eq. (2) provides an initial estimate for Θ

_{1}. The initial estimate for Θ

_{1}is then used to determine

*θ*

_{o1}and then

*ϕ*

_{1}(i.e.,

*ϕ*

_{1}is assumed to be perpendicular to

*θ*

_{o1}). The

*ϕ*

_{1}value is then used to obtain

*n*

_{1}, which is in turn used to find a new value for

*γ*

_{1}. The new value for

*γ*

_{1}is put into Eq. (2) and used to update the estimated value of

*Θ*

_{1}at which point the process repeats. The process is repeated until

*γ*

_{1}converges. The converged value for

*n*

_{1}is used to calculate

*K*

_{Θ}, which is in turn used to determine the first beam’s torsional spring constant

*K*. As all beams in the system are identical, they all possess the same spring constant

*K*. A diagram of this iteration process is provided in Fig. 4 for the readers’ convenience.

*x*-coordinate of the tip of the first beam, (

*x*

_{1}), at its maximum deflected state is

*i*= 1. As shown in Fig. 2, the horizontal distance a beam extends past the next beam’s origin is expressed as

*x*

_{i}>

*s*condition is not met, beam(

*i*) will not contact beam(

*i*+ 1) prior to the first beam passing below the force bar and returning to its vertical position. The deflected angle of all remaining beams(

*i*> 1) can be expressed as

*i*) exerts a force on beam(

*i*+ 1) and that beam(

*i*+ 1) exerts an equal but opposite force on beam(

*i*) (i.e., Newton’s Third Law [18]) supports the following notation

*F*,

_{i}_{i+1}is the force acting on beam(

*i*) induced by the beam (

*i*+ 1). The force vectors acting between beam(

*i)*and beam(

*i*+ 1) are assumed to act perpendicular to beam(

*i*)’s end angle (

*θ*). This assumption is discussed in more detail in Sec. 5.2 Limitations. For each beam in the system,

_{oi}*ϕ*is used to determine

_{i}*n*which is in turn used to update each beam’s

_{i}*γ*(see Eqs. (A4)–(A6) in the Appendix).

_{i}*Θ*

_{i}), we can now proceed to solve for

*F*. To do so, we start by solving for the force on the final beam in the system as it is the only beam with just one force and a torque acting on it. With a moment in the counterclockwise direction being defined as positive, the

_{peak}*x*and

*y*components of all force vectors are defined as

*F*

_{32}(the force on the final beam in the system) from the free body diagram shown in Fig. 3(c) yields

*final beam’s resistive force*can be expressed as

*F*(the force applied to the first beam). This is accomplished by summing the moments acting about the first beam’s torsional spring and rearranging to solve for

_{peak}*F*

_{peak}*EI*of the beams in the system can obtained, where

*E*is Young’s modulus and

*I*is the second moment of area. This requires solving Eq. (15) for

*K*. Dividing both sides of Eq. (15) by

*K*yields

*F*

_{21}

*/K*which is found in the denominator of the right-hand side of Eq. (17) is a known value that has been solved previously using Eq. (14). In other words, both sides of Eq. (14) can be divided by

*K*to yield an expression for

*F*. This expression can then be substituted into Eq. (17) to yield an explicit equation for

_{21}/K*K*. This explicit expression for

*K*is then substituted into the following equation to determine the flexural rigidity of the beams:

The earlier equations were derived from the perspective of the first beam’s maximum deflected state. This approach was taken for simplicity and to provide the context in the upcoming sections. However, the equations can be implemented as a function of the *x* position of the force bar. More specifically, virtually placing the force bar at the beam height *l* initially and incrementally decreasing the force bar height, until it is equal to *h* or until the desired deflection is reached, results in a continuous plot of the reaction force *F* on the force bar versus beam displacement.

*F*as the force bar moves across the beams results in the idealized plot shown in Fig. 5. Notice that the force will reach an initial peak or maximum force (

*F*) when the first beam has reached its maximum deflection (Fig. 6(a)). Immediately after the first beam deflects under the force bar,

_{peak}*F*will sharply drop and then continue to rise until the second beam reaches its maximum deflection (Fig. 6(b)) at which point

*F*=

*F*. This trend will continue until the force bar approaches the final beam in the system and the number of beams in contact with one another is reduced as shown in Figs. 6(c) and 6(d). At this point, the reaction force

_{peak}*F*will be less than the initial

*F*. The reaction force

_{peak}*F*will be equal to

*F*(

_{peak}*m*) number of times where

*t*is the total number of beams and

*u*is the maximum number of beams in contact at the frontmost beam’s maximum deflection.

Note the equations presented above can be used to solve for the reaction force given *EI*, the length of the beams (*l*), spacing distance (*s*), and the height of the force bar (*h*). Alternatively, if given the reaction force (*F*), the length of the beams (*l*), the spacing of the beams (*s*), and the height of the force bar (*h*), the equations can be used to solve for *EI*.

## 3 Methods to Assess the Accuracy of the Multiple Inline, Interacting Cantilever Beam Model

The closed-form Multiple Inline Interacting Cantilever Beam Model presented in Sec. 2.1 was implemented in a python script, and the results were compared to physical and computational experiments to confirm its accuracy. Across the experiments, all input parameters were varied to assess the robustness of the closed-form solution. Details of the physical experiment are presented first followed by a description of the computational finite element model experiments.

### 3.1 Physical Experimental Setup.

The closed-form solution method presented in Sec. 2.1 was validated through physical experiments. In particular, six rectangular sheet metal strips were securely placed directly inline and equally spaced by width *s*. A rectangular aluminum bar (i.e., a force bar) oriented perpendicular to the sheet metal strips at height *h* was attached to a load and displacement sensor. The force bar was then slowly driven across the beams, causing each beam to deflect and eventually pass under the bar. The force–displacement response of the bar was input into the closed-form solution to estimate the *EI* of the beams. Five different sets of beams underwent testing. The estimated *EI* values were then compared to the actual *EI* values for each set of beams. Actual *EI* values were determined from three-point bending experiments as described below.

The geometric and material properties of the beams were acquired using the following protocol. Beam lengths were measured by the same individual with an imperial ruler, while the width and thickness were measured with a set of digital calipers. Four randomly selected rectangular strips cut from the same sheet metal stock underwent individual three-point bending tests using an Instron universal testing machine (Model 5965, Instron Crop., Norwood, MA). The Instron software (Bluehill 3.0) was used for instrumentation control and data acquisition of displacement, force, and the calculation for *E*. All tests limited displacement to prevent yielding or other physical deformation. The supports were spaced 10 cm apart and the loading anvil was lowered at a rate of 0.13 cm s^{−1}. The test stopped at a beam displacement of 0.3 cm. The mean *E* of the four beams was then calculated and assigned to all corresponding beams cut from the same stock. A total of five sets of beams were created. Each set consisted of six beams, and each set was constructed to possess a unique *EI* value.

Each set of beams were subjected to three tests with different beam-to-beam spacings (*s*). The first test had a beam-to-beam spacing of 19.1 mm, the second with a spacing of 24.9 mm, and the third with a spacing of 49.9 mm. The height of the force bar was adjusted so that the beams would not yield but they would come into contact with at least one other beam during the test. Across the 15 experiments (5 sets of beams × 3 spacings), the height of the force bar ranged 5.08–16.51 mm below the top of the beams. Table 1 summarizes the testing conditions for all 15 experimental tests.

Set | EI (N mm^{2}) | s (mm) | l (mm) | h (mm) | # of Force Peaks |
---|---|---|---|---|---|

A_{1} | 41,900 | 24.9 | 180 | 167 | 4 |

B_{1} | 63,000 | 24.9 | 180 | 165 | 4 |

C_{1} | 111,000 | 24.9 | 231 | 219 | 4 |

D_{1} | 154,000 | 24.9 | 205 | 191 | 4 |

E_{1} | 335,000 | 24.9 | 231 | 219 | 4 |

A_{2} | 41,900 | 49.9 | 180 | 167 | 5 |

B_{2} | 63,000 | 49.9 | 180 | 165 | 5 |

C_{2} | 111,000 | 49.9 | 231 | 219 | 5 |

D_{2} | 154,000 | 49.9 | 205 | 191 | 5 |

E_{2} | 335,000 | 49.9 | 231 | 224 | 5 |

A_{3} | 41,900 | 19.1 | 186 | 175 | 3 |

B_{3} | 63,000 | 19.1 | 186 | 175 | 3 |

C_{3} | 111,000 | 19.1 | 237 | 225 | 3 |

D_{3} | 154,000 | 19.1 | 211 | 198 | 3 |

E_{3} | 335,000 | 19.1 | 237 | 225 | 3 |

Set | EI (N mm^{2}) | s (mm) | l (mm) | h (mm) | # of Force Peaks |
---|---|---|---|---|---|

A_{1} | 41,900 | 24.9 | 180 | 167 | 4 |

B_{1} | 63,000 | 24.9 | 180 | 165 | 4 |

C_{1} | 111,000 | 24.9 | 231 | 219 | 4 |

D_{1} | 154,000 | 24.9 | 205 | 191 | 4 |

E_{1} | 335,000 | 24.9 | 231 | 219 | 4 |

A_{2} | 41,900 | 49.9 | 180 | 167 | 5 |

B_{2} | 63,000 | 49.9 | 180 | 165 | 5 |

C_{2} | 111,000 | 49.9 | 231 | 219 | 5 |

D_{2} | 154,000 | 49.9 | 205 | 191 | 5 |

E_{2} | 335,000 | 49.9 | 231 | 224 | 5 |

A_{3} | 41,900 | 19.1 | 186 | 175 | 3 |

B_{3} | 63,000 | 19.1 | 186 | 175 | 3 |

C_{3} | 111,000 | 19.1 | 237 | 225 | 3 |

D_{3} | 154,000 | 19.1 | 211 | 198 | 3 |

E_{3} | 335,000 | 19.1 | 237 | 225 | 3 |

Note: Expected peaks refers to the number of force peaks that occurred during the experiment.

#### 3.1.1 Flexural Rigidity Estimations From Physical Experiments.

For each physical experiment, the reaction forces exerted on the force bar were measured with a load cell and the *x*-displacement of the force bar was recorded. These data were then used to determine the *EI* of the beams using a python script that follows the process outlined in Sec. 2.1. The input parameters to the script were the peak force (*F _{peak}*) from the experimental tests, beam-to-beam spacing (

*s*), beam length (

*l*), and height of the force bar (

*h*). The load cell utilized in the experimental setup only measured the horizontal component of the force exerted on the force bar (i.e., the

*x*component of the force vector). The total force exerted on the force bar (i.e., the

*x*and

*y*components of the force vector) was calculated by assuming the force vector acted perpendicular to the beam end angle,

*θ*. As expected, each experimental test experienced several peak forces due to sequential engagement and disengagement of the beams with the force bar. The peak force (

_{o}*F*) used to estimate

_{peak}*EI*was calculated by averaging the first

*m*number of expected peak forces (see Table 1) from each test according to Eq. (A1).

### 3.2 Computational Finite Element Model Setup and Validation.

The closed-form solution presented in Sec. 2.1 was further investigated by developing a parametric finite element model (FEM) simulation of multiple inline, interacting cantilever beams undergoing large deflections. The model could simulate a wide array of input parameters. The two-dimensional FEM was developed in abaqus/cae 2019 [19–21]. The metal beams were modeled as two-noded linear beam elements, and the force bar was modeled as a discrete rigid. All contact was modeled as frictionless. Material and contact damping was used to aid in force–displacement stability and was confirmed to have a negligible effect as compared to the internal potential energy of the beams. The model was analyzed as a dynamic simulation in Abaqus/Explicit 2019, capturing full contact and nonlinear effects [19–21]. The FEM formulation was validated against physical experiments conducted with one, two, three, and six beams, as shown in Fig. 7.

#### 3.2.1 Finite Element Method Experiments.

Once the FEM was validated, it was used to assess the accuracy of the closed-form solution over a broad range of input parameters. However, as the first step in this process and to create preliminary insight into the closed-form solution’s predictivity, a single in-depth analysis of the reaction force on the force bar was conducted. In particular, the case of six beams of *EI* = 63,000 N mm^{2}, *s* = 24.9 mm, and *l* = 180 mm being deflected by a force bar driven at *h* = 165 mm was analyzed via physical experimentation, via a FEM simulation, and with the closed-form solution. The entire force–displacement response of the system from each method (i.e., from the FEM, physical experiment, and closed-form solution) was then compared.

After the preliminary, in-depth single case experiment explained earlier, a more comprehensive data triangulation [22,23] experiment was conducted. Specifically, one researcher simulated FEMs of eight identical, inline beams (*l* = 180 mm) at 10 different *s* values undergoing deflection due to 10 unique *h* values. Thus, a total of 100 FEMs were produced (10 *s* values × 10 *h* values = 100). For all beams and simulations, *EI* was held constant at 63,000 N mm^{2}. A summary of the input parameters for these FEM simulations is provided in Table 2.

s (mm) | h (mm) | h/l | Max Θ (deg) |
---|---|---|---|

2.00 | 89.81 | 0.50 | 67.2 |

6.22 | 97.36 | 0.54 | 64.0 |

10.44 | 104.91 | 0.58 | 60.4 |

14.67 | 112.46 | 0.63 | 56.8 |

18.89 | 120.01 | 0.67 | 53.2 |

23.11 | 127.56 | 0.71 | 49.4 |

27.33 | 135.11 | 0.75 | 45.4 |

31.56 | 142.66 | 0.79 | 41.1 |

35.78 | 150.21 | 0.84 | 36.5 |

40.00 | 157.76 | 0.88 | 31.3 |

s (mm) | h (mm) | h/l | Max Θ (deg) |
---|---|---|---|

2.00 | 89.81 | 0.50 | 67.2 |

6.22 | 97.36 | 0.54 | 64.0 |

10.44 | 104.91 | 0.58 | 60.4 |

14.67 | 112.46 | 0.63 | 56.8 |

18.89 | 120.01 | 0.67 | 53.2 |

23.11 | 127.56 | 0.71 | 49.4 |

27.33 | 135.11 | 0.75 | 45.4 |

31.56 | 142.66 | 0.79 | 41.1 |

35.78 | 150.21 | 0.84 | 36.5 |

40.00 | 157.76 | 0.88 | 31.3 |

Note: Each simulation consisted of eight identical, inline beams with *l* = 180 mm and *EI* = 63,000 N mm^{2}.

The resulting 100 force response plots (10 *s* values × 10 *h* values) of the 100 FEM simulations and the corresponding system dimensions (*s, l*, and *h*) of each FEM were then given to a second independent researcher. This researcher input these values into the python script described previously (i.e., the closed-form solution) to determine *EI* of the beams in the system. The actual *EI* value used in the FEM was then compared to the EI value obtained via the python script.

### 3.3 Closed-Form Solution Sensitivity Analysis.

To determine how slight measurement errors in beam dimensions or spacing affect the closed-form solutions *EI* predictions, a sensitivity analysis of the closed-form solution was conducted. In particular, the closed-form solution was utilized to calculate *EI* while holding *l* and *F _{peak}* constant but evaluating all possible combinations of

*h*and

*s*with

*h*ranging from 0.5

*l*to 0.99

*l*and

*s*ranging from 0.005

*l*to 0.99

*l*. The numerical derivative of

*EI*was then assessed with respect to

*h/l*and

*s/l*. The derivative values indicate how small changes or errors in input dimensions influence

*EI*calculations.

## 4 Results

The results are presented in four sections. First, results from the comparison between physical experiments and the closed-form solution are provided. Second, an in-depth single case examination comparing the closed-form solution’s force–displacement response to its corresponding physical experiment and FEM simulation is presented. Third, results comparing the closed-form solution to 100 parametric FEMs with various input parameters are presented. Lastly, the results from the sensitivity analysis of the closed-form solution are presented.

### 4.1 Physical Experiments Versus Closed-Form Solution Results.

The closed-form solution showed good agreement with physical experiment data. Figure 8 shows a plot of actual *EI* values of sheet metal beams obtained via three-point bending test versus computed *EI* values of the same beams obtained from the Multiple Inline Interacting Cantilever beam experiments outlined in Secs. 3.1 and 3.2. The *R*^{2} value between the actual *EI* values and computed *EI* values from the closed-form solution was 0.99. However, the slope of the linear regression line was 0.945, indicating that the closed-form solution slightly underpredicts the actual *EI* values of the beams. Figure 8 shows the linear regression line between actual and computed *EI* values in blue as well as an ideal 45-degree line in red.

### 4.2 In-depth Single Case Examination and Force Response Comparison.

As described in Sec. 3.2.1 the force response of a set of six beams was obtained via physical experimentation, a FEM simulation, and with the closed-form solution. Figure 9 shows the force response (*F _{x}*) of the beams plotted against the

*x*-displacement of the force bar, with data from the physical experiment shown in blue, data from the FEM simulation shown in green, and data from the closed-form solution predictions shown in orange. The force response between the physical experiment, FEM, and closed-form solution prediction is in good agreement. However, the closed-form solution slightly overpredicts the force peaks. This is likely because the PRBM approach does not fully capture the complex kinematics of beams slipping under the force bar. A correction factor was developed to overcome this limitation. In particular, we found that multiplying

*h*by 1.0085 yields better agreement with the displacement and force peaks. The force response of the closed-form solution using this correction factor is shown in red in Fig. 9.

### 4.3 Computational Finite Element Model Versus Closed-Form Solution Results.

The closed-form solution was further investigated by comparing its predictions to those of 100 FEM simulations as described in Sec. 3.2.1. In particular, the peak force from the FEM simulations along with the corresponding *l, h,* and *s* values was used as inputs to the closed-form solution to calculate *EI.* The actual *EI* value from the FEM model was then compared to the calculated *EI* from the closed-form solution. When comparing the actual *EI* values from the FEM simulations to the closed-form solution *EI* predictions, a mean percent error of 4.19% with a standard deviation of 4.08% was observed. The percent error from each combination of spacing and maximum deflection is shown in Fig. 10. This figure indicates the closed-form solution becomes less accurate as deflections become larger (i.e., lower *h/l* values). More specifically, at larger deflections, the closed-form solution tends to overpredict *EI* when beams are spaced closely together. At deflections less than 50 deg, the mean percent error is 1.32% (SD = 2.61%).

### 4.4 Closed-Form Solution Sensitivity Analysis.

To determine how slight measurement errors in beam dimensions or spacing affect the closed-form solution’s *EI* predictions, a sensitivity analysis was conducted. The numerical derivative of *EI* with respect to *h/l* and *s/l* is shown in Figs. 11 and 12, respectively. Notice in Fig. 11 that the derivative of *EI* is small unless *h/l* is greater than about 0.9. Thus, errors in *h* and *l* measurements can produce larger errors in *EI* predictions at high *h/l* ratios. Figure 12 shows that the closed-form solution’s *EI* predictions are less sensitive to errors in *s/l* measurements than to errors in *h/l* measurements. Note that the light shading running through the middle of the plot is due to interactions between beams no longer occurring.

## 5 Discussion

Modeling the force–deflection response of multiple inline cantilever beams undergoing large deflections has application to numerous industries and is a classical mechanics problem. The authors are unaware of any other studies which have reported a closed-form solution method to this classical problem. The closed-form solution presented in this study provides an accurate representation of the system’s force–deflection response over a range of system configurations. Results indicated that the closed-form solution is most accurate for deflections less than 50 deg (which corresponds to an *h*/*l* ratio of approximately 0.7). At *h*/*l* ratios greater than 0.9, the solution can be sensitive to small experimental measurement errors in *h* or *l*. Therefore, the closed-form solution is most appropriate for analyzing systems in which the *h*/*l* ratio is between 0.7 and 0.9. This range is typical of numerous natural and synthetic systems. When analyzing systems undergoing extremely large deflections (i.e., greater than 50-deg deflections), the authors recommend using a more time-intensive solution method (e.g., nonlinear finite element analysis (FEA)).

The closed-form solution is most advantageous for applications that require rapid calculations. In particular, it is accurate while requiring substantially less time to evaluate than finite element methods. It took just over 2 min to run the 100 simulations outlined in Table 2 using the closed-form solution. Running the same simulations with the finite element software took more than 1 h. The computation times of both methods are dependent on numerous factors but across multiple experiments we found the closed-form solution to be at least an order of magnitude faster than finite element methods. In addition, large deflection finite element simulations capable of modeling interactions between beams require specialized software and expertise to create and evaluate. It took an experienced researcher approximately 1 h to create the 100 FEA models described in Table 2. Finally, the computational efficiency of the closed-form solution will enable it to be directly embedded in devices used to measure stalk lodging resistance. Electromechanical devices used to measure stalk lodging resistance have limited computational capabilities. Running complex finite element software that requires inverting very large matrices is impractical on these system architectures. However, such system architectures can easily solve the closed-form solution presented in this work.

Another advantage of the closed-form solution is its relative simplicity (i.e., lack of integrals, limits, large matrices, and derivatives). The simplicity of the solution enables it to be used in developing basic intuition regarding the complex mechanical response of systems of closely spaced cantilever beams undergoing large deflections. For example, the authors are currently developing a device to measure the flexural stiffness of agricultural plots of wheat stems. The device will be used to test hundreds of plots a day and to investigate genetic mechanisms underlying wheat stem stiffness and stalk lodging resistance. The authors are using the closed-form solution method to develop intuition regarding influential parameters and inform design decisions. For example, the effect of wheat stem height and spacing (planting density) on the force response of the system is easily understood by analyzing the closed-form solution method. The method may also be applied to study other natural systems or to aid in the design of synthetic systems. Finally, the solution presented in this work can serve as a stepping stone to developing more complex solution methods to systems that violate some of the assumptions made in this work.

### 5.1 Limitations.

The Multiple Inline Interacting Cantilever Beam Model was developed from the perspective of an ideal scenario. Recall the model assumed constant spacing between beams and assumed all beams were rectangular, prismatic, and had identical lengths and flexural stiffnesses. In some systems (e.g., biological or natural systems), these assumptions may not hold. As one example, the beams in many natural systems have circular cross sections and thus experience out of plane bending, torsional deflections, and sliding between adjacent beams. The solution method described above places beams directly inline and is 2-dimensional in nature. It is therefore currently unable to model out of plane deflections that occur when beams with circular cross sections contact one another. Furthermore, the solution is not applicable to systems which possess nonlinear stress–strain responses. These limitations remain the subject of future work. The purpose of this study was to take the first step and lay the foundation for that future work.

Several assumptions were made in developing the closed-form solution method. For example, it was assumed that the force bar applies a force that is perpendicular to the beam it is contacting. This assumption represents an ideal case in which no friction is present. The error arising from this assumption appears small as the closed-form solution showed good agreement with physical experiments. The assumption is likely to remain reasonable in systems with minimal friction. The model also utilizes a static perspective to analyze a dynamic process. At higher speeds, the reliability of the model will likely decrease as inertial forces become more influential [24].

Lastly, we assumed that the force vector acting between beam(*i*) and beam(*i* + 1) acts perpendicular to the end angle of beam(*i*). However, the precise direction of the force vector acting between the beams is not expressly known as there is some error in the beam end angle as calculated by the pseudo-rigid body model. A small analysis was therefore conducted to determine which of three potential assumptions would minimize the error in the closed-form solution calculation of *EI*. In particular, the 100 simulations outlined in Table 2 were evaluated using three different assumptions for the direction of the force vector between beam(*i*) and beam(*i* + 1). The percent error in *EI* across all 100 simulations was calculated. These results are summarized in Table 3.

## 6 Conclusion

A novel solution method to represent the force–deflection response of multiple, inline, interacting cantilever beams has been developed. The solution method was predicated upon the *pseudo-rigid body model* concept. Just four input parameters (*l, h, s,* and either *F _{peak}* or

*EI*) are required to solve for the response of the system. The closed-form solution was shown to provide accurate predictions when compared to physical and computational experiments. The solution method is most accurate with deflections less than 50 deg.

Due to its relative simplicity, the model can be applied to build intuition regarding numerous natural and synthetic systems (agricultural crops, forests, brushes, etc.). Future work should be conducted to improve the accuracy and effective range of the solution method. In particular, multiple pins and torsional springs could be added for increased accuracy [25], non-rectangular beam geometries should be examined, and a solution for systems which possess nonuniform beams should be explored.

## Acknowledgment

This work was supported by a Seed Grant from the University of Idaho Office of Research and Economic Development, from the National Science Foundation Award #1826715 and the United States Department of Agriculture-National Institute of Food and Agriculture award #2016-67012-28381. Any opinions, findings, conclusions, or recommendations expressed in this publication are those of the author(s) and do not necessarily reflect the view of the U.S. Department of Agriculture or of the National Science Foundation.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*h*=height of force bar

*l*=length of cantilever beams

*m*=number of maximum force peaks

*s*=beam-to-beam spacing

*t*=total number of beams in system

*u*=maximum number of beams in contact at frontmost beam’s maximum deflection

*E*=Young’s modulus (modulus of elasticity)

*F*=applied force acting on cantilever beam

*I*=cross-sectional moment of inertia

*K*=torsional spring constant

*P*=horizontal component of

*F**T*=torque

*c*=_{o}parametric angle coefficient

*d*=_{i}horizontal distance that beam(

*i*) extends past beam(*i*+ 1)*F*=_{peak}Maximum reaction force exerted on the force bar before the beam slides beneath the force bar

*F*=_{x peak}*x*-component of*F*_{peak}*K*_{Θ}=nondimensionalized torsional spring constant

*EI*=flexural stiffness

- force =
bar rigid body which deflects cantilever beams, see Fig. 1

- FEM =
finite element model

*nP*=vertical component of

*F*- PRBM =
pseudo-rigid body model

*SD*=standard deviation

- (
*α*^{2})=_{t}nondimensionalized transverse load index

*β*=90 deg—Θ

*γ*=characteristic radius factor

- Θ =
PRBM angle of deflection

*θ*_{o}=beam end angle

*ϕ*=angle of

*F*with respect to undeflected axis

### Appendix

#### Single Cantilever Beam Model

Note that the following equations are derived in detail in Ref. [16]. That derivation is briefly reported here with slight adaptations for the reader’s convenience.

*l*with a force applied at its free end as depicted in Fig. 13(a). Friction forces and the weight of the cantilever beam are considered negligible for simplification. The beam is assumed to undergo large deflections, and the stress–strain response of the material is linear elastic. The total force applied at the free end

*F*has both a horizontal component

*P*and a vertical component

*nP*. For the vertical force component, a positive

*n*describes a compression force acting on the undeflected beam. Thus,

*γl*where

*γ*is the characteristic radius factor. As shown in Ref. [16],

*γ*varies as a function of

*ϕ*, the angle of the force applied, and is defined with respect to the undeflected axis (vertical axis in this case). From the PRBM shown in Fig. 13(b))

*θ*) and the PRBM angle, Θ, expressed as

_{o}*c*, the parametric angle coefficient, is a function of

_{o}*n*. This allows

*ϕ*to be corrected from

*γ*as a function of

*n*:

*K*

_{Θ}, a nondimensionalized torsional spring constant. This stiffness coefficient is related to the transverse force, which causes the link to deflect and produce a torque at the characteristic pivot. The nondimensionalized transverse load index can be written as

*F*is the traverse force which produces a torque at the characteristic pivot. The relationship between (

_{t}*α*

^{2})

*and*

_{t}*Θ*is nearly linear, allowing the force–deflection relationship to be described as

*K*

_{Θ}can be described as a function of

*n*

*K*

_{Θ}and

*γ*are estimations with accuracy limits that are in terms of the maximum Θ, as described in Ref. [16].

*F*produces a torque at the characteristic pivot, expressed as

_{t}*K*or the flexural rigidity,

*EI*, expressed below respectively