Simple approximate formulas for postbuckling deflection of heavy elastic columns

Columnar buckling is a ubiquitous phenomenon that occurs in both living things and man-made objects, regardless of the length scale ranging from macroscopic to nanometric structures. In general, analyzing the post-buckling behavior of a column requires the application of complex mathematical methods because it involves nonlinear problem solving. To complement these complex methods, this study presents simple analytical formulas for the large deflection of a heavy elastic column under combined loads. The analytical formulas relate the concentrated load acting on the tip of the column, the column's own weight, and the deflection angle of the column through a simple mathematical expression. This can assist in obtaining an overall picture of the post-buckling behavior of heavy columns from an application point of view.


I. INTRODUCTION
When a vertical slender column is axially compressed, it may suddenly flex archwise and lose its ability to withstand the load. This deflection, called columnar buckling, is a crucial consideration in structural design because the column becomes unstable under compressive stress considerably lesser than the failure strength of the constituent materials 1 . In addition to the field of structural design, columnar buckling is relevant in a wide range of fields in which slender column structures are involved; plant science, nanomaterial engineering 2 , and robotics 3 are only a few to mention. In plant science, for example, columnar buckling has been presumed to impose an upper limit on the growth height of upright stems and trunks 4-6 . Plants are subjected to self-weight-induced downward load and concentrated loads caused by e.g., the petals, leaves, or rain water accumulated near the tip; therefore, plant growth is considered to be regulated for avoiding buckling under these combined loads. Further, it has been suggested that columnar buckling occurs in nanometric tubes and wires [7][8][9][10][11][12] , whose mechanical stability is crucial for improving the performance of nanoelectromechanical devices [13][14][15] .
From a mathematical point of view, columnar buckling is the bifurcation of the solution of a nonlinear differential equation. This nonlinear equation governs the static equilibrium states of an elastic column under axial compression. Bifurcation signifies the transition of the equilibrium column shape from a vertically standing state to a deflected (postbuckled) one. This transition occurs at critical values of the variable coefficients involved in the nonlinear equation, where the coefficients characterize the compression strength. Due to the nonlinearity, this differential equation cannot be solved in an elementary way, unless it is approximated that the column deflection is exceedingly small. Therefore, various methods have been proposed for solving it to realize high-accuracy analysis of columnar buckling behaviors under different mechanical conditions, including the geometry, loading conditions, and boundary conditions 2,  . Most meth-ods have been based on either the Bessel function, elliptic integral, or other nonelementary functions/integrals, ensuring high accuracy in the analysis of column buckling behavior even with large deflections; those based on generalized hypergeometric functions 29 and the singular perturbation method 31 are only a few to mention. While such the sophisticated techniques are of fundamental importance, the mathematical formulas derived from them are complex and often beyond intuitive understanding. Using a simpler more manageable formula at the expense of a certain numerical accuracy helps in obtaining an overall picture of the buckling behavior, particularly from an application perspective.
In this study, simple analytical relationships between the concentrated load acting on the tip of the column, the column's own weight, and the column deflection angle are formulated. The estimated formulas concisely reflect the contribution of the tip-concentrated load and selfweight to the tilt angle at the column tip, respectively. In addition, it is numerically observed that the tip-tilt angle is almost a square root of the load increment within the postbuckling state. All the results are within the realm of elementary calculus, providing an intuitive perspective of columnar buckling under combined loading.

A. Governing equation
We consider a vertically built heavy elastic column with a free end at the tip. Figure 1(a) shows the configuration of a column subjected to downward concentrated load at the tip as well as uniformly distributed weight. Under the combined load condition, the local bending moment, M (ℓ), satisfies the relation of Here, F is the concentrated downward load applied at the tip, ρ is the self-weight of the column per axial length, and L is the column length; ℓ is the arc length measured from the bottom end, and θ is the local tilt angle. Angle θ is defined such that θ = 0 in the vertical direction and θ = π/2 in the horizontal; see Fig. 1 The column is sufficiently slender such that the local curvature is proportional to the local bending moment M (ℓ) as where E is the Young's modulus and I is the second moment of area of the column. For convenience, we introduce three dimensionless parameters 18 : where α and β represent the relative importance of the concentrated load F and self-weight ρL, respectively, to the flexural rigidity EI. The normalized arc length s is defined such that s = 0 at the tip and s = 1 at the bottom. Using these three dimensionless parameters, the governing equation is derived from Eqs. (1) and (2) as and the associated boundary conditions are θ(s) = 0 at s = 1 (i.e., at the bottom), (5) dθ(s) ds = 0 at s = 0 (i.e., at the tip).

