Mechanical Stress Analysis of Tree Branches

Veiw figure View Figure

The component of the distributed weight that acts perpendicular to the branch and parallel to the cross sectional area is the shearing distributed load, w_s(x). The component of the distributed weight that acts parallel to the branch and perpendicular to the cross sectional area is the axial distributed load, w_a(x). The shearing distributed load contributes to the normal stress due to bending, while the axial distributed load contributes to the normal compressive stress on the branch. Both stresses were needed to be considered in order to calculate the total normal stress acting on a given cross section.

General equations were derived to calculate stress for each of the six cases. All of the equations derived assume a uniform MOE throughout the branch. To calculate the axial compressive stress acting on a cross section of a tree branch,the general equation for the axial distributed load,the load per unit length Equation (1) was

(1)

where m is the total mass of the branch, L is the total length of the branch, g is the gravitational constant, x is the distance from the tip of the branch to the where the stress analysis is desired, ρ is the density, V(x) is the volume of the branch up to the specified point x, A(x) is the branch cross sectional area and θ is the angle of the tree branch with respect to vertical. If the branch is perfectly horizontal then θ will be 90 degrees, making the axial distributed load zero. Using Equation (1) the axial force acting at any location from the tip of the branch, F_a(x),was determined from Equation (2) ^[4].

(2)

The normal axial stress due to the axial component of the weight is ^[4]

(3)

The normal stress due to bending is caused by the component of weight acting perpendicular to the branch. Shearing distributed load, w_s(x), which acts perpendicular to the length of branch was calculated using

(4)

The variables are the same as for the distributed axial load. Shear force V_s(x), is the component of weight that acts perpendicular to the length of the branch and parallel to the cross sectional area. Shear force was determined from the shearing distributed load ^[4].

(5)

The bending moment in terms of distance from the tip of the branch (Equation 6) was obtained by integrated the shear force, V_s(x), with respect to distance from the tip of the branch, x ^[4].

(6)

The normal stress due to bending was expressed as^[4]

(7)

where c(x) is the radial distance from the center of the cross section to the location in the radial direction where the stress analysis is desired. For instance, to obtain the bending stress at the top of the cross section c(x) would be the radial distance from the center of the cross section to the top of the cross section. I(x)is the moment of inertia of the cross section. Normal stress due to bending is compressive (negative) or tensile (positive) depending on the location of the cross section being examined. Tops of branches will be in tension and bottoms will be in compression. The neutral axis is where the stress passes from positive to negative and does not experience the effects of bending. The left and right sides of the cross section lie along the neutral axis and therefore have a bending stress of zero.

Total normal stress acting on a cross section is the sum of the stress due to bending and axial stresses. When adding these two stresses together the sign convention of each stress must be considered.

Equations (1) through (7) were general equations used to perform the stress analysis for each of the cases. For some of the cases additional analyses were needed. For case 1 the generic procedure described above can be followed exactly.

Case 2 was a fixed elliptical cross section. The equations described above were used by substituting in the cross sectional area of an ellipse, , and the moment of inertia for an ellipse, ^[4].

Cases 3 and 4 were tapered; therefore additional calculations have to be considered. Case 3 was tapered and had a circular cross sectional area (Figure 2).

The angle of taper was calculated using ^[1]

(8)

where R_ois the radius at the base of the branch. Rearranging for the radius of the branch in terms of distance from the tip of the branch, R(x), gives Equation (9) ^[1].

(9)

Figure 2. Case 3 tree branch model, tapered with a circular cross section. Where α is the angle of taper

Download as

Veiw figure View Figure

The cross sectional area in terms of distance from the tip of the branch then became ^[1]

(10)

The moment of inertia also varied with x and had to be calculated using R(x). By making these adjustments to the generic equations the stress analysis was performed for case 3.

Case 4 modeled the branch as a tapered elliptical beam. For case there were two angles of taper, one in the x-y plane and another in the x-z plane (Figure 3).

