Effect of Dimensional Variables on the Behavior of Trees for Biomechanical Studies

Abstract

The dimensional variables of trees play an important role in biomechanical studies that seek to estimate the risk of falls, since they influence their biomechanical behavior in relation to the forces to which they are subjected, and therefore, their safety factor. The aim of this research is to evaluate the effect of dimensional variables of trees on their mechanical behavior. A finite-element model was used to perform linear static analysis. The wood from the tree was considered clean, and the architectural model was based on dimensional variables of species commonly used in urban afforestation in São Paulo, Brazil. Different slenderness, tapering, height, and load level were used to analyze the tree mechanical behavior. The numerical-simulation model facilitates the evaluation of the influence of dimensional parameters of trees on deflections and stresses. The behavior of the deflections varies according to height, diameter, and loading level. Since the model considers the geometric variations of the section, the stresses show smooth variations along the trunk. The maximum module values of positive and negative stresses are not equal, and can undergo sudden variations in position along the trunk when local maximum stresses become global maximums.

1. Introduction

In biomechanics studies, it is customary to apply known theories of mechanics developed for conventional structures, including many limitations to represent the plant’s behavior [1,2]. Some of these limitations are the use of equations that are only valid for small deflections and static equations for analyzing the wind load, which is dynamic [1,3]. Wood and most plant materials are viscous and therefore have a nonlinear behavior, making theory and reality not coinciding [2]; even so, the studies are valid, so it is possible to obtain approximations. Another significant limitation is the isotropy of the material, as cited by [1]. This can be minimized with the use of models that allow the application of the orthotropic condition; however, to circumvent this limitation, it is necessary to know the complete elastic parameter of the wood.

The safety factor is another limitation related to applying theories of conventional mechanics in trees. For conventional structures, the safety factor is usually considered the relationship between the load capacity and the acting load [4]. This proportion must be greater than or equal to one for the structure to be reliable. However, many publications indicate that plants are structures that withstand stresses much higher than those expected by the conventional design resistance [5,6,7,8,9] Proposals for measurements of plant safety factors have been developed through studies of the relationship between breakdown and expected service stresses in situations of applied external forces [9]. Thus, the safety factors of trees reported in the literature are generally based on studies applying static loads [8], such as pulling tests. In ref. [10] showed that the safety factors in the trees vary sinusoidally, so that the smallest safety factors occur for peripheral branches and the lower parts of the trunk. As such, the intermediate parts of the tree have significantly greater safety factors. According to the authors, this is a tree strategy so that the branches break first, reducing the crown, and thus reducing the drag forces acting on the trunk, especially in the intermediate parts where the maximum stresses are expected. The literature suggests that trees could withstand stresses three to five times higher than other wooden structures without breaking [5,6,7,8,10,11]. In addition, the wood strength obtained in tests using specimens is not compatible with the stresses supported by the tree trunk, even within the limit of proportionality [12,13], making field trials with whole trees important.

In addition to the limitations related to the theories of mechanics, there are many others related to the material of live trees. Some worth mentioning are the differences in rupture for dry and physiologically inactive wood [14], great capacity to undergo deflections much higher than the elastic-band linear [14], complex geometry, great adaptive ability to respond to external stimuli, as well as changes due to the formation of different types of wood [1].

Although the representation of the wind as static horizontal loads (even though they are dynamic) [15] simplifies the analysis and makes it more accessible, it can introduce some inaccuracies [1,2,3]. In ref. [16] shows that the maximum inclination of the root zone of the structure measured by the static-pulling test was 68% lower than the one obtained with natural wind, demonstrating that the static analysis can underestimate the tree’s response. In ref. [17] also reported that their results show that trees fell with wind speeds 56% lower than the critical speed values obtained through the traction test. The finite-element method (FEM) can be used to study the dynamic behavior of a structure, but this approach is more complex, involving inertial-mass forces, elastic forces, and damping forces, whereas static analysis considers only the elastic forces [1,18]. The dynamic response of the tree is important, but the wind load is also important and quite complex to model, requiring accurate empirical measurements of many parameters of the tree and load conditions to obtain reliable results [1,18]. Despite knowing the limitations and errors, static analyses to simulate the wind are still performed due to the inherent difficulties of dynamics analyses. In addition, it is important to note that the oscillation induced by the wind load alters the effort produced in the trunk, sometimes reducing it, and this favorable condition is not considered when using horizontal loads to simulate the wind [1,3,19]. Furthermore, dynamic loads are short-term loads, while static loads are relatively long-term [1]. These two aspects (canopy movement and short-term load) can play a favorable role in reducing the load transmitted to the trunk and improving the response of the wood, respectively, indicating that static analysis will not always underestimate the results. New results using the digital-image-correlation system (DIC) may be a future way of mapping strains in standing trees, even on complex tree surfaces, allowing for the simplification of dynamic analysis in biomechanical studies on trees [18].

Static-load tests using a cable attached to the tree to apply a controlled tension (pulling test) have been widely used to induce nondestructive bending stress in the linear elastic range and thus extrapolate the critical-bending moment. In ref. [12] shows that the material properties determined in small specimens are less relevant for evaluating tree strength than values derived from the test of whole tree stems, indicating that despite the limitations, static-field methods are important tools to be used in biomechanical studies. However, this test does not consider that plastic deformations frequently occur in live trees, increasing their carrying capacity [1]. This test is also used to estimate the resistance to uprooting by models, relating the bending moments necessary to induce small changes in the slopes of the soil–root plate. These models are used to extrapolate failure loads but have been criticized because the angle of rotation of the root plate at maximum strength varies with the age of the tree, root architecture, and soil structure [1]. In ref. [17] concluded that it is necessary to solve some problems related to the pulling-test equipment and that further research is needed to know the drag coefficient, the pedohydraulic condition of the soil, and the relationship between static and dynamic conditions. Some authors, cited in the review by [1], show that static-load tests will advance by using optical digital 3D-image correlation (DIC) to measure strain and stress in branches, trunks, and roots, allowing further studies in bending, torsion, and fracture in woody plants and further in-depth understanding of strain during static-load tests, allowing consideration of the impact of plastic deformations. Refs. [19,20,21] are examples of the incorporation of DIC in static-pulling tests. Both researchers show that the DIC method can be used successfully to obtain high-quality strain data in tensile tests.

