A Numerical Approach to Estimate Natural Frequency of Trees with Variable Properties

Abstract

Free vibration analysis of a Euler-Bernoulli tapered column was conducted using the finite element method to identify the vibration modes of an equivalent tree structure under a specified set of conditions. A non-prismatic elastic circular column of height L was analysed, taking distributed self-weight into account. Various scenarios were considered: column taper, base fixity, radial and longitudinal stiffness (E) and density (ρ) and crown mass. The effect on the first natural frequency was assessed in each case. Validation against closed form solutions of benchmark problems was conducted satisfactorily. The results show that column taper, base fixity and E/ρ ratio are particularly important for this problem. Comparison of predictions with field observations of natural sway frequency for almost 700 coniferous and broadleaved trees from the published literature showed that the model worked well for coniferous trees, but less well for broadleaved trees with their more complicated crown architecture. Overall, the current study provides an in-depth numerical investigation of material properties, geometric properties and boundary conditions to create further understanding of vibration behaviour in trees.

1. Introduction

Applied wind loads have a significant impact on the ecology, physiology and morphology of land plants [1,2]. Responses of plants to applied wind loads range from minor movement of leaves, branches and stems through to catastrophic failure [3]. In forests (as well as in other crop plant systems), the catastrophic failure of single trees or groups of trees is a source of economic loss and can pose a threat to human safety [1,2].

Previous studies have investigated how wind affects tree growth and development, the mechanics of wind loading on trees and the risk of tree failure [1,2,3,4,5,6]. These studies have shown that the size, shape and internal wood properties of a tree are the result of the applied forces acting on it, and that these developmental responses improve the ability of the tree to resist these applied loads [2,5,6,7]. Several authors have developed models based on engineering principles to investigate the mechanical response of trees to wind loading [8,9]. These models can be used to help investigate the mechanics of tree failure or to provide the biomechanical inputs into investigations of tree developmental responses to applied loads [2,6]. When modelling the interaction between wind and trees, it is important to consider trees as dynamic structures and determine the important parameters, such as the fundamental frequency and damping, that affect their response [10,11,12]. The dynamic behaviour of trees can be affected by their geometric (e.g., mass distribution and cross-section along tree structures) and material (e.g., density, viscous damping and modulus of elasticity) properties. It is also affected by the characteristics of their root systems and the properties of the soil they are growing in [13]. Several previous studies have focused on the effects of these properties on tree dynamic behaviour, either experimentally [14,15,16,17,18] or theoretically [19,20,21,22,23,24,25,26,27].

Several common approaches to estimate the fundamental (lowest) frequency of trees are presented in the literature. One approach is to assume that a tree can be approximated by a single degree of freedom system, so that its fundamental angular frequency, ${\omega }_n,$ is equal to the square root of the ratio of its stiffness (k) to mass (m) [21].

$${\omega }_n=\sqrt{\frac{k}{m}}$$

(1)

Milne [28] presented a method for estimating a whole tree’s fundamental angular frequency which was based on the application of the principle of conservation of energy and sub-dividing it into sections:

$${\omega }_n=\sqrt{\frac{{\displaystyle \sum }0.5{\eta }_i{\theta }_i}{{\displaystyle \sum }0.5x_i^2{m}_i}}$$

(2)

where the subscript i refers to a specific vertical section of the tree stem, x_i is the maximum horizontal displacement, m_i is the mass, η_i is the bending moment in the stem and θ_i is the angular displacement of the stem for the ith vertical section.

Another method proposed was derived from the dynamics of distributed property systems [29]. For this method, it is assumed that the modulus of elasticity and the mass are uniformly distributed, and the cross-section is constant throughout the length of the tree stem.

$${\omega }_n={\alpha }_n^2\sqrt{\frac{EI}{m{L}^4}}\&{\alpha }_n=1.875,4.694,\dots$$

(3)

where E is the modulus of elasticity, I is the second moment of area (or moment of inertia), m is uniformly distributed mass along the length (L) of the stem and α_n are the solutions of the beam equation and depend on the vibration mode. However, the cross-section of a tree stem is not constant along its length, so it is more realistic to represent it as a tapered cantilever beam. The following equation developed by Blevins [30] can be used to find the natural frequency (f_n) of slender and tapered cantilever beam elements:

$$f_n=\frac{{\lambda }^2}{2\pi L^2}\sqrt{\frac{E{I}_0}{\rho A_0}}=\frac{{\omega }_n}{2\pi }$$

(4)

where I₀ is the basal area moment of inertia, A₀ is the area of the cross-section at the base of the cantilever beam (stem), L is the length of the tree stem, E is the modulus of elasticity and ρ is the material density. The unitless parameter, λ, changes due to variations in the physical properties of the tree (such as shape, mass distribution, vertical or horizontal orientation, type of basal support, mode of bending, taper of the beam and so on). Values for trees can be estimated through targeted data collection for a specific species. For example, Milne [28] estimated the average value of λ to be 2.08 (standard deviation = 0.06) for Sitka spruce (Picea sitchensis (Bong.) Carr.) trees.