Figure 3. Case 4 taper in the x-y planes and x-z planes. R_Hoand R_Voare the radii in the horizontal and vertical directions respectively at the base of the tree branch.

is the angle of taper in the horizontal direction. β is the angle of taper in the vertical direction

Download as

Veiw figure View Figure

The angle of taper in the x-z plane was determined from the horizontal radius of the ellipse and was calculated using Equation (11).

(11)

Where ϕ is the angle of taper in the x-z plane, R_Ho is the horizontal radius at the base of the branch, and R_H(x) is the horizontal radius with respect to distance from the tip of the branch. Equation (12) provides the horizontal radius at any location, x, along the branch.

(12)

The angle of taper in the x-y plane was calculated using

(13)

Where β is the angle of taper in the x-y plane, R_Vo is the vertical radius at the base of the branch, and R_V(x) is the vertical radius with respect to distance from the tip of the branch.

(14)

Radii of the branch in the horizontal and vertical directions varied along the length of the branch; therefore, the cross sectional area of the branch needed to be expressed in terms of x, the distance from the tip of the branch as

(15)

The moment of inertia equation for the branch became

(16)

Substituting these equations into the generic equations the stress analysis was performed for case 4.

Case 5 was a fixed circular cross section with a non-uniform material. The MOE varied in the radial direction (Figure 4).

Figure 4. Case 5, a tree branch model of a fixed circular cross section and varying MOE in the radial direction

Download as

Veiw figure View Figure

The cross section of the branch was separated into concentric rings of equal width in which the MOE may vary. The derived equations assumed the branch was made of a uniform material. In order to use these equations the moment of inertia of each ring had to be transformed based on MOE values. Then each branch was analyzed as being composed of a uniform material and the generic equations for the stress analysis was used.

The Matlab program for case 5 was written so that the branch was subdivided into many rings. Therefore based on the outer diameter of the branch and the desired number of rings as inputs, the radius of each ring was calculated. The outer radius of each concentric ring is

(17)

where m represents the total number of rings in the cross section and h represents the number assigned to each ring. The inner-most ring is h=1, while for the outermost ring h=m. R represents the outer radius of the branch. The area of each ring was calculated.

The innermost ring had the smallest MOE value and therefore its moment of inertia was not transformed. All the other rings hadlarger MOE values than the innermost ring. To account for larger MOE values, a ratio of each ring’s MOE value to the MOE value of the innermost ring were used to transform the moment of inertia of each section. For example, if a ring had a MOE that was two times greater than that of the innermost ring, then the moment of inertia of the outer ring was needed to be two times greater as well. After transforming the moment of inertia values, the branch was analyzed as having a uniform material. The ratio of the MOE of an outer ring to the MOE of the innermost section was calculated as ^[4]

(18)

Note for h=1 the ratio became 1. The transformed moment of inertia became ^[4]

(19)

The new moment of inertia of the entire cross section is the sum of the transformed moments of inertia for each ring.

(20)

The stress due to the bending moment on each ring was calculated using the new moment of inertia and the MOE ratio ^[4].

(21)

The axial stress on each ring was calculated by considering the various MOE values of each ring. First the distributed axial load, w_a(x), and the total axial force, F_a(x), were calculated using the total area of the cross section, A_total(x). These variables are independent of the MOE of the branch.

The area and MOE of each ring was considered in all calculations of the axial force acting on its section. The sum of the axial forces acting on each ring was equivalent to the total axial force (Equation 22).

(22)

where F_a(x) is the total axial force due to the weight of the branch. Forces F₁(x) through F_m(x) were the compressive forces acting on each ring. The forces on each ring caused the branch to deform along the axis of the branch, making the branch shorter. This deformation was calculated as ^[4]

(23)

where δ_h is the deformation of the ring, F_his the axial force acting on the ring, A_h is the area of the ring, and E_his the MOE of the ring. Since all of the rings were connected to one another they allhad the same deformation.

(24)

Using equation 24, the axial force on each ring was determined. The axial force on any of the rings in terms of the axial force on ring 1 is