In addition to research aimed at obtaining improvements in static-field tests, researchers have been studying ways to use numerical models to identify parameters that mainly impact the response of trees [22] and estimate the distribution of deformations and stresses along the trunks of trees exposed to gravitational and horizontal loads [23,24]. These models can help evaluate the effects of the architectural characteristics of the roots and soil properties on the distribution of loads around the root system [25] as tools in analyzing tree ruptures [26]. They can also help in the study of the influence of a series of parameters (mechanical properties and geometry of the root system and the stem) on stress, strain, and displacement of the tree model [24]. Numerical models are also being studied to represent the branches [27] better and to simulate the dynamics of the tree [28], and to uproot with wood-breaking behavior [29].

The dimensional variables of the trees (base diameter, diameter at breast height (DBH), height, length of branches, etc.) can be measured directly. Thus, in addition to being an important tool for the forestry industry, the simplicity of its measurement can be a valuable tool to assess the risk of falling urban trees in biomechanical analyses, including through numerical models. Geometry and stem dimensions significantly affect the dynamic response [30] and stresses, strain, and displacements [24] of trees subjected to wind loads than their wood properties. These authors also concluded that nonlinear behaviors are a function of the architecture of the branches. Considering that wind is the most critical load in analyzing fall risk [1,31], this result indicates the importance of tree dendrometry. In predictive models for falls, it is important to know not only the mechanical characteristics of the stem and branches, but also the dimensions and shape of the trees [1,22,31,32,33].

The objective of this research was to evaluate, through numerical simulation, the effect of dimensional tree variables on deflection and normal and shear stresses and to discuss the consequences on the safety factor of trees.

2. Materials and Methods

2.1. Finite-Element Model

The dendrometric characteristics of the trees, such as height, slenderness, and taper, play an essential and complex role in the risk of falling trees. Therefore, to make more viable the analysis of the influence of these dendrometric parameters in the deflection and in normal and shear stresses, this study proposed to create a numerical modeling using the finite-element method (FEM). The model was developed using the Ansys^® Academic Research Mechanical, Release 18.2 software. The system was composed of the trunk, root, and soil bin (Figure 1). The FEM model was created using a mesh with solid and quadratic 3D elements (Figure 1), varying in number according to the size of the tree and the type of element (trunk, root, or soil bin), with a maximum size of 200 mm for the tree and 300 mm for the soil bin. The displacement and stress fields are computed considering the linear static analysis.

The root was simulated as a 0.15 m deep disc, as proposed by [34]. The root diameter varied according to the crown. The crown diameter varied according to the tree’s height and was fixed at 5.5 m in the scenario where the height is fixed. The soil bin proposed by [26] also had a fixed dimension (16 m × 16 m), but in our model, the cubic-shape soil bin was simulated with a dimension proportional to the root (2.62 times). In addition, as proposed by [26], frictions between soil and root were ignored, and for boundary conditions, only the base was considered fixed. The trunk (without crown and branches) was created as conical (larger at the bottom) instead of cylindrical, as proposed by [10] (Figure 1).

