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 [16]. 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.

Fig. 1
System of multiple inline, interacting cantilever beams: (a) undeflected state of three, vertical cantilever beams of equal length fixed directly inline along the x-axis with equal spacing (s); the force bar, oriented parallel to the z-axis, is at a known, constant height (h). (b) Side view showing the deflected state of the beams as the force bar moves in the x-direction, applying a force (F) at a constant height (h).
Fig. 1
System of multiple inline, interacting cantilever beams: (a) undeflected state of three, vertical cantilever beams of equal length fixed directly inline along the x-axis with equal spacing (s); the force bar, oriented parallel to the z-axis, is at a known, constant height (h). (b) Side view showing the deflected state of the beams as the force bar moves in the x-direction, applying a force (F) at a constant height (h).
Close modal

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.

Consider the case shown in Fig. 1, in which three identical, prismatic, cantilever beams of length 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 Fpeak. 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 β, rather than Θ, where
βi=π2Θi
(1)
where the subscript i indicates a measurement for beam(i).
Fig. 2
Side view showing the PRBM at the first beam’s most deflected state. Beams of equal length, l, are placed with equal spacing, s, along the x-axis.
Fig. 2
Side view showing the PRBM at the first beam’s most deflected state. Beams of equal length, l, are placed with equal spacing, s, along the x-axis.
Close modal
Fig. 3
Free body diagrams of the PRBM for (a) the first beam, (b), the middle beam, and (c) the last beam
Fig. 3
Free body diagrams of the PRBM for (a) the first beam, (b), the middle beam, and (c) the last beam
Close modal
To obtain an expression for Fpeak, the deflected angle of the PRBM (Θ) of each beam must be determined. Starting with the first beam
Θ1=π2sin1(hb1γ1l)
(2)
where b, the length of the fixed base of the PRBM (see Fig. 2), is expressed as
bi=(1γi)l
(3)
Solving Eq. (2) requires iteration and an assumption regarding the direction of the applied force. In particular, it is assumed the force bar imparts a force on the beam that acts perpendicular to the beam end angle, θo1. Initially, however, the force bar is assumed to apply a horizontal force to yield values of n1 = 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 n1, 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 n1 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.
Fig. 4
Diagram of iterative solution to Eq. (2)
Fig. 4
Diagram of iterative solution to Eq. (2)
Close modal
The x-coordinate of the tip of the first beam, (x1), at its maximum deflected state is
xi=γilcos(βi)
(4)
where i = 1. As shown in Fig. 2, the horizontal distance a beam extends past the next beam’s origin is expressed as
di=xis(whilexi>s)
(5)
If the xi > 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=π2tan1(γi1lsin(βi1)di1)
(6)
Recognizing that beam(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
Fi,i+1=Fi+1,i
(7)
where Fi,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 (θoi). This assumption is discussed in more detail in Sec. 5.2 Limitations. For each beam in the system, ϕi is used to determine ni which is in turn used to update each beam’s γi (see Eqs. (A4)–(A6) in the  Appendix).
With the geometry of the system determined (i.e., after solving for all the PRBM angles Θi), we can now proceed to solve for Fpeak. 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 x and y components of all force vectors are defined as
Fx=Fsin(ϕ)
(8)
Fy=Fcos(ϕ)
(9)
Solving for F32 (the force on the final beam in the system) from the free body diagram shown in Fig. 3(c) yields
F32=KΘ3cos(ϕ32)d2+sin(ϕ32)γ2lsin(β2)
(10)
or in more general terms, the final beam’s resistive force can be expressed as
Fi,i1=KΘicos(ϕi,i1)di1+sin(ϕi,i1)γi1lsin(βi1)
(11)
Now, we proceed to solve for forces acting on the middle beams in the system. First, we recognize from Eqs. (7)(9) that
F23x=F32x=F32sin(ϕ32)
(12)
F23y=F32y=F32cos(ϕ32)
(13)
Solving the free body diagram for the middle beam(s) (Fig. 3(b)) reveals that the force on all middle beam(s) in the system can be expressed as
Fi,i1=KΘi+Fi+1,i(cos(ϕi+1,i)xi+sin(ϕi+1,i)(γilsin(βi)))cos(ϕi,i1)di1+sin(ϕi,i1)γi1lsin(βi1)
(14)
Once the forces acting on the last beam and all of the middle beams have been determined, we proceed to solve for Fpeak (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 Fpeak
Fpeak=KΘ1+F21(cos(ϕ21)x1+sin(ϕ21)(γ1lsin(β1)))cos(ϕ10)x1+sin(ϕ10)γ1lsin(β1)
(15)
Alternatively, an expression for the average flexural rigidity or 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
FpeakK=Θ1+F21K(cos(ϕ21)x1+sin(ϕ21)(γ1lsin(β1)))cos(ϕ10)x1+sin(ϕ10)γ1lsin(β1)
(16)
Rearranging terms
K=Fpeak(cos(ϕ10)x1+sin(ϕ10)γ1lsin(β1))Θ1+F21K(cos(ϕ21)x1+sin(ϕ21)(γ1lsin(β1)))
(17)
It is important to note that the quantity F21/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 F21/K. This expression can then be substituted into Eq. (17) to yield an explicit equation for K. This explicit expression for K is then substituted into the following equation to determine the flexural rigidity of the beams:
EI=KlγKΘ
(18)

Note that Eq. (18) is derived in the  Appendix.

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.

Plotting the reaction force 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 (Fpeak) when the first beam has reached its maximum deflection (Fig. 6(a)). Immediately after the first beam deflects under the force bar, F will sharply drop and then continue to rise until the second beam reaches its maximum deflection (Fig. 6(b)) at which point F = Fpeak. 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 F will be less than the initial Fpeak. The reaction force F will be equal to Fpeak (m) number of times where
m=(t+1)u
(19)
t is the total number of beams and u is the maximum number of beams in contact at the frontmost beam’s maximum deflection.
Fig. 5
An idealized plot of the force response of a system six inline cantilever beams deflected by a force bar. The deflected state of the beams at points a, b, c, and d are shown in Fig. 6. Note the magnitude of force at point a is equivalent to those at point b.
Fig. 5
An idealized plot of the force response of a system six inline cantilever beams deflected by a force bar. The deflected state of the beams at points a, b, c, and d are shown in Fig. 6. Note the magnitude of force at point a is equivalent to those at point b.
Close modal
Fig. 6
Deflected states of six inline cantilever beams deflected by a force bar. (a) State of the system at the first beam’s maximum deflection. (b) State of the system at the second beam’s maximum deflection. (c) State of the system at the fifth beam’s maximum deflection. Notice the number of beams in contact has been reduced to two. (d) State of the system at the sixth beams maximum deflection.
Fig. 6
Deflected states of six inline cantilever beams deflected by a force bar. (a) State of the system at the first beam’s maximum deflection. (b) State of the system at the second beam’s maximum deflection. (c) State of the system at the fifth beam’s maximum deflection. Notice the number of beams in contact has been reduced to two. (d) State of the system at the sixth beams maximum deflection.
Close modal

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.

Table 1

Conditions for each physical experiment

SetEI (N mm2)s (mm)l (mm)h (mm)# of Force Peaks
A141,90024.91801674
B163,00024.91801654
C1111,00024.92312194
D1154,00024.92051914
E1335,00024.92312194
A241,90049.91801675
B263,00049.91801655
C2111,00049.92312195
D2154,00049.92051915
E2335,00049.92312245
A341,90019.11861753
B363,00019.11861753
C3111,00019.12372253
D3154,00019.12111983
E3335,00019.12372253
SetEI (N mm2)s (mm)l (mm)h (mm)# of Force Peaks
A141,90024.91801674
B163,00024.91801654
C1111,00024.92312194
D1154,00024.92051914
E1335,00024.92312194
A241,90049.91801675
B263,00049.91801655
C2111,00049.92312195
D2154,00049.92051915
E2335,00049.92312245
A341,90019.11861753
B363,00019.11861753
C3111,00019.12372253
D3154,00019.12111983
E3335,00019.12372253

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 (Fpeak) 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, θo. 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 (Fpeak) used to estimate 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 [1921]. 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 [1921]. The FEM formulation was validated against physical experiments conducted with one, two, three, and six beams, as shown in Fig. 7.

Fig. 7
Validation of the finite element model. The experimentally measured (solid line) and finite element model (dashed line) force–displacement curves for (a) one, (b) two, (c) three, and (d) six beam experiments.
Fig. 7
Validation of the finite element model. The experimentally measured (solid line) and finite element model (dashed line) force–displacement curves for (a) one, (b) two, (c) three, and (d) six beam experiments.
Close modal

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 mm2, 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 mm2. A summary of the input parameters for these FEM simulations is provided in Table 2.

Table 2

Model Parameters of the 100 FEM simulations (10 beam-to-beam spacings (s) × 10 maximum deflections = 100 simulations)

s (mm)h (mm)h/lMax Θ (deg)
2.0089.810.5067.2
6.2297.360.5464.0
10.44104.910.5860.4
14.67112.460.6356.8
18.89120.010.6753.2
23.11127.560.7149.4
27.33135.110.7545.4
31.56142.660.7941.1
35.78150.210.8436.5
40.00157.760.8831.3
s (mm)h (mm)h/lMax Θ (deg)
2.0089.810.5067.2
6.2297.360.5464.0
10.44104.910.5860.4
14.67112.460.6356.8
18.89120.010.6753.2
23.11127.560.7149.4
27.33135.110.7545.4
31.56142.660.7941.1
35.78150.210.8436.5
40.00157.760.8831.3

Note: Each simulation consisted of eight identical, inline beams with l = 180 mm and EI = 63,000 N mm2.

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 Fpeak 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 R2 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.

Fig. 8
Comparison of EI values: 3-pt bending versus closed-form solution method. Linear regression analysis (blue line) demonstrated the closed-form solution method is accurate (R2 = 0.99). Error bars represent the standard deviation of EI measured via 3-pt bending experiments. The red line shows the theoretical, ideal case of y = x to indicate accuracy.
Fig. 8
Comparison of EI values: 3-pt bending versus closed-form solution method. Linear regression analysis (blue line) demonstrated the closed-form solution method is accurate (R2 = 0.99). Error bars represent the standard deviation of EI measured via 3-pt bending experiments. The red line shows the theoretical, ideal case of y = x to indicate accuracy.
Close modal

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 (Fx) 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.

Fig. 9
Force–displacement comparisons. Fx plotted against x displacement for six inline cantilever beams (EI = 63,000 N mm2, l = 180 mm, h = 166 mm, and s = 24.9 mm) from a physical experiment (shown in blue), from an FEM simulation (shown in green), and from the closed-form solution (CFS) (shown in orange). The closed-form solution with a correction factor of 1.0085 multiplied to h is shown in (red).
Fig. 9
Force–displacement comparisons. Fx plotted against x displacement for six inline cantilever beams (EI = 63,000 N mm2, l = 180 mm, h = 166 mm, and s = 24.9 mm) from a physical experiment (shown in blue), from an FEM simulation (shown in green), and from the closed-form solution (CFS) (shown in orange). The closed-form solution with a correction factor of 1.0085 multiplied to h is shown in (red).
Close modal

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%).

Fig. 10
Contour plot of percent error (%) in EI between the closed-form solution and FEM as a function of beam-to-beam spacing (mm) and maximum beam deflection (deg)
Fig. 10
Contour plot of percent error (%) in EI between the closed-form solution and FEM as a function of beam-to-beam spacing (mm) and maximum beam deflection (deg)
Close modal

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.

Fig. 11
Contour plot showing the numerical derivative of EI (normalized) with respect to h/l for combinations of h/l and s/l where EI is normalized to the maximum EI found across all simulations
Fig. 11
Contour plot showing the numerical derivative of EI (normalized) with respect to h/l for combinations of h/l and s/l where EI is normalized to the maximum EI found across all simulations
Close modal
Fig. 12
Contour plot showing the numerical derivative of EI (normalized) with respect to s/l for combinations of s/l and h/l where EI is normalized to the maximum EI found across all simulations. Note the difference in colormap scales compared to Fig. 11.
Fig. 12
Contour plot showing the numerical derivative of EI (normalized) with respect to s/l for combinations of s/l and h/l where EI is normalized to the maximum EI found across all simulations. Note the difference in colormap scales compared to Fig. 11.
Close modal

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.

Table 3

Direction of force vector between adjacent cantilever beams

Force toMean error in EI (%)Standard deviation (%)
Beam(i)4.194.08
Beam(i+1)7.244.17
Beam(i)+Beam(i+1)25.613.90
Force toMean error in EI (%)Standard deviation (%)
Beam(i)4.194.08
Beam(i+1)7.244.17
Beam(i)+Beam(i+1)25.613.90

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 Fpeak 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

  •  
  • co =

    parametric angle coefficient

  •  
  • di =

    horizontal distance that beam(i) extends past beam(i + 1)

  •  
  • Fpeak =

    Maximum reaction force exerted on the force bar before the beam slides beneath the force bar

  •  
  • Fx peak =

    x-component of Fpeak

  •  
  • 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.

Consider a single, prismatic cantilever beam of length 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,
F=Pn2+1
(A1)
As shown in Ref. [16], the deflection of the beam pictured in Fig. 13(a) can be represented by the PRBM, shown in Fig. 13(b). In particular, the flexible cantilever beam is represented by two connecting rigid links with a pin joint and a torsional spring at the characteristic pivot. The characteristic pivot is positioned so that the longer link has a length equivalent to the characteristic radius of the approximately circular path of the free end of the beam. The characteristic radius is written as γ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))
ϕ=tan1(1n)
(A2)
Fig. 13
(a) Cantilever beam undergoing large deflection due to a force applied at the free end with (b) its PRBM representation
Fig. 13
(a) Cantilever beam undergoing large deflection due to a force applied at the free end with (b) its PRBM representation
Close modal
Howell [16] shows that there is a nearly linear relationship between the beam end angle (θo) and the PRBM angle, Θ, expressed as
θo=coΘ
(A3)
in which co, the parametric angle coefficient, is a function of n. This allows ϕ to be corrected from
ϕ=π2θo
(A4)
Rearranging Eq. (A2) presents
n=1tanϕ
(A5)
The following piecewise function can be used to define γ as a function of n:
γ={0.8416550.0067807n+0.000438n2(0.5<n<10.0)0.8521440.0182867n(1.8316<n<0.5)0.912364+0.0145928n(5.0<n<1.8316)
(A6)
The PRBM represents the beam’s resistance to deflection with the stiffness coefficient, 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
(α2)t=Ftl2EI
(A7)
where Ft is the traverse force which produces a torque at the characteristic pivot. The relationship between (α2)t and Θ is nearly linear, allowing the force–deflection relationship to be described as
(α2)t=KΘΘ
(A8)
KΘ can be described as a function of n
KΘ={3.024112+0.121290n+0.003169n21.9676472.616021n3.738166n2(5<n2.5)2.649437n30.891906n40.113063n52.6548550.509896×101n+0.126749×101n2(2.5<n1)0.142039×102n3+0.584525×104n4(1<n10)
(A9)
Note that both KΘ and γ are estimations with accuracy limits that are in terms of the maximum Θ, as described in Ref. [16].
The transverse force Ft produces a torque at the characteristic pivot, expressed as
T=Ftγl
(A10)
This torque can also be written as
T=KΘ
(A11)
where K is the torsional spring constant. Setting Eqs. (A10) and (A11) equal and rearranging for Ft yields
Ft=KΘγl
(A12)
Combining Eqs. (A7) and (A8)
Ftl2EI=KΘΘ
(A13)
Inserting Eq. (A12) into the above then gives
KΘlγEI=KΘΘ
(A14)
Equation (33) can then be rearranged for K or the flexural rigidity, EI, expressed below respectively
K=γKΘEIl
(A15)
EI=KlγKΘ
(A16)

References

1.
Barten
,
H. J.
,
1944
, “
On the Deflection of a Cantilever Beam
,”
Quarterly Appl. Math.
,
2
(
2
), pp.
168
171
. 10.1090/qam/10879
2.
Barten
,
H. J.
,
1945
, “
Corrections to My Paper on the Deflection of a Cantilever Beam
,”
Quarterly Appl. Math.
,
3
(
3
), pp.
275
276
. 10.1090/qam/13361
3.
Bisshopp
,
K. E.
, and
Drucker
,
D. C.
,
1945
, “
Large Deflection of Cantilever Beams
,”
Quarterly Appl. Math.
,
3
(
3
), pp.
272
275
. 10.1090/qam/13360
4.
Frisch-Fay
,
R.
,
1961
, “
A New Approach to the Analysis of the Deflection of Thin Cantilevers
,”
ASME J. Appl. Mech.
,
28
(
1
), pp.
87
90
. 10.1115/1.3640472
5.
Navaee
,
S.
, and
Elling
,
R. E.
,
1992
, “
Equilibrium Configurations of Cantilever Beams Subjected to Inclined End Loads
,”
ASME J. Appl. Mech.
,
59
(
3
), pp.
572
579
. 10.1115/1.2893762
6.
Wang
,
T. M.
,
1968
, “
Nonlinear Bending of Beams With Concentrated Loads
,”
J. Franklin Inst.
,
285
(
5
), pp.
386
390
. 10.1016/0016-0032(68)90486-9
7.
Laurance
,
W. F.
, and
Curran
,
T. J.
,
2008
, “
Impacts of Wind Disturbance on Fragmented Tropical Forests: A Review and Synthesis
,”
Austral Ecol.
,
33
(
4
), pp.
399
408
. 10.1111/j.1442-9993.2008.01895.x
8.
Cook
,
D. D.
,
Meehan
,
K.
,
Asatiani
,
L.
, and
Robertson
,
D. J.
,
2020
, “
The Effect of Probe Geometry on Rind Puncture Resistance Testing of Maize Stalks
,”
Plant Meth.
,
16
(
1
), pp.
1
11
. 10.1186/s13007-019-0534-5
9.
Stubbs
,
C.
,
Seegmiller
,
K.
,
McMahan
,
C. S.
,
Sekhon
,
R.
, and
Robertson
,
D. J.
,
2020
, “
Diverse Maize Hybrids are Structurally Inefficient at Resisting Wind Induced Bending Forces That Cause Stalk Lodging
,”
Plant Meth.
,
16
(
1
), pp.
1
15
. 10.1186/s13007-020-00608-2
10.
Berry
,
P. M.
,
Sterling
,
M.
, and
Mooney
,
S. J.
,
2006
, “
Development of a Model of Lodging for Barley
,”
J. Agronomy Crop Sci.
,
192
(
2
), pp.
151
158
. 10.1111/j.1439-037X.2006.00194.x
11.
Stubbs
,
C. J.
,
Baban
,
N. S.
,
Robertson
,
D. J.
,
Alzube
,
L.
, and
Cook
,
D. D.
,
2018
, “Bending Stress in Plant Stems: Models and Assumptions,”
Plant Biomechanics
,
A.
Geitmann
and
J.
Gril
, eds.,
Springer
,
Cham
, pp.
49
77
.
12.
Gardiner
,
B.
,
Byrne
,
K.
,
Hale
,
S.
,
Kamimura
,
K.
,
Mitchell
,
S. J.
,
Peltola
,
H.
, and
Ruel
,
J. C.
,
2008
, “
A Review of Mechanistic Modelling of Wind Damage Risk to Forests
,”
Forestry
,
81
(
3
), pp.
447
463
. 10.1093/forestry/cpn022
13.
Schelhass
,
M. J.
,
Kramer
,
K.
,
Peltola
,
H.
,
Van Der Werf
,
D. C.
, and
Wijdeven
,
S. M. J.
,
2007
, “
Introducing Tree Interactions in Wind Damage Simulation
,”
Ecol. Modell.
,
207
(
2–4
), pp.
197
209
. 10.1016/j.ecolmodel.2007.04.025
14.
Erndwein
,
L.
,
Cook
,
D. D.
,
Robertson
,
D. J.
, and
Sparks
,
E. E.
,
2020
, “
Field-Based Mechanical Phenotyping of Cereal Crops to Assess Lodging Resistance
,”
Appl. Plant Sci.
,
8
(
8
), p.
e11382
. 10.1002/aps3.11382
15.
Cook
,
D. D.
,
de la Chapelle
,
W.
,
Lin
,
T. C.
,
Lee
,
S. Y.
,
Sun
,
W.
, and
Robertson
,
D. J.
,
2019
, “
DARLING: A Device for Assessing Resistance to Lodging in Grain Crops
,”
Plant Meth.
,
15
(
1
), p.
102
. 10.1186/s13007-019-0488-7
16.
Howell
,
L. L.
,
2001
,
Compliant Mechanisms
,
John Wiley & Sons
,
New York
, pp.
135
160
.
17.
Howell
,
L. L.
,
Midha
,
A.
, and
Norton
,
T. W.
,
1996
, “
Evaluation of Equivalent Spring Stiffness for Use in a Pseudo-rigid-Body Model of Large-Deflection Compliant Mechanisms
,”
ASME J. Mech. Des.
,
118
(
1
), pp.
126
131
. 10.1115/1.2826843
18.
Thornton
,
S. T.
, and
Marion
,
J. B.
,
2004
,
Classical Dynamics of Particles and Systems
, 5th ed,
Brooks/Cole
,
Belmont, CA
.
19.
Matthews
,
F. L.
,
Davies
,
G. A. O.
,
Hitchings
,
D.
, and
Soutis
,
C.
,
2000
,
Finite Element Modelling of Composite Materials and Structures
,
Elsevier
,
New York
.
20.
Hibbitt
,
K.
,
Karlsson
,
B. I.
, and
Sorenson
,
E. P.
,
2016
,
ABAQUS/Standard Theory Manual
,
Dassault Systemes Simulia Corp.
,
Johnston, RI
.
21.
Smith
,
M.
,
2016
,
ABAQUS Analysis Manual
,
Dassault Systemes Simulia Corp.
,
Johnston, RI
.
22.
Lawlor
,
D. A.
,
Tilling
,
K.
, and
Smith
,
G. D.
,
2016
, “
Triangulation in Aetiological Epidemiology
,”
Int. J. Epidemiol.
,
45
(
6
), pp.
1866
1886
. 10.1093/ije/dyw314
23.
Nelson
,
N.
,
Stubbs
,
C. J.
,
Larson
,
R.
, and
Cook
,
D. D.
,
2019
, “
Measurement Accuracy and Uncertainty in Plant Biomechanics
,”
J. Exp. Bot.
,
70
(
14
), pp.
3649
3658
. 10.1093/jxb/erz279
24.
Li
,
N.
,
Su
,
H.-J.
, and
Zhang
,
X.-P.
,
2017
, “
Accuracy Assessment of Pseudo-rigid-Body Model for Dynamic Analysis of Compliant Mechanisms
,”
ASME J. Mech. Rob.
,
9
(
5
), p.
054503
. 10.1115/1.4037186
25.
Su
,
H.-J.
,
2009
, “
A Pseudorigid-Body 3R Model for Determining Large Deflection of Cantilever Beams Subject to Tip Loads
,”
ASME J. Mech. Rob.
,
1
(
2
), p.
021008
. 10.1115/1.3046148