There have been several studies that have collected data on the natural frequency of trees in field studies and then fitted regression equations to these data. For example, Mayhead [31] fitted a number of models to the data collected during the 1960s by the British Forestry Commission. He found that sway period (T) was best predicted using the following equation:

$$T=\frac{1}{f_n}=0.86+0.74\frac{L\sqrt{ML}}{DB{H}^2}$$

(5)

where M is the total mass of the tree, L is the height of the tree and DBH is the diameter at breast height (approximately 1.4 m). More recently, Moore and Maguire [11] combined the results of several previous studies and re-analysed the data for the relationship between tree height, tree diameter and tree natural frequency across eight different conifer species (a total of 602 trees). They presented the following empirical formula based on regression analysis:

$$f_n=0.0948+3.4317\frac{DBH}{L^2}-0.7765I_P\frac{DBH}{L^2}$$

(6)

where DBH is the diameter at breast height (cm), L is stem length (m) and I_p is an indicator variable which has a value of 1.0 if the genus is Pinus and 0.0 otherwise. A more recent and comprehensive analysis of tree natural frequency data by Jackson et al. [32] shows that the fundamental natural frequency of conifer trees and broadleaved species with monopodial architecture can be predicted using the equation for the vibration of a tapered cantilever beam. However, for some broadleaved species with more complex architecture (i.e., a more decurrent form), this approximation to a tapered cantilever beam no longer holds.

According to Equation (1), to estimate the modal frequencies of trees it is important to focus on the variation in mass and stiffness matrices. The shape and internal wood properties of trees reflect their life histories and changing biomechanical requirements over time [33,34]. Many wood science, tree mechanics and allometry studies have focused on the variation in wood stiffness (E) and density (ρ) within trees. For example, longitudinal gradients in both these properties were found for black locust (Robinia pseudoacacia L.) by Niklas [35]. This was attributed to the transition from sapwood to heartwood along the tree with age. Similar differences in physical properties have been observed along the stems of a range of species [36,37,38]. However, no single pattern of variation in these properties has been found. For example, in radiata pine (Pinus radiata D Don), an increase in E from the stem base to approximately 20–30% of the stem height was found, followed by a decrease towards the top of the stem [39,40]. Similar trends of increasing E with height up the stem over the lowest portions of the stem were observed in Cryptomeria japonica (L.f.) D. Don [41] and Picea abies (L.) H. Karst [42]. In black locust, Niklas [35] found reduced wood density near the centre of the stem with an increase in wood density with increasing radial distance from the pith as there was a transition from heartwood to sapwood. He also found an approximately parabolic variation in radial stiffness and E/ρ variation across the stem diameter. Similar non-linear radial variation in wood density and stiffness have been found for many, mainly conifer, species [43,44,45]. In broad terms, there appear to be two main strategies in trees for radial stiffness: (1) a stiff ‘crust’, with the outer layers stiffer than the inner, and (2) a more rigid inner core [35].

Despite an awareness of the importance that variations in mass and stiffness within a tree have on its dynamic behaviour, few if any quantitative analyses have explicitly considered these. In this paper, we examine the effects of radial and longitudinal variations in stiffness and density on the fundamental vibration frequency of a tree structure. We also examine the effect of taper, base fixity and crown mass on the vibration modes. This work builds on the recently published study by Dargahi et al. [46], which focused on the effects of material property variation, tree taper and base fixity on the buckling behaviour of trees under self-weight loading.

2. Methods

2.1. Theoretical Basis for Vibration Analysis

A vibration analysis generally follows four steps. First, the structure or system of interest is identified, its boundary conditions are estimated and its interfaces with other systems are identified. Second, the natural frequencies and mode shapes of the structure are determined by analysis or direct experimental measurement. Third, the time-dependent loads on the structure are estimated. Fourth, these loads are applied to an analytical model of the structure to determine its dynamic response. The crucial steps in the vibration analysis are the identification of the structure and the determination of its natural frequencies and mode shapes [30].

Any structure with mass and elasticity will possess one or more natural frequencies of vibration. The natural frequencies are the result of cyclic exchange of kinetic and potential energy within the structure. The kinetic energy is associated with velocity of the structural mass, while the potential energy is associated with storage of energy in the elastic deformations of a resilient structure. The rate of energy exchange between the potential and kinetic forms of energy is the natural frequency [30].

Figure 1 shows a straight-sided cantilever column without damping subjected to an external variable load, P. The governing equation for transverse vibration of this column can be written as:

