Mathematical modelling to determine the greatest height of trees.

This study aimed to analyse the critical height of a column whose weight varies vertically in order to obtain a simple scaling law for a tree where the weight distribution considered. We modelled trees as cantilevers that were fixed to the ground and formulated a self-buckling problem for various weight distributions. A formula for calculating the critical height was derived in a simple form that did not include special functions. We obtained a theoretical clarification of the effect of the weight distribution of heavy columns on the buckling behaviour. A widely applicable scaling law for trees was obtained. We found that an actual tree manages to distribute the weight of its trunk and branches along its vertical extent in a manner that adequately secures its critical height. The method and findings of this study are applicable to a wide range of fields, such as the simplification of complicated buckling problems and the study of tree shape quantification.

Introduction

Over the course of evolution, living organisms have acquired a remarkable degree of wisdom. For example, standing plants in nature have acquired mechanisms to efficiently resist various external forces, such as bending, shear, and buckling[1,2]. Many previous studies have focused on the bamboo, which has a peculiar structure, investigating the rationality of this structure[3,4]. A wild bamboo has an extremely long and slender body. The shape as such possesses a poor resistance to buckling. However, previous studies found that the bamboo reduces self-weight by incorporating a hollow structure, such that its buckling resistance is effectively improved by adjusting the node interval and the vascular bundle distribution.

Based on those prior results, this study focused on trees, which are self-standing plants similar to bamboo. The bamboo consists of a hollow cylinder with light branches and leaves. In contrast, trees have a solid cross-section and a heavy body. Therefore, the buckling risk under self-weight highly increases as trees grow taller. Nevertheless, trees in the wild naturally acquire appropriate heights and concomitant levels of mechanical stability that are adapted to harsh natural environments. This implies that mechanical strategies for avoiding self-buckling are incorporated in the forms adopted by trees.

Greenhill[5] conducted a study on the self-buckling of trees. He modelled a tree as a cylindrical cantilever and derived the exact solution for the critical height for self-buckling as follows:where is the modulus of elasticity , is the unit volume weight , is the radius , and is a constant (). Equation (1) shows that the critical height of a tree is proportional to the 2/3 power of the radius. McMahon proved that this law is correct for wild trees[6]. This scaling law has been applied widely in forest science and ecology[7-10] owing to its simplicity. von Karman and Biot[11] solved the governing differential equation by using a series solution and derived a formula for the critical height, which is almost equivalent to Greenhill's equation ().

Other factors that determine critical height are hydraulic conditions[12,13], wind force effects[14,15], and genetic factors[16]. The scaling laws found by Greenhill correspond to the findings of McMahon[6], which indicated that the equation formulating the critical height for self-buckling includes scaling laws that are widely applicable to trees. Most of the prior studies of the scaling laws associated with trees considered that the critical height is determined by the self-buckling condition[17-19].

Greenhill's formula is widely applied in the fields of forestry science and ecology, and it is also the starting point for research on self-weight buckling in civil and mechanical engineering. Grishcoff[20] investigated a column that is loaded simultaneously with its self-weight and concentrated loads. Wang and Drachman[21] extended Grishcoff's study to an inverted column. Whereas these solutions were approximate, Chai and Wang[22] and Duan[23] used different methods to obtain the exact buckling stresses under self-weight and concentrated loads. Zorica[24] extended this problem to a diagonal cantilever beam.

The self-buckling problem of columns has been extended to tapered columns, which have recently been actively used in civil engineering structures[25]. A study of tapered columns was initiated by Smith[26]. He formulated the buckling problem for linearly tapered columns, which ignores the self-weight. Subsequent studies extended the problem scope to heavy columns with arbitrary taper ratios, using various approaches[27-29]. An investigation of curvilinearly-tapered columns[30], another investigation that introduces a model that includes concentrated loads at arbitrary positions[31], theoretical consideration regarding post-buckling behaviours of heavy columns[32-34], and studies employing other considerations[35,36] have been performed. However, these studies are highly complicated, and none of them expresses a simple law such as Greenhill's formula.