B. Approximation solution based on series expansion
Our immediate task is to derive an approximate solution for Eq. (4) with respect to θ(s), depending on α and β. To realize this, we expand θ(s) using a Maclaurin series up to the n-th order as follows: where θ (k) 0 indicates that d k θ(s)/ds k at s = 0. It is shown below that all the derivatives θ (k) 0 on the right of Eq. (7) can be expressed using functions of α, β, and θ 0 . First, the boundary condition at s = 0 is Next, from Eq. (4), For the higher-order derivatives, we use Eq. (4) to reduce the order of differentiation; for instance, implying that Repeating this procedure up to the ninth order, we obtain In the following discussion, the series expansion is truncated at the ninth degree for securing the numerical accuracy in engineering application. This is because, when α is sufficiently smaller than the unit, the terms of θ (7) 0 and θ (8) 0 hardly contribute to θ(s) as they have a common factor of α. Similarly, when β is nearly equal to zero, the terms of θ 0 , as well as θ give only minor contributions to θ(s) as they have a common factor of β. Therefore, the series expansion up to the ninth degree is necessary to include the effects of higher-order terms even if either α or β is much smaller than the unit.
It is to be noted that the boundary condition at s = 1 requires θ as proved by imposing Eq. (5) into Eq. (7). All θ hence, Eq. (18) relates the tip-tilt angle θ 0 to the loading conditions determined by an α and β pair. Equation (18) may have several roots of θ 0 depending on the values of α and β. Among them, the least non-negative root is the physically relevant tip-tilt angle of the column. Through actual calculation, we solved Eq. (18) using the bisection method to find the least root. After obtaining θ 0 , θ(s) can be calculated using Eq. (7) for any value of s; thereby, the deflection curve of the column can be drawn for a given α and β. If the origin of the x-y coordinate plane is set at the bottom end of the column, the deflection curve in the x-y plane can be parametrically expressed as The integration was numerically performed according to Simpson's 3/8 rule.

III. RESULT
A. Relationship between the combined loads and tip-tilt angle Figure 2 depicts the typical deflection behaviors of the column for various α and β. When α is fixed at α = 1.0, as presented in Fig. 2(a), the column remains vertical until β reaches a critical value of β c = 4.75. When β exceeds β c , the column begins to deflect and the tiptilt angle θ 0 increases monotonically with β. Note that the threshold value β c changes, if α is fixed at a different value, as demonstrated later. A similar deflection process is observed when β is fixed; see Fig. 2(b). For instance, when β is fixed at β = 1.0, the critical value of α is α c = 2.16 beyond which θ 0 increases with α. Again, the value of α c depends on the choice of β.
The contour map in Fig. 3 illustrates the variation of the tip-tilt angle θ 0 with the changes in the values of α and β. The upper limit of the contour line is set to 1.6 (unit: radian); data above this limit are omitted in the plot. It can be clearly observed that θ 0 = 0 within the left-bottom triangular region (gray colored). This indicates that for relatively small α and β, the vertically standing state is energetically favored; hence, the column is not deflected but merely compressed in the axial direction. In the rainbow-colored region, in contrast, θ 0 > 0, implying that deflection will occur. The boundary curve that separates the two equilibrium states θ 0 = 0 and θ 0 > 0 is found to be a nearly linear slanted line in the α-β plane, which has intersections with the α-axis and β-axis, respectively, at β c = 7.815 and α c = 2.468. The latter value is close to the normalized Euler buckling load π 2 /4 with deviation of 0.024%. If the series expansion in Eq. (7) were truncated at the seventh or fifth degree, for example, the deviation from teh Euler value would be enlarged to 0.097% or 2.78%, respectively. It is also noted that our values of α c and β c obtained by the ninthdegree expansion are nearly equal to those evaluated in the existing work with differences less than 1% 29,31 . An important observation in Fig. 3 is that all the contour lines are nearly parallel to each other, approximately obeying the following relationship: The series of power-law curves are the edges of the cutting plane that appear when the three-dimensional θ 0 (α, β) surface spanning within the α-β-θ 0 space (highlighted by the rainbow-colored portion in Fig. 3) is cut along a plane parallel to the α-axis (and perpendicular to the page). All the curves are well fitted to the following expression: particularly in the θ 0 ≤ 1.0 range. The threshold α c decreases monotonically with the increase in β, as observed in the contour map in Fig. 3. In contrast, the values of the coefficient A and exponent p α are found to be insensitive to the change in β; they are approximately estimated as A ≃ 1.6 and p α ≃ 0.48, with slight dependence on β, as shown in Fig. 4(c). Another class of θ 0 power-law behaviors be deduced when α is fixed and β is tuned. Figure 4(b) depicts the power-law behavior of θ 0 as a function of β for fixed α.
The fitting curves obey the following expression: with B ≃ 0.9 and p β ≃ 0.48. Similar to the previous case, both B and p β are insensitive to the variation in the fixed value of α, as shown in Fig. 4(d).
Emphasis should be placed on the fact that the optimal values of p α and p β , both nearly equal to 0.5, are robust against the variations of α and β that characterize the loading conditions. This implies that the tip-tilt angle θ 0 is almost a square root of the increment, if either the concentrated or distributed load is amplified gradually. These nearly square-root behaviors of θ 0 with the load increment can be approximately expressed as Using these approximate formulas, the observed increase in θ 0 with the increase in load can be estimated easily by measuring the initial value of θ 0 before increasing the load. These square root formulas are the second main finding of this study.