$${\displaystyle \sum }{F}_y=0\Rightarrow \frac{\partial \nu \left(x\right)}{\partial x}=P\left(x,t\right)-m\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}$$

(7)

where y is the lateral direction, x is the elevation from the base, u(x,t) is the motion of the column in y direction, m(x) is mass per unit length and P(x,t) is the external transverse force.

We know that:

$$\nu \left(x,t\right)=\frac{\partial M\left(x,t\right)}{\partial x};M\left(x,t\right)=D\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2};m\left(x\right)=\rho \left(x\right)A\left(x\right)$$

(8)

where ν(x,t) is transverse shear, M(x,t) is the bending moment, ρ(x) is density and A(x) is the cross-sectional area. Substituting Equation (8) into Equation (7) gives:

$$\rho \left(x\right)A\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}+\frac{{\partial }^2}{\partial x^2}\left[D\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]=P\left(x,t\right)$$

(9)

where D(x) = E(x)I(x) is flexural rigidity, E(x) is Young’s modulus and I(x) is the area moment of inertia.

This is the general governing differential equation for vibration of a non-prismatic Bernoulli beam (column) under variable transverse force without considering damping. Setting P(x,t) = 0, the free vibration equation is obtained.

In the above system of equations, the deflection due to the shear stress in the column was ignored. The analysis of beam (column) vibration including both the effects of rotational inertia and shear deformation is called Timoshenko beam theory, which is outside the scope of this paper.

Considering the motion represented by harmonic vibration, the transverse motion of the column is obtained using the following relation:

$$u\left(x,t\right)=\phi \left(x\right)q\left(t\right)=\phi \left(x\right).e^{i\omega t}$$

(10)

where ϕ(x) and ω are the mode shape function and angular natural frequency of the beam, respectively. Substitution of Equation (10) into Equation (9) leads to:

$$\frac{d^2}{d{x}^2}\left[D\left(x\right)\frac{d^2\phi \left(x\right)}{\partial x^2}\right]-{\omega }^2\rho \left(x\right)A\left(x\right)\phi \left(x\right)=0,0\le x\le L$$

(11)

For simplicity, for a column with a non-prismatic circular cross-section, we can define a non-dimensional coordinate ξ = x/L, then:

$$I\left(\xi \right)=I_0{(1-\alpha \xi )}^4$$

(12)

where α is the taper ratio [(D₀-D_t)/D₀] and D₀ and D_t are the diameters of the tree at the base and top of the tree stem, respectively.

In general, we can rewrite the flexural rigidity and distributed mass as:

$$D\left(\xi \right)=E\left(\xi \right)I\left(\xi \right)=E_0{I}_{0}g\left(\xi \right)m\left(\xi \right)=\rho \left(\xi \right)A\left(\xi \right)={\rho }_0{A}_{0}f\left(\xi \right)$$

(13)

where I₀, E₀, ρ₀ and A₀ are the area moment of inertia, Young’s modulus, density and cross-section area at the base of the stem, respectively.

Consequently, Equation (11) can be transformed into a normalised form:

$$\frac{d^2}{d{\xi }^2}\left[g\left(\xi \right)\frac{d^2\phi \left(x\right)}{\partial {\xi }^2}\right]-\left(\frac{{\omega }^2{H}^4{\rho }_0{A}_0}{E_0{I}_0}\right)f\left(\xi \right)=0,0\le \xi \le 1$$

(14)

Therefore, we can define natural frequency parameter (λ) as:

$$\frac{{\omega }^2{H}^4{\rho }_0{A}_0}{E_0{I}_0}={\lambda }^4\Rightarrow \omega =\frac{{\lambda }^2}{L^2}\sqrt{\frac{E_0{I}_0}{{\rho }_0{A}_0}}orf=\frac{{\lambda }^2}{2\pi L^2}\sqrt{\frac{E_0{I}_0}{{\rho }_0{A}_0}}$$

(15)

The column base condition is bounded by two cases: clamped (fixed) and pinned. In general, we can consider intermediate cases using a rotational spring restraint at the base of the column. The deflection, u, is therefore related to the rotation angle, θ, by:

$$\theta =\frac{\partial u\left(x,t\right)}{\partial x}=\frac{\partial u\left(\xi ,t\right)}{\partial \xi }$$

(16)

The boundary conditions applied to this column at the spring-supported end are:

$$\begin{gathered} u\left(0,t\right)=0 \\ EI\frac{{\partial }^{2}u\left(0,t\right)}{\partial x^2}=k\frac{\partial u\left(0,t\right)}{\partial x} \end{gathered}$$

(17)

where K is the rotational spring stiffness and the spring is assumed to be massless.

2.2. Calculation of Natural Frequency Using Numerical Analysis