(25)

The total axial force on the branch in terms of the force on ring 1 became

(26)

Equation 26was rearranged to solve for the axial force acting on ring 1 (Equation 27).

(27)

This procedure was repeated to determine the axial forces acting on each ring. Using the axial force acting on each ring the axial stress on each ring was calculated as

(28)

The total stress for case 5 is the sum of the axial and bending stresses.

Figure 5. Case 6: Curvy tree branch model, where y_t is the height of the branch on top of the axis,y_b is the height below the axis

Download as

Veiw figure View Figure

Case 6 accounted for the curves of a branch in the x-y plane, seen from the side view. The curviness of the branch was determined by measuring the height of the branch above and below the axis of the branch at incremental distances along its length. The axis line was the line that connects the center of the cross section at the base of the branch to the center of the cross section at the tip of the branch (Figure 5).

The distance from the axis to the top of the branch is y_t. The distance from the axis to the bottom of the branch is y_b_. The diameter of the branch is the sum of y_t and y_b. The branch was assumed to have a circular cross section in this case because the varying heights of the branch were only considered in the x-y plane (side view), not the x-z plane (top view).

To perform a stress analysis for case 6, stresses were calculated at the base of the branch, where x=L. At this location the axis of the branch coincided with the center of the branch cross section. The axis was defined earlier as the line that connects the centers of the cross sections at the base of the branch and at the tip of the branch. Performing the stress analysis at the base provided the maximum stresses experienced by the branch.

To account for the curviness of the branch a best fit polynomial equation was determined for y_t and y_bin terms of x. Because the branch is curved, the centroid of the branch may not lie on the branch axis. Therefore the axial component of weight applied at the centroid must be moved to the axis to perform the stress analysis. To account for moving the forces an additional bending moment must be added.

The centroid of the area in the x-y plane was determined by first calculating the centroid of the areas above and below the axis separately. Then the two centroids were combined to obtain the centroid of the entire branch. The areas of the top and bottom parts of the branch must be calculated separately, by integrating the y_t and y_bequations to find the area under the curve. Using the areas the x-coordinate of the centroid for the top (Equation 29) and bottom (equation 30) portions of the branch were calculated using ^[4]

(29)

(30)

where y_t and y_b are both functions of x, y_i=f_i(x). To determine the y-coordinate of the centroid, two additional best fit polynomial equations were determined for x in terms of y_t and y_b, giving x_i=f_i(y). The y-coordinates of the centroids above (Equation 31) and below (Equation 32) to the axis line are ^[4]

(31)

(32)

Next the x and y coordinates of centroid, relative to the branch axis, were calculated for the entire branch (Equations 33 and 34) ^[4].

(33)

(34)

To perform the stress analysis, the axial force, F_a, was moved from the centroid to the center of the branch’s cross section. The perpendicular distance that F_a, must be moved was. Therefore the additional bending moment that must be added to account for moving the axial force was

(35)

The bending moment that was due to the shear force (Equation 36) was the product of the shear force and the perpendicular distance between the x-coordinate for the centroid and the base of the branch.

(36)

Recall that was measured from the tip of the branch, therefore it was subtracted by the length to get the distance between the base and the centroid.

The total bending moment was the sum of these two bending moments where their signs were considered. After the total bending moment was calculated, the stress due to bending was determined. The axial stress was calculated using Equation 3 after the force was moved to the center of the cross section. The total normal stress acting on the cross section at the base was the sum of normal stress due to bending and the axial stress in which the signs were considered.

Figure 6. Image of the Graphical User interface

Download as

Veiw figure View Figure

The equations derived for each of the cases were used to write codes in Matlab that can perform a stress analysis on any branch. The Matlab code provided an analytical solution for the stress analysis. To make the code more user friendly a Graphical User Interface (GUI) was generated to determine the stresses (Figure 6).

2.2. Finite Element Analysis