Previously, we theoretically derived a simple scaling law for a heavy tapered column for the first time[37] with the aim of applications in engineering as well as in ecology and forest science. We found that the critical height of a linearly tapered column is inversely proportional to the 1/6 power of the ratio of the radius of the free end to the radius of the fixed end. We also found that the Greenhill power law holds even in a tapered column.

As aforementioned, the critical buckling height of trees and the self-buckling of columns have been investigated for various purposes and for various fields of application. However, the effects of branches and leaves, which are considered important factors for buckling resistance, were not considered in previous studies. A clarification of the relationship between the critical height and the weight and distribution of branches and leaves will have significant engineering implications, such as how to distribute weight in the design of tower structures and how to control or reduce the occurrence of buckling under self-weight. Furthermore, a simple scaling law for the buckling length and the leaf weight of trees, similar to Greenhill’s equation, would be of great value in forestry science, ecology, and other fields of research.

The purpose of this study is to derive the critical height of a column whose weight varies in the vertical direction. We modelled the trees as cantilevers that were fixed to the ground. We formulated a self-buckling problem that took account of various weight distributions by using the series solutions of von Karman and Biot[11]. In considering real trees, the weight distribution of the trunk of a tree was assumed to be constant, and the branches were assumed to be variable quantities in the vertical direction. The formula for the critical height was derived in a simple form that did not include special functions, similar to the Greenhill formula. The effect of the weight distribution of heavy columns on the buckling behaviour was theoretically clarified, and a widely applicable scaling law was obtained.

The methods and results presented herein provide the buckling length of a heavy column considering its weight distribution. They also simplify the treatment of conventional complex buckling problems in engineering[29,30]. Furthermore, they can be used to investigate the rationality of tree morphology, to study tree shape quantification in forest science[38,39], as well as in a wide range of other applications.

Methods

Density function

In previous studies, the density of a tree was treated as constant in the vertical direction (). However, in this study, the density was treated as a function of height to clarify the relationship between the weight distribution and the critical height of a tree. We considered the following density functions:where is the density at the reference point and is a parameter that controls the shape of the density distribution. An example of the distribution of each model is shown in Fig. 1. The coordinate axis was aligned with the neutral axis, with at the free end and at the fixed end.

Figure 1

Density functions for models A, B, C, and D, respectively.

Model A is a model in which the density at the free end is constant and independent of parameter () whereas the density at the fixed end is times the density at the free end () . Model B in Eq. (3) is the inverted form of Model A, i.e., the density at the fixed end is constant ) whereas the density at the free end is times the density at the fixed end (). For the same value of in Models A and B, both models have the same total weight. Based on the observation that the total weight of the branches rarely exceeds the total weight of the trunk[40,41], we assumed that the range of in Models A and B was . Models C and D in Eqs. (4) and (5) are models in which the total weight of the system is independent of the parameter if the height is invariant. The difference between Models C and D is whether the distribution is linear or curved. In these models, the density distribution is constant in the vertical direction, as in Greenhill’s equation with . Moreover, the weight is concentrated at the top when and at the bottom when . A summary of these results appears in Table 1. Not only to clarify the effect of weight distribution on the critical height, but also to ensure future expandability, we used the above four density distribution functions as they are geometrically simple and can represent various distribution shapes.

Table 1

Calculation models and related numerical conditions.

Model	Function	Range of n	Constant-density state
A	ρx=n-1l_cx+1ρ_0	0≤n≤3	n=1
B	ρx=1-nl_cx+nρ_0	0≤n≤3	n=1
C	ρx=2nl_cx+1-nρ_0	-1≤n≤1	n=0
D	ρx=1-ncosπl_cxρ_0	-1≤n≤1	n=0

The calculation model was a cantilevered elastic column with a radius and a length (see Fig. 2). Greenhill's formula[5] that was verified by McMahon[6] and has been broadly used in various research, it had been derived by modelling trees as non-tapered cylinders. Furthermore, the geometrically simple models are desirable to obtain the wide applicability and extensibility. For the above reasons, we modelled trees as non-taper cylinders. The weight from the upper end to any point is given by:where is the cross-sectional area . When the column buckles under its self-weight, the shear force at any point is given by:where is the deflection angle. Note that the deflection angle is not included in the integration of Eq. (6) with respect to because this is only a volume integration to determine the weight . If we assume that is very small, can be approximated as . Therefore, Eq. (7) can be written as