The first natural frequency parameter (λ) of a column was found with the generalised eigenfrequency approach using the finite element analysis method. The commercial software SAP2000 (Computers and Structures Inc, Berkeley, CA, USA) and COMSOL Multiphysics (COMSOL Inc, Burlington, MA, USA) were used to investigate the effect of the different parameters on the natural frequency of non-prismatic circular columns. A free end–fixed base case was modelled, with rigid support for torsional moments and lateral translations. The column was discretised (sub-divided) into finite elements, with the properties varying linearly across the elements. The effect of variations in selected geometric and material properties on the vibration behaviour was investigated: column taper, end fixity, radial and longitudinal stiffness/density variation and relative crown mass and location.

2.2.1. Effect of Base Fixity

We did not attempt to model the behaviour of the whole root–soil system explicitly. Since the effect of the mass of the root–soil system on the overall behaviour has been found to be less significant with higher centres of mass of the above-ground components [18], this has been ignored. In addition, increases in the rooting depth will tend to increase the stiffness of the root–soil volume, counteracting any changes in the root–soil mass [47]. Hence, we have only included a basal rotational spring stiffness term (K) in the numerical models (assumed to be identical in the two orthogonal horizontal directions) to test the effect of the overall root–soil system rigidity on tree vibration behaviour. Each rotational stiffness was normalised by the stiffness magnitude to create a non-dimensional spring stiffness (K’) and the results are presented herein using the non-dimensional rotational stiffness K’ = KL/E₀I₀. The value of K’ will approach infinity for an equivalent clamped (fixed) case (θ$\to$0) and zero for a pinned condition (where M (base bending moment)$\to$ 0). Validation of the models was conducted by comparison with the previous analyses of Blevins [30].

2.2.2. Effect of Longitudinal Stiffness Variation

To model the longitudinal stiffness gradient variations for trees described in the literature, three cases (see Figure 2a, (i) to (iii)) were chosen for these analyses: (i) elastic modulus decreasing linearly from the bottom (E_b) to the top of the stem (E_t), with E_t/E_b = 0, 0.25, 0.5, 0.75 and 1, (ii) elastic modulus increasing linearly from the bottom (E_b) to one-third of the stem height (E_m), followed by a linear decrease to the top of the stem (Et), with ϕ = 1/3, E_t/E_b = 0.25 and E_b/E_m = 0.25, 0.35, 0.5, 0.75 and 1 and (iii) E_t/E_b = 0.25 and E_b/E_m = 0.5 and varying the position of the gradient change (point m) with L/3, L/2, 2L/3 and 5L/6. In each case, the density (ρ) has been held constant, thus the value of E/ρ changes in a similar manner to that described in Equation (18).