Using an approximate solution would allow branches with a wide variety of geometries to be analyzed, the software Abaqus was selected to perform the stress analysis so that branches with various geometries could be modeled. Abaqus uses Finite element analysis (FEA) to compute an approximate solution for stress. Therefore the method selected to model the branch in the program would have an effect on the final stress calculated values. Two separate methods for modeling tree branches were considered, a 2-D wireframe and a 3-D wireframe. Stress values for each option were compared.

A simple branch with a constant cross section, and material was modeled using each of the two methods(2D and 3D wireframe) described above in Abaqus. For each branch, the base end was fixed; preventing deflection and rotation in all directions, and gravity was applied. Once the model was completed, Abaqus performed the stress analysis and provided a contour image displaying stress, strain or displacement throughout the entire branch (Figure 7 and Figure 8).

Figure 7. 2D “Simple Beam” branch model

Download as

Veiw figure View Figure

Figure 8. 2D Model containing secondary branches

Download as

Veiw figure View Figure

The results from Abaqus and our analytical solution were identical (zero error). The 2D wireframe and the 3D wireframe both had the same approximate solutions. We selected 3D wireframe for all subsequent simulations. The 2D wireframe alternative was possible but the 3D wireframe was more appropriate since branches are three dimensional.

After verifying FE simulation, a more comprehensive simulation was developed. The model simulated 3D geometry for branches with side branches and considered the presence of leaves and fruits. Each branch was modeled in Abaqus as a 3D wire frame. The geometry was modeled as connected datum points with 3D wires. Each branch segment (between adjacent points) was assigned a material and a diameter. The material assignment was based on whether the segment was categorized as old, medium or new growth. The diameter of the segment was the average diameter of the two ends of the segment.

After the branch was created, loads and a boundary condition were applied to the branch. A load was applied to each point where a leaf or a fruit was present, and a gravitational force was applied over the entire branch. As stated above, a boundary condition was applied to the base of the branch. The boundary condition was fixed, preventing deflection and rotation in any direction. The stress analysis was then performed. Figure 9 shows a sample simulation on a real tree branch.

Figure 9. 3D Finite Element analysis on a real tree branch (Sassafras albidium)

Download as

Veiw figure View Figure

Once the model was created, the GUI is developed which required data of geometry and mechanical properties within an excel file and outputted an executable file in Abaqus. Thereafter, the program was initiated in Abaqus to perform stress analysis. This process enables users, even with limited knowledge on using Abaqus, to perform a 3D stress analysis on tree branches.

Next, the developed analytical analysis on the 6 cases examined trends. For cases 1 through 5 all of the input variables were kept constant except for one in order to examine how changing that single variable affects the results. The variables considered included branch length, diameter, density and aspect (or angle) with respect to the vertical. For case 5 changes in MOE were also considered. We also analyzed the results from FE model.

3. Results and Discussion

For each of the cases tested, all of the cases produced similar trends. With an increase in branch length, mass increased since density was constant. However the ratio of length to mass remained constant since they both increased at the same rate. Therefore, as branch length increased, maximum stress increased (Figure 10).

As the length increased, maximum stress increased at a greater than linear rate possibly due to an increased bending arm. As for the tree branch models, cases 1, 2, and 5 (non-tapered branches), all produced nearly identical results. Note that the data points for case 1 in Figure 10 coincide with those of case 2. According to the stress analysis, the benefits of elliptical model versus the circular cross sections were negligible. Tapered branches (cases 4 and 5) had smaller stress values than non-tapered branches. As branch length increased, differences between tapered and non-tapered branches were apparent. Non-Tapered branches were larger, with larger stress values.

Figure 10. Length of branch versus maximum stress

Download as

Veiw figure View Figure

When branch diameters were varied among samples, masses also varied since densities were constant. However, diameters and masses did not increase at the same rates. When examining changes in diameters solely, increases in diameters produced decreases in stress (Figure 11). However, increases in diameters per unit mass lead to higher stresses (Figure 12).

Figure 11. Diameter versus stress

Download as

Veiw figure View Figure