Figure 2

Calculation model.

From the elastic curve equation of the beam, the bending moment at any point can be obtained aswhere is the deflection. Because the deflection angle is very small, Eq. (9) can be approximated as:

Governing differential equation

From the relationship between the shearing force of Eq. (7) and the end bending moment of Eq. (10), the governing differential equation is obtained as

Equation (11) is a second-order ordinary differential equation that has an independent variable and dependent variable . We integrated the equations that are obtained by substituting the density functions in Eqs. (2)–(5) for in Eq. (11) and then applied the following transformation:where is a constant.

The governing differential equations were obtained in the following form:

If we consider the governing equations in their simplest form, the constant for all the models can be defined as

Critical height formula

The general solution of the governing equations was obtained by the series solution method using Mathematica. For example, in Model C, the series solution of Eq. (15) is obtained as follows:where and are arbitrary constants.

We applied the following boundary conditions of the cantilever to the general solution of the governing differential equations, deriving the following algebraic equation for the critical height:

By applying the boundary condition at the free end, we found that in Eq. (18). Moreover, by applying the boundary condition in the fixed end, Eq. (18) can be rewritten as

By considering the condition for finding a non-trivial solution (), we obtained the following algebraic equation with respect to :

When , Eq. (21) corresponds to the equation of von Karman and Biot[11]. The critical height can be determined for any value of by calculating the smallest positive real number that satisfies Eq. (21) and substituting that value of into the following equation (which follows from Eqs. (11), (15), and (17)):where is Greenhill’s constant (). Equation (22) is expressed in the form of the product of the coefficient varying with the density distribution and the critical height in the constant-density column introduce by Greenhill. Equation (21) has only two variables, and . In other words, varies only with . This property is the same for all density models, and the critical height can be obtained using the same method.

Buckling mode

In this section, we describe the calculation of the buckling modes. First, we integrated Eq. (18) with respect to , in order to obtain the deflection :where is an arbitrary constant. By considering the boundary condition at the fixed end (), is given by:

Next, we substituted Eq. (24) into Eq. (23) and applied the condition , The deflection was obtained as follows:

Because Eq. (25) includes the unknown coefficient we used the deflection ratio to eliminate it. The maximum deflection appears at , and the deflection ratio to draw the buckling modes is given by

Therefore, by substituting into Eq. (26), the buckling modes under self-weight could be obtained. We obtained the buckling modes for other density models with the same method.

Example: application to real trees

In this section, we consider the effect of weight distribution on critical height in real trees by using the critical height formula derived by the aforementioned methods. Based on the investigation on the weight distribution of real trees[40,41], we used the following functions to add the trunk weight to Models C and D, respectively:where is the density of the branches [kg/m3], and is the density of the trunk [kg/m3]. An example of the distributions in Eqs. (27) and (28) is shown in Fig. 3. In this study, we expressed the density ratio of the branches and the trunk, , as a weight ratio because the density ratio of the branches and trunk is equal to their weight ratio. By using the same method as for the other models, the governing differential equations were obtained as follows:

Figure 3

Examples of density distributions representing real trees.

Calculation

Method for solving the critical height equation

Algebraic solving of high-order equations such as Eq. (21) is impossible. Therefore, we used numerical computation to find the smallest positive real-valued solution that satisfies the high-order equation presented in this study. For the numerical method, we adopted the secant method, which does not require differentiation to allow for future extensions.

We rewrote Eq. (18) as follows:

By increasing the value of , we explored the interval where . Using these values of and as initial values, we iterated using the following equation:

We used the following equation as a convergence criterion:

Convergence of series solutions

To determine the expansion order of the series solution in the high-order equation, we defined the convergence judgement function as follows:where is the increment of the expansion order and is the minimum positive real solution in the high-order equation (Eq. (18)) when the expansion order is . In many studies, the height of a tree is measured in centimetres[42,43]. Therefore, we used the following judgement function to determine the expansion order:

We increased by until the above condition was satisfied. The smallest expansion order that satisfied Eq. (35) was the expansion order used in this study.