To propose the tree model (architecture) and the range of the dimensions of the trunk and branches, we used some results from [35] as a reference. The study [35] was carried out with six species of trees (Handroanthus pentaphylla (Ipê Rosa), Cenostigma pluviosum (Sibipiruna), Tipuana tipu (Tipuana), Schinus molle (Aroeira Salsa), Caesalpinia ferrea (Pau Ferro), and Schinus terenbinthifolius (Aroeira pimenteira), which are all commonly planted in the state of São Paulo, Brazil.

As the proposed model uses only the conical trunk (without canopy), it was necessary to consider, in a simplified way, the effect of the branches and leaves. Thus, the weight of branches and leaves was considered as a concentrated load at the end of the branch, allowing the calculation of a moment, which was applied to the trunk. To obtain reference data for calculating the weight of branches, we also used the results (diameter, length, tapering, and density of tree branches) obtained in previous research [3,35]. We added 4.0% of the weight of the branch to consider the weight of the leaves, as proposed by [36]. The mean value of calculated moments (969 N·m) was applied at three levels, representing the most common tree architecture obtained in the previous study [35]. The height levels varied according to the total height of the tree (

Table 1. Scenarios with total tree height (H), slenderness (λ = H/DBH), taper (T = (Db-Dt)/H), diameter at breast height (DBH), diameter at the base (Db) and at the top (Dt) of the trunk, and position (Y) of the moments (YM₁, YM₂, YM₃) and of the horizontal forces (YFH) applied in the numerical simulations.

Scenario	H m	λ	T	DBH m	Db m	Dt m	YM1 m	YM2 m	YM3 m	YFH m
1	5	22	0.030	0.23	0.27	0.116	2.00	3.00	4.00	3.50
	10	22	0.030	0.45	0.49	0.194	2.00	4.67	7.33	6.0
	15	22	0.030	0.68	0.72	0.271	2.00	6.33	10.67	8.50
	20	22	0.030	0.91	0.95	0.348	2.00	8.00	14.00	11.00
	25	22	0.030	1.14	1.18	0.425	2.00	9.67	17.33	13.5
2	7.5	8	0.030	0.94	0.98	0.752	2.00	3.83	5.67	4.75
	7.5	16	0.030	0.47	0.51	0.283	2.00	3.83	5.67	4.75
	7.5	24	0.030	0.31	0.35	0.127	2.00	3.83	5.67	4.75
	7.5	32	0.030	0.23	0.27	0.048	2.00	3.83	5.67	4.75
	7.5	40	0.030	0.19	0.23	0.002	2.00	3.83	5.67	4.75
3	7.5	22	0.010	0.34	0.35	0.279	2.00	3.83	5.67	4.75
	7.5	22	0.020	0.34	0.37	0.217	2.00	3.83	5.67	4.75
	7.5	22	0.030	0.34	0.38	0.155	2.00	3.83	5.67	4.75
	7.5	22	0.050	0.34	0.41	0.031	2.00	3.83	5.67	4.75
Validation	9.4	39	na	0.24	0.32	na	3.00	5.13	7.23	6.20

Y measured from base to top; moments (M₁, M₂, and M₃) to simulate the branches in three levels. na = not available. Scenario 1: slenderness 22; taper 0.030; and total tree height ranging between 5 m, 10 m, 15 m, 20 m, and 25 m. Scenario 2: tree height 7.5 m; taper 0.030; and slenderness ranging between 8, 16, 24, 32, and 40. Scenario 3: tree height 7.5 m; slenderness 22; and taper ranging between 0.010, 0.020, 0.030, and 0.050.

). The first level was always at the first fork (the point where the canopy begins), which was set at 2.0 m (average value obtained in [35]). From that point, the canopy height (H-h) was divided into three parts to place the remaining levels of branches (Table 1). For each level, four branches were adopted around the trunk: two in the X direction, with moments around Z; and two in the Z direction, with moments around X (Figure 1). Moments on opposite faces were considered with opposite signs (X1 positive and X2 negative and Z1 positive and Z2 negative).

A reference density (994 kg·m⁻³) for trunk wood, obtained in [35], was introduced in the FEM model to calculate the trunk self-weight.

Three scenarios were proposed to evaluate the effect of dimensional variables on the behavior of trees subjected to moments, gravitational and horizontal loads. In the first scenario, total trunk height (H) varied while slenderness (l) and taper (T) values were fixed (Table 1). In the second scenario, l varied while H and T values were set, and, finally, in the third scenario, T varied while H and l values were fixed (Table 1).

2.2. Horizontal Load

For simplification purposes, the dynamic wind loads given by the wind pressure (q) were transformed into a horizontal load (P_HV) (Equation (1)), applied to the center of gravity of the crown. This simplification is common in the literature [3].

$$P_{HV}=q{A}_c$$

(1)

where A_c is the area of the front surface of the crown.

Wind pressure (q) depends on air density (ρ) and wind speed (V) according to Equation (2).

q = ½ ρ V ²

(2)

Assuming that the air density is approximately 1225 kg·m⁻³ (at sea level and 15 °C) and considering the acceleration of gravity (m·s⁻²), the value of the wind pressure (q in N·m⁻²) is obtained (Equation (3)).

$$q=0.613 V^2$$

(3)

Therefore, horizontal loads depend on wind pressure and canopy area (Equation (1)). The wind pressure depends on the wind speed, while the crown area depends on the size and architecture of the crown. As [35] shows that the most common tree architecture was the monopodial shaped like a globose (round) crown, the crown was simplified to a circumference (Figure 1) with a diameter equal to the height of the tree minus the height of the first fork assumed in this work as the average value obtained by [35], which was two meters.

Assuming a range of wind speeds, we also obtain a range of wind pressures (Equation (3)) and, therefore, horizontal loads (Equation (1)). We started the procedure testing a range of wind speeds covering speeds below, within, and above basic wind speeds (peak gust speed averaged over a short time interval of about 3 s duration and measured at 10 m above ground) used in the Brazilian standard, [37] which are 30 m·s⁻¹ to 50 m·s⁻¹. We calculated the wind pressure (Equation (2)) using this wind-speed range (from 30 m·s⁻¹ and 50 m·s⁻¹), and the horizontal load (Equation (1)) was calculated, with frontal canopy area (Ac in

Table 2. Canopy parameters (diameter = D_c and area = A_c) of trees with different heights (H), adopted in the three scenarios, with respective values of wind pressure (q = P_HV/A_c) calculated for the different levels of horizontal load (P_{HV =} 5 kN, 10 kN and 15 kN) and associated wind speeds (V = 0.613 V²).

				5 kN		10 kN		15 kN
Scenario	H m	Dc m	Ac m2	q N/m2	V m/s	q N/m2	V m/s	q N/m2	V m/s
1	5	3	7.1	707	34.0	1415	48.0	2122	58.8
	10	8	50.3	99	12.7	199	18.0	298	22.1
	15	13	132.7	38	7.8	75	11.1	113	13.6
	20	18	254.5	20	5.7	39	8.0	59	9.8
	25	23	415.5	12	4.4	24	6.3	36	7.7
2 and 3	7.5	5.5	23.8	210	18.5	421	26.2	631	32.1
Validation	9.4	6.4	32.2	155	15.9	311	22.5	466	27.6

Scenario 1: slenderness 22; taper 0.030; and total tree height ranging between 5 m, 10 m, 15 m, 20 m, and 25 m. Scenario 2: tree height 7.5 m; taper 0.030; and slenderness ranging between 8, 16, 24, 32, and 40. Scenario 3: tree height 7.5 m; slenderness 22; and taper ranging between 0.010, 0.020, 0.030, and 0.050. Initial crown height fixed at 2.0 m (crown heights = H − 2.0).

). As we need to apply the same horizontal load to all scenarios (Table 1), the second analysis was to use the horizontal-load range (obtained for the different areas of canopy) to calculate the respective wind pressure (Equation (1)) and so the respective wind speeds capable of producing this level of load (Equation (3)). Since below 5 kN and above 15 kN the wind speed was, respectively, far below (for the largest canopy) and far over (for scenario 1) the proposed wind-speed range, we choose the three load levels (5 kN, 10 kN, and 15 kN) to be applied to all scenarios. According to the Beaufort scale, at speeds from 5.5 to 7.9 m·s⁻¹ the wind moves small branches; from 8 to 10.7 m·s⁻¹ it moves large branches and small trees; from 13.9 to 17.1 m·s⁻¹ it moves large trees; from 17.2 to 20.7 m·s⁻¹ it can break branches; from 20.8 to 24.4 m·s⁻¹ it can damage trees; and from 24.5 m·s⁻¹ it can uproot trees. The chosen load levels allow us to simulate wind speeds from 4.4 m·s⁻¹ to 58.8 m·s⁻¹ (Table 2), but with the exception of trees with H = 5 m where the speeds are greater than 30 m·s⁻¹, all other conditions are simulating at speeds below 22 m·s⁻¹ (Table 2), i.e., in the range where large tree damage is not expected by the Beaufort scale. Furthermore, for shorter trees, greater stability is expected due to their lesser slenderness, and, therefore, these short trees will likely be better able to withstand greater wind forces. The same wind-speed level was used by [10], for example, who simulated wind speeds of 10, 20, and 50 m·s⁻¹.

These horizontal loads were applied to the tree model in a position coincident with the gravitational center of the crown (Table 1).

2.3. Material Properties

Soil was considered isotropic, and the values obtained by [34] were adopted (Density = 1730 kg·m⁻³, Modulus of Elasticity = 13.5 MPa and Poisson’s ratio = 0.4).

The trunk and root wood were considered orthotropic, and the values adopted in the simulations were the average physical and mechanical properties in the saturated condition, studied in the research group via ultrasound test in polyhedral specimens. This methodology [38,39] allows us to obtain, with ultrasound, the stiffness matrix and the inversion, the compliance matrix; therefore, the nine elastic constants of the wood (three longitudinal moduli, three shear moduli, and six Poisson’s ratio). Thus, for the numerical simulation, considering the material’s orthotropy, three elasticity moduli (longitudinal (E_L), radial (E_R), and tangential (E_T)); three shear moduli in the planes (G_RT, G_LT, and G_LR); and three Poisson’s ratio in the planes (ν_LR, ν_LT, and ν_RT) are inserted. In general, elastic constants of wood determined by elasticity theory consider the longitudinal direction as X, the radial direction as Y, and the tangential direction as Z. We use this same terminology for the local coordinate system. In the complete characterization of this material, six Poisson’s ratio were obtained (ν_RL, ν_TL, ν_LR, ν_TR, ν_LT, and ν_RT) (

Table 3. Average elastic parameters considered in the numerical simulation.

Elastic Parameters	Directions or Planes
Longitudinal modulus of elasticity	Longitudinal	$E_L$ [MPa]	10387
	Radial	$E_R$ [MPa]	1908
	Tangential	$E_T$ [MPa]	1290
Shear module	Radial/Tangential	$G_{RT}$ [MPa]	515
	Longitudinal/Tangential	$G_{LT}$ [MPa]	1173
	Longitudinal/Radial	$G_{LR}$ [MPa]	1573
Poisson’s ratio	Radial/Longitudinal	${\nu }_{RL}$	0.09
	Tangential/Longitudinal	${\nu }_{TL}$	0.08
	Longitudinal/Radial	${\nu }_{LR}$	0.49
	Tangential/Radial	${\nu }_{TR}$	0.46
	Longitudinal/Tangential	${\nu }_{LT}$	0.64
	Radial/Tangential	${\nu }_{RT}$	0.78

). However, the FEM only allows the incorporation of three values, as it uses reciprocal relationships internally (Equation (4)):

$$\frac{{\nu }_{TL}}{E_T}=\frac{{\nu }_{LT}}{E_L};\frac{{\nu }_{RL}}{E_R}=\frac{{\nu }_{LR}}{E_L};\frac{{\nu }_{RT}}{E_R}=\frac{{\nu }_{TR}}{E_T}$$

(4)

Given that the ν_RL and ν_TL Poisson’s ratio are the lowest (Table 3) and in practice are obtained with less experimental precision, the calculations were performed, providing the program with the values of ν_LR (0.49) and ν_LT (0.64). In the case of the RT plane, we followed the same logic, using the highest value (ν_RT = 0.78).

Although the model allows the use of different properties for root and trunk, the same values (Table 3) were used for both because we had no properties available for root.

2.4. FEM Simulation

With the simplified topology adopted (Figure 1) and the wood parameters (Table 3), for each scenario (Table 1), the deflection and normal (parallel to fiber) and shear (LR plan) stresses (positives and negatives), as well as the position of these values along the trunk, were obtained, using the FEM (Ansys R18.2). The FEM simulation model used in this paper considers the deflection of all systems (root-plate and stem) (Figure 2).

2.5. Safety Factors

For all scenarios, safety factors were analyzed using the traditional expression for structures (relation between load-bearing capacity and active load).

The active loads were obtained in FEM simulations (maximum normal and shear stresses), and the load-bearing capacity using average strength data in saturated condition obtained by associating the strength classes of the European Standard [40] to E_L, obtained as mentioned previously (Table 3). For this association, we corrected the E_L obtained in saturated conditions (Table 3) to the reference moisture content (12%) using the European Standard [41] (1% for each 1% variation around the reference value up to 20%) and used the value of corresponding strength (f_m,k = 39 MPa). The relations proposed in the European standard [41] were also used to infer the characteristic shear strengths (f_v,k = 0,2 (f_m,k)^0.80 = 3.8 MPa). The [42] was used to determine the design values (Equation (5)), considering the modification coefficient (k_mod) for solid wood, class of service 3 (totally exposed directly to rainwater and climatic conditions), and short-term loading class (permanent + wind), resulting in k_mod = 0.70. For the partial safety coefficient, γ_m = 1.3 (solid wood) was used.

$$f_{w,d}=k_{mod}\frac{f_{w,k}}{{\gamma }_m}$$

(5)

where: f_w,d = design resistance; k_mod = modification coefficient; and γ_m = partial safety coefficient.

Thus, considering engineering parameters applied to conventional wooden structures, the design resistance, using strength in the saturated condition, was inferred in bending (f_m,d = 21 MPa) and shear (f_v,d = 2 MPa). With the design resistances considered as load-bearing capacity in bending and in shear as well as the acting forces (normal and shear stress) obtained from the simulations, the safety factors for all scenarios were calculated and discussed based on the specific behavior of the trees presented in the literature.

3. Results and Discussion

3.1. Deflection

3.1.1. Comparison with Literature

Comparing deflection values using literature is complex because different authors use different measurement methodologies. Some of them consider only the inclination of the root plate [15,17], and others consider the inclination of the root plate and the deflection of the stem together [26,33,43], as in this study (Figure 2). Furthermore, the deflection is influenced by geometric tree parameters such as height, diameter, and taper [24], and also by mechanical properties of the trunk and root wood, as modulus of elasticity [15]. Therefore, the comparison with literature is focused on the analysis of the general order of magnitude of the values obtained and the coherence of behavior, considering geometric parameters of the trees and properties of the wood adopted by different authors.

In ref. [43] uses two trees with the same height (H = 20 m) and diameters at breast height (DBH = 24.5 cm and 24.8 cm), both same species with density 300 kg·m⁻³ and modulus of elasticity in bending 3584 MPa. The horizontal pulling force, with a magnitude much lower (maximum 2.25 kN and 1.65 kN) than those used in this research (5 kN to 15 kN), was applied at 12.23 m and 12.20 m (about 60% H). The closest condition to compare our data is in scenario 1, where the tree with H = 20 m (Table 1) was subjected to a horizontal load of 5 kN at 55% H (

Table 4. Values and positions (Y = %H from base to top) of maximum deflection (δ_M) and deflection in the horizontal load (δ_L) position (mm), maximum negative (σ₋) and positive (σ₊) normal stresses (MPa), and maximum negative (σ₋) and positive (σ₊) shear stresses (MPa) for the height (H), slenderness (λ), and taper (T × 102) values adopted in each of the simulated scenarios. Horizontal load 5 kN.

H	λ	T	δM	Y	δL	Y	σ−	Y	σ+	Y	τ−	Y	τ+	Y
5	22	3	13	100	10	70.0	2.70	40	2.72	40	0.54	40	0.47	40
10	22	3	58	100	57	60.0	0.85	19	0.70	19	0.50	46	0.33	19
15	22	3	152	100	152	56.7	1.68	0	0.74	70	0.78	0	1.02	0
20	22	3	293	100	293	55.0	2.51	0	2.00	0	1.99	0	1.52	0
25	22	3	440	100	440	54.0	3.45	0	2.51	0	2.19	0	2.14	0
7.5	8	3	29	100	29	63.3	0.74	0	0.17	75	0.29	0	0.27	0
7.5	16	3	27	100	27	63.3	0.66	63	0.37	26	0.29	26	0.27	75
7.5	24	3	29	100	27	63.3	2.20	26	1.79	26	0.56	26	0.58	51
7.5	32	3	34	100	30	63.3	4.40	26	4.70	26	1.01	26	0.88	26
7.5	40	3	60	100	42	63.3	9.32	26	10.38	26	2.50	26	2.94	26
7.5	22	1	28	100	27	63.3	0.98	26	0.60	26	0.32	26	0.30	26
7.5	22	2	28	100	27	63.3	1.03	26	1.30	26	0.29	26	0.25	26
7.5	22	3	28	100	27	63.3	1.13	26	1.12	26	0.38	26	0.32	26
7.5	22	4	28	100	27	63.3	1.51	63	2.18	26	0.51	26	0.49	51

). For a horizontal load of 1.65 kN, the average displacement at the point of load application obtained by [43] was 1.5 m. If we apply a composition of differences of DBH (24.7 cm/91 cm = 0.27) and modulus of elasticity (3584 MPa/10387 MPa = 0.34) calculated by multiplying these differences (0.09), the expected displacement (0.14 m) will be of the same order of magnitude as those obtained in this paper (0.293 m in Table 4).

In ref. [34] simulated a tree with H = 8 m, DBH = 0.40 m, and horizontal force applied at 5 m (62.5% H). The simulated horizontal forces ranged from 5 kN to 75 kN. There is no information about the modulus of elasticity of the stem and root used in the model. The closest condition in this research is in scenario 2 for slender = 16, where H = 7.5 m and DBH = 0.47 m (Table 1) with horizontal load applied at 63.3% H. For this condition, the deflection in the horizontal-load position was 0.027 m (horizontal load = 5 kN (Table 4)), 0.028 m (horizontal load = 10 kN (

Table 5. Values and positions (Y = %H from base to top) of maximum deflection (δ_M) and deflection in the horizontal load (δ_L) position (mm), maximum negative (σ₋) and positive (σ₊) normal stresses (MPa), and maximum negative (σ₋) and positive (σ₊) shear stresses (MPa) for the height (H), slenderness (λ), and taper (T × 10²) values adopted in each of the simulated scenarios. Horizontal load 10 kN.

H	λ	T	δM	Y	δL	Y	σ−	Y	σ+	Y	τ−	Y	τ+	Y
5	22	3	22	100	16	70.0	5.35	40	5.48	40	1.08	40	0.93	40
10	22	3	59	100	58	60.0	1.58	19	1.53	19	0.50	46	0.33	19
15	22	3	153	100	153	56.7	1.69	0	0.74	12	0.78	0	1.02	0
20	22	3	294	100	294	55.0	2.51	0	2.00	0	1.99	0	1.52	0
25	22	3	440	100	440	54.0	3.46	0	2.51	0	2.19	0	2.14	0
7.5	8	3	29	100	29	63.3	0.89	63	0.20	63	0.29	0	0.27	0
7.5	16	3	29	100	28	63.3	1.30	63	0.80	26	0.29	26	0.27	75
7.5	24	3	34	100	30	63.3	4.49	26	3.51	26	1.08	26	1.12	51
7.5	32	3	51	100	38	63.3	8.74	26	9.45	26	2.02	26	1.77	26
7.5	40	3	111	100	69	63.3	18.72	26	20.74	26	5.02	26	5.80	26
7.5	22	1	31	100	29	63.3	1.81	26	1.37	26	0.44	0	0.47	0
7.5	22	2	31	100	29	63.3	2.03	26	2.54	26	0.50	26	0.51	26
7.5	22	3	31	100	29	63.3	2.22	26	2.30	26	0.77	26	0.67	26
7.5	22	4	33	100	29	63.3	3.03	63	4.32	26	1.03	26	0.95	26

)), and 0.029 m (horizontal load = 15 kN (

Table 6. Values and positions (Y = %H from base to top) of maximum deflection (δ_M) and deflection in the horizontal load (δ_L) position (mm), maximum negative (σ₋) and positive (σ₊) normal stresses (MPa), and maximum negative (σ₋) and positive (σ₊) shear stresses (MPa) for the height (H), slenderness (λ), and taper (T × 10²) values adopted in each of the simulated scenarios. Horizontal load 15 kN.

H	λ	T	δM	Y	δL	Y	σ−	Y	σ+	Y	τ−	Y	τ+	Y
5	22	3	31	100	21	70.0	8.00	40	8.24	40	1.62	40	1.38	40
10	22	3	62	100	59	60.0	2.31	46	2.35	19	0.50	19	0.44	19
15	22	3	155	100	153	56.7	2.03	0	0.77	0	0.78	0	1.02	0
20	22	3	295	100	294	55.0	2.64	0	2.00	0	1.99	0	1.52	0
25	22	3	441	100	441	54.0	3.53	0	2.51	0	2.19	0	2.14	0
7.5	8	3	29	100	29	63.3	1.15	0	0.31	0	0.28	0	0.29	0
7.5	16	3	30	100	29	63.3	1.95	26	1.23	26	0.41	26	0.30	26
7.5	24	3	41	100	33	63.3	6.78	26	5.23	51	1.59	26	1.65	51
7.5	32	3	71	100	48	63.3	13.08	26	14.20	26	3.03	26	2.66	26
7.5	40	3	164	100	99	63.3	28.12	26	31.09	26	7.55	26	8.67	26
7.5	22	1	36	100	31	63.3	2.68	26	2.14	26	0.66	26	0.71	26
7.5	22	2	36	100	31	63.3	3.03	26	3.78	26	0.76	26	0.78	26
7.5	22	3	36	100	31	63.3	3.32	26	3.48	26	1.16	51	1.02	26
7.5	22	4	39	100	32	63.3	4.56	26	6.46	51	1.55	26	1.43	51

)). For these load levels, [34] obtained 0.02 m, 0.03 m, and 0.04 m, with the same order of magnitude as our results.

In ref. [33] modeled, using the structural theory of a cantilever beam, the behavior of eight trees with H ranging from 9.6 m to 13.2 m subjected to large deflections using a pulling test. The modulus of elasticity of the trees’ wood was 2000 to 6400 MPa, and the diameter (D) was measured at 2 m from the base of 10 cm to 19.2 cm. The horizontal load was applied at about 70% H. The deflection of the trees ranged from 2 m to 6 m and showed positive dependence on the modulus of elasticity and diameter. The mean deflection of two trees (H = 9.6 m, DBH = 15.5 cm, and modulus of elasticity = 4000 MPa) tested by [33] is 2.5 m. We calculated a factor to compare the tree in our scenario 2 (H = 10 m, DBH = 45 cm and modulus of elasticity = 10,387 MPa and horizontal load applied at 60% H) by multiplying the diameter ratio (0.34) and modulus of elasticity ratio (0.39) and the expected deflection is 0.33 m for these trees. This value is on the same order of magnitude as that obtained by the model proposed in this paper (0.57 m).

The results of the comparisons show that although at first glance, the deflections obtained with the proposed model are very small, the values obtained by adjusting the diameter and modulus of elasticity are within the expected order of magnitude.

3.1.2. Behavior of Trunk Deflection with Tree Height

The deflection calculated by the model (Table 4, Table 5 and Table 6) is complex because it involves root-plate and trunk deflection under the interactions of self-load, horizontal load (simulating wind), and moments (simulating the branches’ effect). The trunk self-load only causes compression stress if the trunk is in a vertical position, but with the deflection of the trunk, this self-load will also have a horizontal component that will add bending. Moments representing branches also increase bending above the first bifurcation, and part of these moments are acting above the horizontal load position. Horizontal load causes bending of the trunk below this point, which is proposed as the center of gravity of the canopy. In a simplified way, until the break an increasing linear deflection with the height from the base to the top of the trunk is expected, but as the deflection is due to all load components in addition to the variation in dendrometric characteristics (conical trunk), this was not the behavior obtained for all scenarios.

In scenario 1, only the smaller tree (H = 5 m) showed increasing deflection with the tree’s height to the top (from 0%H to 100%H) (Figure 3). Independent of the load level, deflection values in scenario 1 for trees with a height of 10 m to 25 m show deflections at the top (dM) very similar (≥95%) or even equal to those at the point of horizontal-load application (dL), maintaining the line practically constant for the section above the point of application of the horizontal load (Figure 3). For the smaller tree (H = 5 m), the differences between the deflection at the top and at the horizontal-load point increase as the load level increases (horizontal-load position 76.9% of the top position for 5 kN horizontal load, 72.7% for 10 kN horizontal load, and 67.7% for 15 kN horizontal load) (Table 4, Table 5 and Table 6). This behavior is probably related to the association of the magnitude of the top diameter and the concentration of moments arising from the branch effect. In taller trees, the diameters are greater than about 0.2 m (Table 1), and as the position of the first bifurcation is fixed at 2 m, the height of the canopy increases with increasing tree height, and thus increasing the distance from the points of the trees where moments are applied, reducing its concentration in the canopy.

In scenario 2, the differences between the deflection in the horizontal-load position and the top position increase with increasing slender, and are influenced by the load level (Figure 4). For the first load level (5 kN), the only tree with l = 40 has the deflection growing from the base to the top. As the load level increased to 10 kN and 15 kN, this behavior started to be observed in the tree with l = 32 and l = 24, respectively (Figure 4). The slender is related to the height of the tree and the diameter at breast height, which is associated with the diameter of the top (conical trunk). In this scenario, the tree’s height is fixed (7.5 m), so the distance of the moments from branches is also fixed. The most significant differences of deflection at the top and the horizontal load point came from top diameters less than 0.127 m (above slender = 24, Table 1).

In scenario 3, as in scenario 2, the height of the tree is fixed (7.5 m); all conditions have the same moment concentrations, but the top diameter decreases as the taper increases, and as in the other two scenarios, more significant differences are obtained for the smaller top diameter. Up to 10 kN load level, the behavior is the same for all taper levels, keeping the deflection line practically straight for the section above the horizontal-load-application point (Figure 5). The behavior in scenario 3 started to change to load = 15 kN.

In ref. [33] presents graphics of stem deflection at different heights for eight trees showing variation in the slope of the curve from the point of application of the horizontal load, which in some trees is also straighter.

In the European Standard [44] focused on timber design, the maximum instantaneous deflection allowed for the cantilever condition, used to represent the condition of trees, is the length of the piece (L) by 150 or 250. For heights above 15 m in scenario 1 for all loading levels (Table 4, Table 5 and Table 6), the deflection was not compatible with this limit. However, it is important to remember that these limits, set by timber design standards, are only a reference, as the plants have an adaptive capacity that makes them capable of deforming well beyond the usual elastic limits of structures [7].

3.2. Normal Stresses

For maximum normal stresses (positive and negative) (Figure 6), the influence of slenderness is the most significant, followed by tapering. Regarding tapering, ref. [24] found stress influence for trees with heights bigger than four meters, corresponding to all of our scenarios. Height is the least influencing factor. Height variation had a small impact mainly on trees with up to 10 m (Figure 6). Statistically, height had less influence on normal negative stresses (nonsignificant regression model at all load levels), and its effect on normal positive stresses depended on the load level (

Table 7. Regression models considering maximum (negative and positive) normal stresses (σ) as a function of height (H), slenderness (λ), and taper (T × 10²) for each level of horizontal load adopted.

Variables	Regression Model	Load Level kN	p-Value	R2 (%)
Negative Normal Stresses (σ−)
σ− × H	-	5	0.2498 *	-
	-	10	0.2121 *	-
	-	15	0.0634 *	-
σ− × λ	σ−= 0.54 − 0.006 λ2	5	0.0047	95.1
	σ−= 0.54 − 0.006 λ2	10	0.0034	96.1
	σ−= 0.54 − 0.006 λ2	15	0.0032	96.2
σ− × T	σ−= 0.54 − 0.006 T2	5	0.0014	99.7
	σ−= 0.54 − 0.006 T2	10	0.0017	99.7
	σ−= 0.54 − 0.006 T2	15	0.0020	99.6
Positive Normal Stresses (σ+)
σ+ × H	-	5	0.5016 *	-
	σ+ = sqrt (−6.4 + 65/H)	10	0.0474	77.9
	σ+ = sqrt (−20.7 + 411/H)	15	0.0209	86.9
σ+ × λ	σ+ = −1.25 + 0.07 λ2	5	0.0042	95.4
	σ+ = −2.56 + 0.013 λ2	10	0.0039	95.6
	σ+ = −3.83 + 0.020 λ2	15	0.0039	95.6
σ+ × T	σ+ = sqrt (0.35 + 0.0017 T2)	5	0.0409	92.0
	σ+ = sqrt (1.61 + 0.0066 T2)	10	0.0329	93.5
	σ+ = sqrt (3.81 + 0.015 T2)	15	0.0303	94.0

* not significant with 95% confidence level.

).

In scenario 1, for the trees with heights ranging from 15 m to 25 m, the maximum normal stresses (positive and negative) always occurred at the base (Table 4, Table 5 and Table 6). However, for the shortest trees (5 m and 10 m), the maximum normal stresses rose to higher trunk positions, reaching more central regions (Table 4, Table 5 and Table 6). Since our model considered the tapered trunk, the geometric parameters that go into calculating the stresses varied, even though the moments generated by the horizontal load are greater at the base. In ref. [45], using mathematical models, showed that the variation in dimension and mass along the trunk influences the stress distribution. In ref. [10] simulated a 13 m tree composed of three cylindrical elements with successive diameter reduction, without branches, subjected to horizontal-static loads to simulate wind speed of 10, 20, and 50 m·s⁻¹. The authors’ results showed that although the bending moments increased from the top to the bottom of the tree and the maximum moment was located near the base, the normal stresses varied along the stem following a third-degree polynomial. The results obtained by [33] show that the maximum normal stresses occur between 19% and 63% of the total height of the tree. The same study also shows that shortest trees had the highest maximum stress position (63% H), and the stress distribution along the trunk generally followed a fourth-degree polynomial. In ref. [46] found that uniformly tapered stems (i.e., those of Cryptomeria japonica trees) subjected to wind and snow loads had breaks located between 50% and 70% of the tree’s height, while In ref. [43] showed trunk break at 13.5% H. The distribution of stresses along the stem is influenced not only by the horizontal load considering the wind, but also by the moments applied considering the branches and by the horizontal components of the self-weight, which appear due to the deflection of the trunk. Furthermore, stresses are affected by diameter variation, which impacts the moment of inertia of the transversal section of the trunk. These loads are composed of the horizontal load, changing the behavior of the stresses along the trunk [47].

In scenario 2, to have slenderness variation keeping the height (7.5 m) and the taper (0.030) fixed, it was necessary to vary the DBH from 0.19 m to 0.94 m, reducing the maximum stresses in smaller l conditions that correspond to the highest DBH. A compilation of the results of pulling tests with more than 4500 trees carried out in ten different countries shows that, as expected, taller trees are subjected to higher wind loads (Wessolly, 2004 cited by [11]), with an increase in the resulting moment in the base. However, the stem diameter strongly influences the resulting stresses (Wessolly, 2004 cited by [11]). The maximum normal stresses (positive and negative) in this scenario range from 0% H to 63% H (Table 4, Table 5 and Table 6). Studies by [33] and [10] attribute the variation of stress along the stem to the slenderness, with a deflection from the base (where the maximum moments occur due to the action of wind forces) for the upper regions of the stem. These authors also reject the hypothesis of constant stress distribution along the stem, and conclude that such distribution is strongly dependent on the relationship between height and diameter.

In scenario 3, with defined height (7.5 m) and slenderness [22] and varied taper, normal stresses (positive and negative) grow as the tapering increases, and the position of these maximums varies from 23% H to 63% H (Table 4, Table 5 and Table 6).

3.3. Shear Stresses

The behavior of the position for maximum shear stresses along the trunk, simulated in Scenario 1 (trees of different heights), was similar to that of normal stress (Figure 7); that is, maximum stresses (positive and negative) at the base for trees ranging from 15 m to 25 m high, and in higher positions for lower trees (H = 5 m and 10 m) (Table 4, Table 5 and Table 6). In scenario 2, the maximum shear stress occurs at the base only for the minimum λ (correspondent to the maximum DBH) for all levels of load (Table 4, Table 5 and Table 6). In scenario 3, the maximum shear stresses (positive and negative) increase with taper growth, while the position of the maximum varies from 26% H to 51% H. Statistically, height had an influence on the negative shear stresses only for the lowest horizontal load (5 kN) and on the shear positive stress for the two highest horizontal-load levels (5 and 10 kN) (

Table 8. Regression models considering maximum (negative and positive) shear stresses (τ) as a function of height (H), slenderness (λ), and taper (T × 10²) for each level of horizontal load adopted.

Variables	Regression Model	Load Level kN	p-Value	R2 (%)
Negative Normal Stresses (σ−)
τ− × H	τ− = −0.31 − 0.0032 H2	5	0.0142	89.8
	-	10	0.0631 *	-
	-	15	0.2136 *	-
τ− × λ	τ− = 0.07 − 0.0014 λ2	5	0.0135	90.2
	τ− = 0.42 − 0.0031 λ2	10	0.0084	92.8
	τ− = 0.72 − 0.0046 λ2	15	0.0072	93.5
τ− × T	τ− = −0.29 − 0.00009 T2	5	0.0321	93.7
	τ− = −0.44 − 0.00024 T2	10	0.0332	93.5
	τ− = −0.67 − 0.00037 T2	15	0.0330	93.5
Positive Normal Stresses (σ+)
τ+ × H	τ+ = sqrt (−0.42 + 0.008 H2)	5	0.0027	96.6
	τ+ = sqrt (−0.11 + 0.007 H2)	10	0.0148	89.6
	-	15	0.0810 *	-
τ+ × λ	τ+ = sqrt (−1.55 + 0.0016 λ2)	5	0.0418	97.4
	τ+ = sqrt (−1.49 + 0.0021 λ2)	10	0.0024	95.6
	τ+ = sqrt (−1.38 + 0.0023 λ2)	15	0.0022	94.2
τ+ × T	τ+ = sqrt (0.054 + 0.00007 T2)	5	0.0409	92.0
	τ+ = sqrt (0.17 + 0.00029 T2)	10	0.0329	99.5
	τ+ = sqrt (0.41 + 0.00065 T2)	15	0.0303	99.6

* not significant with 95% confidence level.

). Slender and tapper had a statistical influence on the shear stress for all load levels (Table 8).

3.4. General Behavior of the Stresses along the Trunk

Due to the characteristic of the applied load, which consists of some forces and moments concentrated along the trunk, the stress/strain field in the model has some local maximum stresses/strains. Therefore, as the geometric parameters change, some of these local maximum stresses/strains can become the maximum global value. Therefore, although Table 4, Table 5 and Table 6show some sudden changes in the position of the maximum value, the stress/strain field in the model shows smooth variations due to changes in the geometric parameters.

3.5. Safety Factors

If safety factors below 1.0 were used as an indication of the risk of falling, one could conclude that a horizontal load of 10 kN would already be sufficient to break trees with H = 25 m in scenario 1, and trees with l greater to 40 in scenario 2 (

Table 9. Safety factors (relationship between load capacity and acting load) of normal and shear stresses for three horizontal load levels, three scenarios, and conditions (C) in each scenario.

Scenario	C	Force of 5 kN				Force of 10 kN				Force of 15 kN
		σ−	σ+	τ−	τ+	σ−	σ+	τ−	τ+	σ−	σ+	τ−	τ+
1	1	7.8	7.7	3.9	3.8	2.6	2.5	3.7	4.3	1.9	2.2	1.2	1.4
	2	24.8	29.9	13.3	13.7	9.1	8.9	4.0	6.1	4.0	6.1	4.0	4.6
	3	12.5	28.2	12.4	28.2	10.3	27.1	2.6	2.0	2.6	2.0	2.6	2.0
	4	8.4	10.5	8.4	10.5	7.9	10.5	1.0	1.3	1.0	1.3	1.0	1.3
	5	6.1	8.4	6.1	8.4	5.9	8.4	0.9	0.9	0.9	0.9	0.9	0.9
2	1	28.4	120	23.6	104	18.2	66.9	7.0	7.5	7.0	7.5	7.0	7.0
	2	31.8	56.6	16.1	26.2	10.8	17.1	7.0	7.5	7.0	7.5	4.9	6.7
	3	9.5	11.8	4.7	6.0	3.1	4.0	3.6	3.4	1.9	1.8	1.3	1.2
	4	4.8	4.5	2.4	2.2	1.6	1.5	2.0	2.3	1.0	1.1	0.7	0.8
	5	2.3	2.0	1.1	1.0	0.7	0.7	0.8	0.7	0.4	0.3	0.3	0.2
3	1	21.5	35.3	11.6	15.3	7.8	9.8	6.3	6.6	4.6	4.2	3.0	2.8
	2	20..4	16..1	10..4	8..3	6..9	5..6	7.0	8.0	4.0	3.9	2.6	2.6
	3	18.6	18.8	9.4	9.1	6.3	6.0	5.2	6.2	2.6	3.0	1.7	2.0
	4	14.0	9.7	6.9	4.9	4.6	3.2	4.0	4.1	1.9	2.1	1.3	1.4

Scenario 1: slenderness 22; taper 0.030; and height conditions: 5 m, 10 m, 15 m, 20 m, and 25 m. Scenario 2: tree height 7.5 m; taper 0.030; and slenderness conditions: 8, 16, 24, 32, and 40. Scenario 3: tree height 7.5 m; slenderness 22; and taper conditions: 0.010, 0.020, 0.030, and 0.050.

). In scenario 3, there was no broken tree with the tested load levels. The levels of wind speed that should correspond to these load levels depend on the size of the crown that we adopted in numerical simulations (Table 2).

Unlike wooden structures, trees can withstand maximum stresses without breaking [5,6,7,8,10,11]. In ref. [7] concluded that trees can reach stresses up to five times greater without breakage. The safety factor calculated for trees reached levels from 0.20 (5 times greater than a conventional wooden structure) to 0.33 (3 times greater than a conventional wooden structure) and is within the limit expected by the literature (3 to 5 times higher than that of structures). Considering the results, our simulation indicates that trees with l = 40 subject to 15 kN horizontal forces are at risk of falling.

4. Conclusions

-

The deflections’ behavior showed that the linearity with height is valid up to the point of application of the horizontal load, but from that point onwards, it depends on the diameter and, in some cases, also on the loading level.

-

Since the model considers the geometric variations of the section (deflections and diameter variations), the stresses show smooth variations along the trunk. However, the maximum values positive and negative, in module, are not equal and can undergo sudden variations in position along the trunk when local maximum stresses become global maximums.