Figure 12. Diameter to mass ratio versus maximum stress

Download as

Veiw figure View Figure

In the above figures cases 1 and 2 yielded nearly identical results. For all of the cases as branch diameters increased, stresses decreased. Therefore, when branch mass was not considered, smaller diameter branches had less stress. Among the branches tested, larger diameter branches had more stress. Relationships among tree branch models were the same in which the tapered models had less stress than the fixed models, and the shape of the cross section had a negligible effect on the stress. Thicker branches resisted more stress.

Figure 13. Mass versus maximum stress

Download as

Veiw figure View Figure

Figure 14. Aspect from vertical versus maximum stress

Download as

Veiw figure View Figure

For each of the cases the branch density and branch aspect relative to vertical were also analyzed. The density was varied by keeping the dimensions of the branch constant and varying mass. As branch density increased, the stress increased (Figure 13) and as branch aspect increased (became more horizontal), stress increased (Figure 14).

There were noticeable differences among the various cases. Cases 1 and 5, the two cases with fixed circular cross sections, produced higher stress than the other cases. Cases 1 and 5 also produced very similar results. Case 2 had lower stress values than case 1, which indicated that branch aspect and mass the elliptical shape reduced stress. This result was further supported by the significantly lower stress values produced by case 4, the tapered elliptical branch compared with case 3, the tapered circular branch. Again the tapered branches had less stress than non-tapered branches.

For case 5, the effects of varying MOE values were analyzed. Case 5 was divided into three rings for this analysis, and MOE values were selected for each ring based on results from previous research (unpublished results). The percent increase in MOE values between the first and second rings was 113%, while the percent increase between the second and third rings was 25%. To study the effect of change of MOE in overall stress, the MOE of all three rings was increased or decreased by 5% until there was either an additional positive or negative 20% from the reference percent (Figure 15).

Figure 15. Percent change in MOE from reference percent versus maximum stress

Download as

Veiw figure View Figure

As the percent change in MOE increased between concentric rings, stress increased. Thus, branches composed of materials with more uniform MOE throughout experienced less stress than branches with large differences of MOE (in the radial direction).

Case 6 accounted for the curviness of branches. Various dimensions were entered into the program to create branches with varying degrees of curves. For consistency all of the branches had the same base diameter and the same tip diameter. First, stress on a branch with no curves was compared to case 3. Both should produce similar results because they both had the same base diameter and were very narrow at the branch tips. However, there was a 33% difference between the two branches (Table 1).

Table 1. Comparison of the stresses in a non-curved branch as calculated by cases 3 and 6. The percent error calculation suggests that case 3 is more accurate

Download as

View current table in a new window

Tables index

Veiw figure View Table

The differences in stress values may result from the use of a best fit polynomial equation even when the branch may be more accurately represented by a straight line. Based on these results case 6 should only be used when the branch is curved, or when there is an abnormality in the branch which would shift the centroid of the branch away from the branch axis.

Figure 16 shows the shapes of two branches developed by the best fit polynomial equations and their corresponding maximum stress values.

Figure 16. Branch Shapes generated by polynomial equations and their corresponding stress values

Download as

Veiw figure View Figure

The tree branch generated in Figure 16a has small curves. The stress value calculated for this branch more closely aligns to the stress calculated using case 3 displayed in Table 1. The branch in Figure 16a has small curves so it resembles case 3. However, since curves were present, a polynomial a polynomial equation provided a better relationship than a linear equation. Figure 16b has more prominent curves, and produced a larger stress value. The Matlab program was written so that the polynomial fit equation will pass through all of the input points. Therefore, whenever a branch has a drastic change in height all changes should be entered. In addition, input points may be added at locations between changes in order to guide the path of the line. The FE analysis was also used to study the stress patterns on main branches with and without leaves. Figure 17 shows that while stresses increased when leaves were present, the overall stress pattern was the same as when leaves were not present.

Figure 17. Stress pattern with and without leaves

Download as

Veiw figure View Figure