Relationship between critical height and parameter

In this section, we clarify the relationship between the critical height and the density distribution parameter . The critical height ratio was defined as follows:where was the critical height in the variable density model, and was the critical height in the model with constant density. If the specification of each model is the same, the critical height is larger than the critical height of the density constant state when is larger than 1. In contrast, the critical height is smaller than the critical height of the state with constant density when is smaller than 1. Because and are given by Eqs. (22) and (1), respectively, the critical height ratio was obtained as follows:

In order to clarify the relationship between the critical height ratio and the density distribution parameter , we obtained the regression curve (discretely calculated) from . By using this method, the critical height ratio was expressed as a simple function of . By applying it to Eq. (22), the relationship between the critical height and the weight distribution was clarified. In Sect. 2.5, the critical height ratio was defined as follows:where is the critical height when the trees have no branches, which is the same as Eq. (1). The quantity is the critical height when the density distribution parameter is , and the weight ratio is . In other words, Eq. (38) replaces with in Eq. (37).

We assumed the following models as regression models; each parameter was determined by regression analysis using the non-linear least-squares method. Furthermore, the significance level was :

The exponential model and the power-law model in Eqs. (39) and (42) are simple models expressed in monomial form. The linear and polynomial models in Eqs. (40) and (41) are more complicated than the above models, as they are both expressed in a polynomial form. In particular, the polynomial model in Eq. (41) is a quadratic polynomial and is the most complicated among the four models. The power-law model in Eq. (42) cannot be applied to data with . Therefore, we use only the data with to perform regression analysis in the case of the power-law model.

Results and discussion

In this study, based on the physical characteristic values found in[43] and[44], we used an elastic modulus , a standard density , a gravitational acceleration , and a radius .

Convergence judgement

The results of convergence judgement in each model are displayed in Fig. 4. Values for the initial expansion order and the increment . Were used. The values of the convergence judgement function are indicated on the vertical axes in the figure and the values of are indicated on the horizontal axes. In all models, the calculated values converged as the expansion order increased. However, the convergence of Model D was inferior to the other models that had linear density distributions, i.e., Models A, B, and C. The specific convergence judgement function (indicated by the dashed line in each panel) satisfied Eq. (25) in Models A, B, and C. In contrast, this function did not satisfy Eq. (25) in Model D. The specific convergence judgement function (represented by the single-dotted line in each panel) satisfied Eq. (25) for all the models. Therefore, we selected as the expansion order.

Figure 4

Convergence judgement.

Effect of the density distribution on the critical height

The relationship between the critical height and the density distribution parameter is illustrated in Fig. 5. The values of the critical height are indicated on the vertical axes in the figure, and the values of are indicated on the horizontal axes.

Figure 5

Relationship between critical height and parameter

In Model A, the critical height slowly decreased with increasing , whereas the critical height showed a sharp decrease with increasing in Model B. Comparing the two models in more detail, the critical in Model B was higher than that in model A in the interval , whereas the two models showed agreement at . In the interval , the critical height in Model A was higher than that in Model B. The reason for these results was that Model A concentrated more weight at the top of the tree in the interval , whereas Model B concentrated more weight at the top in the interval . This result indicated that increasing the weight of the upper portion of the tree significantly reduced the critical height, as compared to the result for an increase in the weight of the lower portion.

In Models C and D, which have the same total weight regardless of , the critical height displayed a curvilinear increase . Retaining the equal total weights, but with the weight distributed such that the upper and lower parts were lighter and heavier, respectively, the critical height showed a dramatic increase. The difference between Model C (straight line) and Model D (curved line) was negligible when , the value for Model C being slightly higher. The difference between the two models became apparent at higher values of . The values for Model D were higher than those for Model C when .