IV. DISCUSSION
We determined that the phase boundary curve that separates the vertically standing phase (i.e., θ 0 = 0) and deflection phase (θ 0 > 0) in the α-β plane can be approximately expressed by a simple linear equation (Eq. (21)).
To determine the critical values of α and β with higher accuracy, the following formula can be used 43 : Using this accurate formula of Eq. (25), we can obtain the α-intercept of the boundary curve at α c = π 2 /4, which is consistent with our approximate value of α c = 2.46. Similarly, the β-intercept derived from the accurate formula is β c = 7.83735, which is again proximate to our result of β c = 7.82. This comparison of the critical values confirms the practicality of our simple expression for the phase boundary, given by Eq. (21). One of the interesting ways of using the approximate formulas derived in this study is to diagnose the mechanical stability of upright plants, such as large trees and tall bamboos. Understanding the stability of the tree trunk under self-weight and the applied loads caused by the leaf and branches is crucial when developing management strategies to reduce the risk of damage by abiotic agents. In addition, bamboo is a unique plant that does not break, exhibiting considerable flexure, even when subjected to load 44 ; it is extensively used as building material in areas such as Southeast Asia. Note that when using mathematical expressions in actual stability assessment, rapid estimation of the stability using a simple formula is desired, rather than exact numerical data based on expensive numerical simulation. Our simple formula for heavy column deflection facilitates an approximate perspective of the stability analysis of plants in a postbuckling state. However, to investigate the mechanical response of a plant more quantitatively, it is necessary to consider the variation in both the mechanical and geometrical properties of the plant from the ground to the top as well as the effect of initial imperfection. In this context, we need to extend our formulation considering the spatial variation in the cross-sectional geometry 33,45 , that in Young's modulus 38 , and the presence of cracks in the initial state 35 ; I intended to address these interesting problems in future.
The effect of boundary condition variation on the postbuckling behavior is another interesting topic. In this study, our focus was restricted to the deflection of a column in which one end is fixed and the other is free. Another important case to consider is a column in which both ends are fixed because it is frequently used in structural applications, in both macroscopic and nanoscopic [13][14][15] systems. Imposing different boundary conditions on a column can inherently change its mechanical response to the applied loads. Therefore, as in the case of this study, it is interesting to examine whether an approximate expression for the deflection can be derived from a column with fixed conditions at both ends. Our results can pave the way for exploring simple approximate formulas describing column deflection under various end conditions.

V. CONCLUSION
Simple approximate formulas that describe the postbuckling behavior of columns under a tip-concentrated load F and uniformly distributed load ρ∆ℓ were derived in this study. The main findings include the 0.3-rule for the balance of the two load contributions to the deflection, given by Eq. (21), and the nearly square root evolution of the deflection angle, given by Eq. (24). These simple formulas relate the degree of deflection to each combined load through a basic function, enabling fast estimation and intuitive understanding of the increase in the deflection angle.

ACKNOWLEDGMENTS
I would like to thank A. Takashima for the technical support in manipulating numerical data. Extensive discussions with K. Ishikawa, S. Tsugawa, A. In-oue, Y. Umeno, and M. Sato have been illuminating for coming up with the idea of this present study. This work was supported by JSPS KAKENHI Grant Numbers 18H03818, 19H02020, 19H05359, and 19K03766.