$$\begin{gathered} E=E_b+\left(\frac{x}{l}\right)\left(E_t-E_b\right),0\le x\le lCase\left(i\right) \\ \{\begin{array}{c}E=E_b+\left(\frac{x}{\varphi .l}\right)\left(E_m-E_b\right),0\le x\le \mathsf{\varphi }.l\\ E=E_m+\left(\frac{x-\varphi .l}{l-\varphi .l}\right)\left(E_t-E_m\right),\mathsf{\varphi }.l\le x\le l\end{array}Case\left(ii\right)and\left(iii\right) \end{gathered}$$

(18)

2.2.3. Effect of Radial Stiffness Variation

To consider the effect of radial stiffness (E) variation on the vibration behaviour, a composite section was defined with two concentric rings (see Figure 2b). Different stiffness moduli were assigned to each ring to represent the two different tree strategies (stiff core and stiff outer ring) outlined in the literature. The modelling herein assumes the presence of two rings of different radii, with the same area to represent the outer-wood and core-wood of the tree (i.e., A_o = A_c = 50% of the total area of the base cross-section). The stiffness of each ring was then varied according to the ratios: E_c/E_o of 10⁻³, 0.1, 0.25, 0.5, 1 and 1.5. The same density was assumed for both the outer-wood and core-wood. To calculate the values of first natural frequency parameter (λ), outer-wood properties at the base were used.

2.2.4. Effect of Longitudinal Density Variation

To model the longitudinal density gradients reported in the literature, three cases (see Figure 2c, (i) to (iii)) were chosen for analysis: (i) density varying linearly from the bottom (ρ_b) to the top of the stem (ρ_t), with ρ_t/ρ_b = 0, 0.25, 0.5, 0.75,1 and 1.25, (ii) density decreasing linearly from the bottom (ρ_b) to half of the stem height (ρ_m), followed by a linear increase to the top of the stem (ρ_t), with ϕ = 1/2, ρ_t/ρ_b = 0. 5 and ρ_m/ρ_b = 0.25, 0.5, 0.75, 1 and 1.25 and (iii) ρ_t/ρ_b = 0.5 and ρ_m/ρ_b = 0.25 and varying the position of the gradient change (point m) with L/3, L/2, 2L/3 and 5L/6. In each case, the stiffness (E) has been held constant, as shown in Equation (19).

$$\begin{gathered} \rho ={\rho }_b+\left(\frac{x}{L}\right)\left({\rho }_t-{\rho }_b\right),0\le x\le LCase\left(i\right) \\ \{\begin{array}{c}\rho ={\rho }_b+\left(\frac{x}{\varphi .L}\right)\left({\rho }_m-{\rho }_b\right),0\le x\le \varphi .L\\ \rho ={\rho }_m+\left(\frac{x-\varphi .L}{L-\varphi .L}\right)\left({\rho }_t-{\rho }_m\right),\varphi .L\le x\le L\end{array}Case\left(ii\right)and\left(iii\right) \end{gathered}$$

(19)

2.2.5. Effect of Radial Density Variation

A similar approach to that used for radial stiffness variation was used to model the effect of radial density (ρ) variation on vibration behaviour. A composite section was defined with two concentric rings (see Figure 2b). Different densities were assigned to each ring to explore the two different strategies (light core and heavy outer ring) proposed in the literature. Two rings of different radii were defined, with the same area to represent the core-wood and outer-wood of the tree (i.e., A_o = A_c = 50% of the total area of the base cross-section). The density for each ring was then varied according to the ratios: ρ_c/ρ_o of 0, 0.25, 0.5, 1 and 1.25. The same elastic modulus was assumed for both the core-wood and outer-wood. To calculate the values of first natural frequency parameter (λ), outer-wood properties at the base were used.

2.2.6. Effect of Crown Mass

The effect of crown mass on the natural frequency was investigated by creating an equivalent “crown” using an additional mass (W_c) at different locations along the stem (see Figure 3). (i) In the first case, the additional crown mass was applied to the top of the stem, (ii) in the second, the same crown-to-stem mass ratios were applied to the upper two-thirds (by volume) of the stem and (iii) in the third case, they were applied to the upper half (by volume) of the stem. In all cases, we define different crown-to-stem mass ratios (W_c/W_s) at these particular locations.

2.3. Estimating Natural Frequency for Coniferous and Broadleaved Trees

Fundamental vibrational frequencies (f_n) were collated from the literature on field studies of 703 broadleaf and conifer trees growing at different sites and under different conditions, spanning height and DBH ranges of 5.9–35.4 m and 6.9–82.8 cm, respectively [11,12,15,28,31,48,49,50,51,52,53]. These field data were used to validate the tapered cantilever beam model and provide further fundamental understanding of the sway behaviour of conifer and broadleaf trees. Results from these studies were used to estimate the natural frequencies and damping ratios of trees with the presence or absence of foliage, and for full and de-branched crowns, enabling the effects of the crown architecture to also be assessed. The field sway data covers thirteen different species: (a) conifers—Sitka spruce (Picea sitchensis), Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco), Norway spruce (Picea abies), white spruce (Picea glauca (Moench) Voss.), Lodgepole pine (Pinus contorta Dougl.), Scots pine (Pinus sylvestris L.), Corsican pine (Pinus nigra Arnold) and red pine (Pinus resinosa Ait.), and (b) broadleaves—red maple (Acer rubrum L.), sugar maple (Acer saccharum Marsh.), shagbark hickory (Carya ovata (Mill) K. Koch), red oak (Quercus rubra L.) and lime (Tilia europaea L.). Data on the stiffness (E), density (ρ), taper (α) and crown/stem mass ratio (W_c/W_s) were found for each species in the literature and are summarised in

Table 1. Average values of wood properties (modulus of elasticity of green wood, E_g, and density ρ), tree taper (α) and the ratio of crown mass to stem mass (Wc/Ws) for each species used in the analysis to estimate fundamental natural sway frequency. The number of trees (n) for which sway data were obtained from the literature is also given.

Code	Tree Species	n	Eg (GPa)	ρ (kg/m3)	Eg/ρ	Taper (α)	Wc/Ws
CP	Corsican pine (Pinus nigra)	57	8.70	657	0.013	0.85	0.34
DF	Douglas-fir (Pseudotsuga menziesii)	17	9.83	583	0.017	0.77	0.16
LP	Lodgepole pine (Pinus contorta)	40	6.90	487	0.014	0.83	0.33
NS	Norway spruce (Picea abies)	32	6.23	598	0.010	0.90	0.32
RP	Red pine (Pinus resinosa)	300	8.80	410	0.021	0.81	0.22
SP	Scots pine (Pinus sylvestris)	20	7.33	700	0.010	0.90	0.29
SS	Sitka spruce (Picea sitchensis)	175	7.53	447	0.017	0.84	0.50
WS	White spruce (Picea glauca)	6	7.40	466	0.016	0.91	0.34
	Averages		7.77	551	0.015	0.85	0.31
LT	Lime (Tilia europaea)	18	11.7	530	0.022	0.84	0.18
AR	Red maple (Acer rubrum)	7	9.6	524	0.018	0.67	0.22
QR	Red oak (Quercus rubra)	11	9.9	665	0.015	0.75	0.32
CO	Shagbark hickory (Carya ovata)	5	10.8	649	0.017	0.76	0.39
SM	Sugar maple (Acer saccharum)	15	10.7	560	0.019	0.79	0.18
	Averages		10.54	585.6	0.018	0.76	0.26

(see the Electronic Supplementary Material, Tables S1 and S2). While these properties can all vary substantially within a species due to genetic and environmental factors, we only had species-level average values as data for the individual trees were not available. Instead, we assumed that these properties were normally distributed with a standard deviation that was equal to 5% of the mean value. We then ran 100 simulations of the model to test the sensitivity to variations in these properties. The data were split into four groups for analysis: (1) conifers with full crowns, (2) broadleaves with full crowns, (3) debranched conifers and (4) leafless broadleaves. The natural frequency parameter (λ) was assumed to be 1.5, except for the case of conifers with their branches removed, when it was assumed to have a value of 2.5. Ordinary least squared regression was used to investigate the relationship between tree fundamental frequency and DBH/L². This relationship was compared among the four analysis groups described above. Due to the limited overlap in values of DBH/L² for several of the broadleaved species, it was not possible to test whether this relationship differed among species. The level of agreement between the predictions from the tapered cantilever beam model and the field observations of fundamental natural frequency for these groups was assessed using Lin’s concordance correlation coefficient [54,55].

3. Results

3.1. Effect of Base Fixity

The relationship between the first natural frequency parameter, λ, and the non-dimensional rotational stiffness, K’, for tapered (cylindrical) columns is found to be non-linear (Figure 4a). There is very little increase in natural frequency for large changes in rotational stiffness above a value of K’ = 10. However, once the non-dimensional rotational stiffness (K’) falls below 10, the natural frequency begins to drop rapidly. Comparison between our numerical modelling and Blevins’ solutions show differences of less than 1%. For columns with different rates of taper (α), the variation in first natural frequency parameter, λ, with the non-dimensional rotational stiffness, K’, is shown in Figure 4b. The non-dimensional rotational stiffness was varied from K’ = 0.01 (equivalent pinned condition) to 100 (approaching a clamped condition). For more tapered stems, basal rotational restraint had a greater effect on the natural frequency.

3.2. Effect of Longitudinal Stiffness Variation

For cylindrical columns with constant E, the first natural frequency parameter, λ, increased with increasing column taper ratio (α) (Figure 5a–c). With a linear decrease in E from the bottom to the top of the stem, the natural frequency parameter decreased compared with the case of constant E. The effect of this change in stiffness was relatively small for cylindrical columns (up to 7%), but more significant for highly tapered columns (up to 12%) (Figure 4a). This is due to the more rapid changes in stiffness compared to the sectional area. Under the assumption of a peak in stiffness at half stem height, the natural frequency increased relative to the uniform E case (up to 17%), with the greatest increase occurring for the more tapered stems (Figure 4b). Under the case of the varying changepoint height for E, the largest increases in natural frequency occurred at the lowest change point and for the most highly tapered stems. However, the magnitude of increase was very similar for all changepoint heights (up to 10%).

3.3. Effect of Radial Stiffness Variation

The value of the natural frequency parameter decreased with reducing stiffness ratio (i.e., outer-wood E > core-wood E), up to 7% for both cylindrical and tapered columns (Figure 6). According to the model results, the elastic ratio of 10⁻³ is a lower bound for the behaviour. Beyond this point, the behaviour of the stem approaches that of a hollow section. For the stiffness ratio of E_c/E_o = 1.5, an increase in the natural frequency of up to 3% was found.

3.4. Effect of Longitudinal Density Variation

For a cylindrical column with no taper, a reduction in density from the base to the top resulted in an increase in λ (up to 50%) (Figure 7a). The effect of this reduction in density on λ was less significant for highly tapered columns (up to 40%). This is due to the less rapid changes in mass compared to the sectional area. For the case of a change in density at the midpoint height of the column, the natural frequency increases (up to 23%), and this change was more evident for the highly tapered columns (Figure 7b). Even when the density at the midpoint height was greater than at the base, λ still increased as we assumed that density at the top of the column was 50% of the value at the base. For the case where the height of the minimum value varied up the column, the natural frequency increased, and the largest increase occurred with the highest change point (Figure 7c). However, the magnitude of the increase was very similar for all the changepoint heights (up to 27%).

3.5. Effect of Radial Density Variation

The value of the natural frequency parameter decreased with reducing ratio of core-wood density to outer-wood density ratio, for both cylindrical and tapered columns (Figure 8). According to the model results, the density ratio of 0 is an upper bound for the behaviour. This condition approaches that of a weightless core section.

3.6. Effect of Crown Mass

The natural frequency reduces as the relative crown weight increases (by up to 60%) (Figure 9). When there is no additional crown mass (i.e., W_c/W_s = 0), taper has a noticeable effect on the natural frequency parameter. However, as the ratio of crown mass to stem mass increases, the natural frequency becomes invariant to changes in taper, particularly when the crown mass is added at 2L/3 and L/2 (Figure 9b,c). This suggests that the crown weight is dominating the taper effect to some extent. However, adding a crown mass to the very top of a tapered column reduces the natural frequency while increased taper ratio in the absence of a crown mass increases natural frequency (Figure 9a).

3.7. Estimating Natural Frequency for Coniferous and Broadleaved Trees

There was a strong relationship between tree natural frequency and the ratio of diameter at breast height to total tree height squared (Figure 10). However, the slope of this relationship differed significantly among the four different categories (p < 0.0001). The slope was lowest for broadleaved trees with full crown and was steepest for conifers with their branches removed. There was some indication that the relationship may be non-linear, especially for the broadleaved trees with full crown; however, this is mostly due to the dataset for lime trees which has several high leverage points (i.e., trees with high values of DBH/L² and comparatively low natural frequency). The fundamental natural frequencies predicted by the model proposed in this paper were in general agreement with field observations for conifers with full crown (concordance correlation coefficient = 0.79) (Figure 11). The model overpredicted the natural frequencies of coniferous trees with their branches removed and had much poorer performance for broadleaved trees, particularly those with full crown (concordance correlation coefficient < 0.3). This latter result was mostly driven by the dataset for lime trees where there was a weak relationship between tree natural frequency and DBH/L². With the exception of these data, all the simulations produced results that were in general agreement with field observations.

4. Discussion

As well as the stem geometry and material properties, the crown and branch architecture has been shown to influence the sway frequency, damping and dynamic amplification factor of trees under dynamic loading [15,22,25,32]. The crown shape of trees can vary according to the species and is due to a combination of the inherent form and their growth environment. McCurdy el al. [56] identified six distinct crown shapes: oblong, round, oval, vase, pyramidal and weeping. In forest settings, crown shape is often influenced by tree position relative to other trees in the canopy (although the basic shape is characteristic for all species). The primary factor determining these characteristic crown shapes is apical (epinastic) dominance (i.e., the upward growth of a leader at the expense of the lateral shoots) [57]. Tree species with strong apical dominance usually have height to width dominance and a single dominating trunk and leader. Examples of tree species of this type are sweet gum (Liquidambar styraciflua L.), Yellow tulip poplar (Liriodendron tulipifera L.), cottonwoods (Populus section Aigeiros), alders (Alnus spp.) and most conifers. These are termed excurrent crowns and favour oblong and pyramidal shapes [58]. In contrast, decurrent trees with weak apical dominance have crown widths that grow nearly as fast as their height (especially when open-grown) and this creates large spreading crowns. These species tend to have no predominant leader and can have similar multiple forks in the central trunk. These trees have crowns that tend to favour round or oval shapes and typical species are broadleaf trees such as oaks (Quercus spp.) and maples (Acer spp.) [58].

The tapered cantilever beam sway model proposed in this study adequately predicts the natural frequency for conifers, but not for broadleaves. As other researchers have found, excurrent trees have natural frequencies which are proportional to the product of DBH/L² [11,32]. Interestingly, the presence of a crown seems to reduce the natural frequency and improves the fit of the model compared with the situation where the model is applied to a stem with the crown removed. Broadleaves tend to have a more complex geometry, with multiple close natural frequencies and are therefore more dependent on the branch/crown architecture and presence of leaves [17,32,53]. For the broadleaves with crowns, the model fit is generally quite poor, and as would be anticipated, this improves with the removal of the leaf mass. As occurs with excurrent trees, the natural frequency of decurrent trees also reduces with the removal of the leaf mass from the crown, although this change is relatively small compared to that for excurrent trees. These findings are supported by field studies [8,11,15,28], where pruning and seasonal changes have more profound effects on excurrent trees. The change in natural frequency with pruning have mostly been associated with mass change, but numerical studies by Kerzenmacher and Gardiner [19] showed that mass change alone could not completely account for the change in natural frequency, as trees went from a fully-branched state to a de-branched state. They hypothesized that branch swaying also had an effect on overall dynamic behaviour of trees.

In decurrent trees, the branch sway has much greater importance with relation to the natural frequency [53], leading to the divergence from the proposed mechanical sway model. The distribution of mass in conifers is also quite different, with a greater proportion of mass in the upper half of the crown, compared to broadleaves with a more distributed branch mass. Broadleaves may also have what is known as multiple resonance damping [17,24]. Transfer of mechanical energy can occur due to the overlap of the resonance spectra of the large branch structures, which increases the structural damping. The frequency of these large branches can be close to the natural frequency of the entire tree and the action of the primary and higher-order branches therefore combine to produce the overall sway response [59,60,61]. The energy exchange between these different branches is related to their amplitudes, frequencies and the phase relations, and this creates a coupled damped oscillation system with characteristics very different from the proposed tapered cantilever beam model. It may be feasible that the simpler pendulum response and crown architecture of excurrent trees may make them more vulnerable to dynamic wind damage and higher wind-induced forces.

Previous studies modelling tree dynamic behaviour have generally assumed uniform wood properties in both the radial and longitudinal directions [19,25,50]. However, this ignores the typical patterns that are observed in practice. In many wood science studies, radial and longitudinal patterns of wood density and stiffness are measured, typically at a standard moisture content [33,43,44,45]. In standing trees, these gradients may differ due to moisture content differences between sapwood and heartwood. In our simulations, we accounted for both moisture content differences and gradients in wood properties through varying the ratio of core-wood and outer-wood properties and incorporating different longitudinal gradients. These encompassed the full range of patterns that are expected to exist in nature. Our results showed that varying the internal wood properties’ distributions has a significant effect on tree dynamic behaviour. As expected, placing the stiffer material on the periphery of the stem increases natural frequency as this material is located further from the neutral axis of the stem. While the application of the model to predict natural frequencies of conifers and broadleaves trees used published values of wood density and modulus of elasticity from wood science literature, there is large variation in these properties among individuals within a species and across geographical gradients [45,62,63]. Differences in values of density-specific stiffness (E/ρ) may explain some of the differences between observed and predicted fundamental natural frequencies.

The proposed model also enabled the effect of base fixity on natural frequency to be accounted for. Previous studies that have used the equations of motion for a taper cantilever beam as the basis to study tree oscillations have either assumed the tree to be rigidly anchored and adjusted the effective stem stiffness or represented the whole root–stem system as a spring with a given rotational stiffness. This is because there is a significant body of study is available for the vibration analysis of fixed base-free end columns. Much less work has been done on columns with reduced basal rotational restraint. This is a situation encountered in civil engineering for columns connecting to foundations with relatively ‘soft’ bases and electro-mechanics where nanostructures are fixed to elastic substrates. This problem has been modelled by introducing rotational spring restraint in place of the fixed (clamped) column base condition. For trees, this represents a scenario where the rotational restraint of the root plate is relatively low, e.g., from loss of soil, reduced soil strength/stiffness or root damage [64]. The effect of changing the root system stiffness may be approximately assessed by inspection of the classical equation for rotational stiffness (K) of a rigid circular plate located on the surface of an elastic half-space [65]:

$$K=\frac{8G{R}^3}{3\left(1-\nu \right)}$$

(20)

where G is the composite shear modulus of the soil and roots, R is the radius of the root plate and ν is the Poisson’s ratio of the soil and roots. This model is often applied to foundation problems, but in the current case, G and ν are related to the composite behaviour of the soil and roots in the root plate. This suggests that quite significant changes in soil packing (density) or moisture content are required for the rotational root plate stiffness to affect the vibration behaviour, particularly for trees with high taper ratios.

Overall, the modelling approach employed here has enabled the relative effect of variation in tree taper, internal wood properties, crown mass distribution and the base fixity on tree dynamic behaviour to be investigated. It provides an important advancement over earlier models that make much simpler assumptions for these factors and enables the various biomechanical strategies employed by trees to be investigated. Although the proposed model appears to predict the behaviour of conifers adequately, modelling the mechanical behaviour of broadleaved trees can still be undertaken using the finite element method, but will require a more detailed description of the tree topology. New technologies such as terrestrial laser scanning systems provide a means for capturing the data required to generate realistic descriptions of tree architecture and the resulting mass and stiffness distributions needed to predict their dynamic behaviour [32].

5. Conclusions

Finite element-based vibration analysis was used to investigate the behaviour of free end–fixed base columns. A non-prismatic elastic circular column of height L was analysed, taking self-weight into account. Various scenarios were considered: column taper, base fixity, stiffness (E) and density (ρ) variation along the height and cross-section, and crown mass. The effect on the sway natural frequency was assessed in each case. Validation against benchmark problems was conducted satisfactorily. The results indicate that column taper, base fixity, ellipticity and E/ρ ratio are particularly important for this problem. Comparison with published literature results from actual field-based tree swaying tests shows that the model is able to adequately represent the behaviour of conifer trees but is less suited to broadleaved trees with their more complex architecture. Due to the simple presentation, use of the results contained in the paper should improve estimates of natural frequency for trees and allow more subtle characterisation and understanding of the problem.