For all models, when the density was constant (in agreement with Greenhill's model), we obtained a value of , which is almost the same as Greenhill's solution. This indicates that when the governing differential equation is solved by the series solution method, the solutions are as accurate as those obtained by solving exactly if the number of terms is used. The solutions of von Karman and Biot[11] were derived using .

Derivation a simple scaling law by regression analysis

The relationship between the critical height ratio and the density distribution parameter , and the regression curves obtained by non-linear regression analysis are displayed in Fig. 6. The values of the critical height ratio are indicated on the vertical axes in the figure, and the values of are indicated on the horizontal axes. Detailed results of the regression analysis are listed in Table 2. All the parameters of all the regression models () are valid because the p-values are consistently smaller than the significance level .

Figure 6

Critical height ratio and regression curve.

Table 2

Non-linear regression analysis results.

Regression model	Exponential		Linear		Polynomial			Power-law
Parameter	P_1	P_2	P_3	P_4	P_5	P_6	P_7	P_8	P_9
Model A
Estimated value	1.081	− 0.073	− 0.071	1.077	0.010	− 1.002	1.091	0.978	− 0.067
p-value	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	4.2 × 10–15
AIC	− 2.311 × 102		− 2.123 × 102		− 3.407 × 102			− 1.431 × 102
Model B
Estimated value	1.381	− 0.254	− 0.230	1.315	0.109	− 0.556	1.473	0.972	− 0.210
p-value	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	1.1 × 10–11	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16
AIC	− 7.314 × 101		− 5.800 × 101		− 1.078 × 102			− 1.115 × 102
Model C
Estimated value	1.016	0.190	0.191	1.023	0.068	0.191	0.999	1.207	0.085
p-value	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	5.6 × 10–8
AIC	− 2.154 × 102		− 1.912 × 102		− 3.196 × 102			− 7.314 × 101
Model D
Estimated value	1.025	0.243	0.245	1.036	0.110	0.245	0.998	1.282	0.115
p-value	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	2.0 × 10–16	1.3 × 10–7
AIC	− 1.753 × 102		− 1.513 × 102		− 2.620 × 102			− 5.823 × 101

In Model A, the critical height ratio (the ratio of the actual critical height to the value corresponding to the constant-density case) varies from 1.09 to 0.88 (depending on the parameter ). Moreover, all regression models, except for the power-law model, are consistent with the theoretical solution. Based on the Akaike Information Criterion (AIC), the polynomial model is the most accurate model. However, for the purpose of obtaining a simple scaling law for trees, the simplicity of the regression model is very important. Therefore, the exponential model is a good regression model because it is superior in terms of both accuracy and simplicity.

In Model B, the critical height ratio varies from 1.60 to 0.73. The polynomial and power-law models are consistent with the theoretical solution. In contrast, the exponential and linear models are not accurate in the neighbourhood of and . Notably, the power-law model is not valid at From the standpoint of AIC, the most accurate model is the power-law model. This model has good accuracy and simplicity.

In Model C, the critical height ratio varies from 1.25 to 0.87. All regression models, except for the power-law model, are consistent with the theoretical solution. The power-law model breaks down when , and its error with respect to the theoretical solution is larger than that of the other models. Based on the AIC, the most accurate model is the polynomial model. For the purpose of obtaining a simple scaling law for trees, the exponential model is a good regression model because it is superior in both accuracy and simplicity.

In Model D, the critical height ratio varies from 1.38 to 0.85. Regarding the regression models, all models except for the power-law model are consistent with the theoretical solution. The power-law model cannot account for the critical height when , and its error with respect to the theoretical solution is larger than that of the other models. The exponential model shows good accuracy and simplicity.

To derive a simple scaling law such as Greenhill's, we considered the use of simple fractions to express the critical heights. Simple expressions of the regression models considered optimal for each density model in terms of appropriate fractions are shown in Table 3. The limits of application for each formula are also listed in the table.

Table 3

Simplified formulas for critical height.

Model	Formula	Applicable limit
A	l_c=1110e^-3n/40CEγrl2^1/3	0≤n≤3
B	l_c=n^-1/5CEγrl2^1/3	0<n≤3
C	l_c=e^n/5CEγrl2^1/3	-1≤n≤1
D	l_c=e^n/4CEγrl2^1/3	-1≤n≤1

Application to real trees

The results of the application to a real tree, in which the balance of branch and trunk weights and their distribution on the critical height were investigated, are displayed in Fig. 7. The results for Models E and F are shown in panels (a) and (b), respectively.

Figure 7

Relationship between the critical height and the branch distribution (contained in the weight ratio).

In both models, the critical height increased with increasing . Focusing on , it was possible to increase the critical height by making the top lighter, even when the same number of branches and leaves were used. Comparing Models E and F, the quadratic distribution was more favourable when the top was lighter, and the linear distribution was more favourable when the top was heavier.

Based on King et al.’s study of the weight ratio[40] and Hirata’s study of the weight distribution of trees[41], the value ranges we used for and for real trees were: and . We found the following maximum and minimum values of in these ranges:

The physical meaning of the quantities and in Eqs. (43) and (44) may be expressed as "the maximum and minimum values of the critical height ratio with respect to the case when all branches and leaves are cut off". Even if a real tree has the most unfavourable configuration (), the mechanical critical height is reduced by only about 15%.

Niklas reported that the critical height obtained by the Greenhill equation has a safety factor of approximately () with respect to the actual height of the tree, [45]. Combining the results of Eqs. (43) and (44) with the results of Niklas et al., we can estimate that the safety factor for the critical height equation, considering the weight distribution of trees and the weight ratio of the trunk to the branches, lies approximately in the following range:

Actual trees have highly diverse weight distributions. However, by performing calculations using a density distribution function that can represent such a wide range of weight distributions and then refining the results using actual measurement data, it is possible to estimate the critical height more accurately than using the critical height equation which ignores existing branch and leaf weights.

Buckling modes under self-weight

The buckling modes under self-weight for all the density models are represented in Fig. 8. The vertical axis in each panel represents the dimensionless coordinate , and the horizontal axis in each panel represents the dimensionless deflection . The buckling modes for are shown for models A and B. The buckling modes for are shown for models C, D, E and F.

Figure 8

Buckling modes.

There were no significant differences in the buckling modes for any of the density models. The results are similar to those of the first-order mode in Euler buckling, where the deflection uniformly increases from the fixed end to the free end. No significant difference in the mode shapes was observed when the density distribution parameters and the weight ratios were changed.

Conclusions

In this study, to clarify the effect of the weight distribution of trees on their critical height, the critical height for self-weight buckling was formulated for cylindrical models with various weight distributions and the critical height equations that include the influence of branch weight was derived for the first time. Regression analysis was performed on the discrete theoretical solutions, and a simple relationship between the weight distribution and the limit height was derived for each density model. Furthermore, using the obtained critical heights, the mode shapes under buckling under self-weight for various weight distributions were obtained. As a result, the following new findings were obtained:

When a certain number of branches and leaves are distributed in the height direction, the risk of buckling due to self-weight can be reduced if the upper part of the trees is lighter and the lower part is heavier. This method of weight distribution can increase the critical height by a maximum of approximately 1.25 times for a linear distribution and by a maximum of approximately 1.38 times for a curvilinear distribution, compared with a constant weight distribution in the vertical direction. In addition, the total allowable weight is approximately 2 times larger in the case of a linear distribution and approximately 2.6 times larger in the case of a curved distribution.

The buckling mode when self-weight buckling occurs is similar to the first buckling mode in the case of general long columns, regardless of the weight distribution in the height direction. In addition, the mode shape hardly changes even if the distribution form changes in the same density model.

Based on the weight distributions of real trees and the measurements of branch weights, it can be said that trees distribute their branches and leaves in a very rational way that has little effect on the critical height. Even if we consider the most unfavourable condition of a real tree (), the critical height is only reduced by approximately 15% from the condition with no branches or leaves at all. Based on the findings of Niklas[45], it can be estimated that the safety factor of the actual trees is approximately 3.4 to 3.9 in relation to the theoretical critical height of the present study, considering the branches. The method in this study can obtain more accurate critical height than previous methods that ignore branch weight because a safety factor of approximately 4 was reported in previous research[45].

In future work, based on the methods and findings of this study, the shape law and rationality of plant morphology will be further explored from a theoretical point of view. Currently, we plan to investigate the effects of crown weight and position using more complex density function models and to extend the method to hollow cylindrical models such as bamboo and to solid and hollow models with a taper. In the future, we intend to summarise these findings and propose the next generation of material-saving and high-performance structural designs and rational new design methods.