this is for holding javascript data
Nicholas Davies edited method.tex
over 7 years ago
Commit id: 89368a122dc6a215d17679a768bc059e49ea02a9
deletions | additions
diff --git a/method.tex b/method.tex
index 5318809..98d6ac5 100644
--- a/method.tex
+++ b/method.tex
...
\section{Method}
Stems were approximated as truncated cone meshes, represented using quadratic tetrahedral elements. Each (fully grown) mesh was divided into elements covering one fifth the radius, one fifth the height and one sixteenth of the circumference. Data presented in \citet{waghorn_influence_2007} was used to create a regression to estimate stem radius (\(r_s\)) from stocking rates based on the slenderness ratios. The regression presented in Equation \ref{eq:slenderness} was calculated from reported slenderness ratios for stocking rates of between 275 and 1457 stems per hectare and assuming a stem height of 15 m.
\begin{equation}\label{eq:slenderness}
r_s = e^{−0.000328S−1.868}
\end{equation}
Where \(S\) is the stocking rate, assumed to be either open grown, one stem per hectare or 741 stems per hectare.
As trees grow, the newest cells are superimposed on the existing structure, creating a time evolving stress state. In order to accommodate these time dependent changes, the mesh needs to be 'grown', here time is treated in discrete steps. Wind effects are not accounted for during growth, however, self weight and growth stress are. Growth stresses develop as cells divide in the cambium, resulting in the new cells being in a non-stressed state on the surface of a pre-stressed structure. During secondary wall formation there is a tendency for cells to contract longitudinally, for more details see \citep{archer_growth_1986}. As the pre-existing structure restrains further contraction a stress profile of longitudinal compression at the centre of the stem and tension at the outer edge forms.
The pre-stressed state of the cambium from gravity and growth stresses produces the need to successively ’grow’ the mesh which is to be used so that each mesh addition is added in a non stressed state on top of a stressed surface. Mesh growth is achieved by first defining an initial state, this state can be thought of as a seedling. The seedling is subjected to gravitational forces, deformation occurs and an updated mesh is created. The positions of the nodes and vertices of the deformed mesh are used to calculate the positions of new vertices and nodes to be added in order to represent cambial growth in a non stressed state onto the stressed surface along with apical growth to a predefined height. This procedure is repeated for five time steps imposing gravity (and in required cases growth stress), regardless of the final age, i.e. a tree with wind imposed at 5 years utilises 1 year steps, a tree with wind imposed at 15 years has 3 year time steps.
The material constants for \textit{Pinus radiata} reported by \citet{Davies_2016} were used to parametrise the model. Four sets of constants representing extreme values of green density and stiffness for the species were used. Note the abbreviations of the sample type used here, the first letter defines the density as either \(H\) (High density) or \(L\) (Low density) and the second letter defines the MFA as either \(W\) (Wide) or \(S\) (shallow).
Linear interpolation was used to approximate material constants between two samples to provide a gradient of material properties from the corewood to the outerwood. The TRP of cellular properties was investigated using the gradient of elastic constants and Poisson ratios from different samples by applying simulated wind loads to stems with different radial profiles.
While the TRP of \textit{Pinus radiata} usually follows a non-linear increase for density and a non-linear decrease for MFA \citep{burdon_juvenile_2004}, during the first 15 years the trend can be approximated as linear, after this however the non-linear behaviour becomes more prevalent as the properties stabilise. It should be noted that there is much variation between individuals and this is not universal. The TRP was investigated by calculating the appropriate elastic and shear moduli along with Poisson ratios for each point in the stem. The linear interpolation was calculated between two samples chosen to represent corewood and outerwood, Equation \ref{eq:constants} shows the general form.
\begin{equation}\label{eq:constants}
\lambda = \frac{-r}{r_t}(\lambda_i-\lambda_o)+\lambda_i
\end{equation}
Where \(r\) is the stem radius at the current point, \(r_t\) is the total radius at 15 years, \(λ_i\) and \(λ_o\) are the values of the material properties obtained from \citet{Davies_2016} at \(r = 0\) and \(r = r_t\) respectively.
Due to the micro-structure of wood the native coordinate system for describing material properties is not the same as the global coordinates system used to impose external forces on the stem. Transformation between the two systems is needed. The \(rtl\) local system used for experimental work presented in \citet{Davies_2016} is used in combination with the interpolations presented in Equation \ref{eq:constants} and presented as the stiffness matrix in Equation \ref{eq:stiffness_mat} using Voigt (engineering) notation at any point in the stem. Because the stiffness matrix is calculated in the local coordinates it needs to be converted into an \(xyz\) system in order to apply wind loadings in a sensible fashion, Equation \ref{eq:trans_mat} provides the transformation matrix. The use of Voigt notation allows for transformation of the elastic constants into the global system, via Equations \ref{eq:cm_trans} to \ref{eq:stress_tans}. It is assumed that there is no spiral grain occurrence within the stem, the local \(l\) axis is always parallel with the global \(z\) axis. Material constants in the longitudinal direction are parallel with the \(z\) axis, i.e. there is no correction angle applied to account for taper.
\begin{equation}
S_l = \begin{bmatrix}
\frac{1}{E_r}& -\frac{v_{tr}}{E_t}& -\frac{v_{lr}}{E_l}& 0& 0& 0\\
-\frac{v_{rt}}{E_r}& \frac{1}{E_t}& -\frac{v_{lt}}{E_l}& 0& 0& 0\\
-\frac{v_{rl}}{E_r}& -\frac{v_{tl}}{E_t}& \frac{1}{E_l}& 0& 0& 0\\
0& 0& 0& \frac{1}{2G_{tl}}& 0& 0\\
0& 0& 0& 0& \frac{1}{2G_{lr}}& 0\\
0& 0& 0& 0& 0& \frac{1}{2G_{rt}}
\end{bmatrix}
\end{equation}
\begin{equation}\label{eq:stiffness_mat}
C = S^{-1}
\end{equation}
\begin{equation}\label{eq:trans_mat}
G = \begin{bmatrix}
a_r^xa_r^x& a_r^ya_r^y& a_r^za_r^z& a_r^ya_r^z& a_r^xa_r^z& a_r^xa_r^y\\
a_t^xa_t^x& a_t^ya_t^y& a_t^za_t^z& a_t^ya_t^z& a_t^xa_t^z& a_t^xa_t^y\\
a_l^xa_l^x& a_l^ya_l^y& a_l^za_l^z& a_l^ya_l^z& a_l^xa_l^z& a_l^xa_l^y\\
2a_t^xa_l^x&2a_t^ya_l^y &2a_t^za_l^z &a_t^ya_l^z+a_t^za_l^y &a_t^xa_l^z+a_t^za_l^x & a_t^xa_l^y+a_t^ya_l^x\\
2a_r^xa_l^x&2_r^ya_l^y &2a_r^za_l^z &a_r^ya_l^z+a_r^za_l^y &a_r^xa_l^z+a_r^za_l^x &a_r^xa_l^y+a_r^ya_l^x \\
2a_r^xa_t^x&2a_r^ya_t^y &2a_r^za_t^z & a_r^ya_t^z+a_r^za_t^y &a_r^xa_t^z+a_r^za_t^x & a_r^xa_t^y+a_r^ya_t^x
\end{bmatrix}
\end{equation}
Where \(a^j_i\) is the directional cosine from \(j\) to \(i\) (global to local), \(j = x, y, z\) and \(i = r, t, l\).
\begin{equation}\label{eq:cm_trans}
\boldsymbol{C = G^TC_lG}
\end{equation}
\begin{equation}
\boldsymbol{G^T = D^{-1}}
\end{equation}
\begin{equation}\label{eq:stress_tans}
\boldsymbol{\sigma_l = D\sigma}
\end{equation}
The stiffness matrix \(\boldsymbol{C}\), stress vector \(\boldsymbol{\sigma}\) and strain vector \(\boldsymbol{\epsilon}\) are in the global coordinate system while \(\boldsymbol{C_l}\),\(\boldsymbol{\sigma_l}\),\(\boldsymbol{\epsilon_l}\) are in the local system.\\
The forces being applied in the global system (through the transformation of the elastic constants) results in stresses being calculated in the global system. As the experimental work was conducted in the local system the proportionality limit stresses are in the local system. Equation \ref{eq:stress_tans} was used in order to transform the stresses in the global system (calculated by the model) back into the local system in order to evaluate the state of failure.
The forces due to gravity from self weight of the stem and canopy were applied to appropriate domains. The force due to gravity from self weight of the stem was applied as a body force to the whole domain \(\omega\) (i.e. the whole stem) in Equation \ref{eq:Bs}. Green density \(\rho\) is calculated at each point using Equation \ref{eq:constants} and g is gravity (\(9.81 m/s \)). The canopy is assumed to take the geometry of the upper half of an ellipsoid described in Equation \ref{eq:Vc}, where \(z_0\) is the height at the bottom of the current element and \(z_1\) is the height at the top of the current element, with the top of the ellipsoid at the top of the stem, height \(h\). The start height \(S_c\) of the canopy was estimated using Equation \ref{eq:Sc} which was derived from data presented by \citet{waghorn_influence_2007}. In order to estimate crown radius, \(r_c\) , Equation \ref{eq:rc} is used, which was derived from the assumption that the maximum radius which can be achieved by the crown is half the distance between two trees evenly spaced in the stand. The weight of the canopy on the stem is applied as a body force described in Equation \ref{eq:Bc1}.
\begin{equation}\label{eq:Sc}
S_c = −123.4r_s + 20.6
\end{equation}
\begin{equation}\label{eq:rc}
r_c = 34.4r_s − 2.07
\end{equation}
\begin{equation}\label{eq:Bs}
B_s = \rho g
\end{equation}
\begin{equation}\label{eq:Vc}
V_c = \pi r^2_c \Bigg(z_1 \Bigg( \frac{1-z^2_1}{3(h-S_c)^2}\Bigg) - z_0\Bigg(\frac{1-z^2_0}{3(h-S_c)^2}\Bigg)\Bigg)
\end{equation}
\begin{equation}\label{eq:Vs}
V_s = \frac{\pi}{3}\Big(r^2_{\|_{z_0}}z_0 - r^2_{\|_{z_1}}z_1\Big)
\end{equation}
\begin{equation}\label{eq:Bc1}
B_c|_{\boldsymbol{\Omega_1}} = V_c \rho_c g \frac{1}{V_s}
\end{equation}
\begin{equation}\label{eq:O1}
\boldsymbol{\Omega_1} = \Big\{(x,y,z) \in \boldsymbol{\Omega} : z > S_c \Big\}
\end{equation}
Where \(r_s\) is the radius of the stem (from Equation \ref{eq:slenderness}), \(r_{z_0}\) is the radius of the stem at height \(z_0\) , \(r_{z_1}\) is the radius of the stem at height \(z_1\) and \(g\) is gravity. \(\rho_c\) is the canopy density of \(5.6 kg/m^3\) estimated from data in \citet{Beets_1996} at a stocking of 741 stems per hectare. The canopy force due to gravity is only applied to the sub-domain \(\boldsymbol{\Omega_1}\) , defined in Equation \ref{eq:O1}, where \(z\) is the vertical coordinate of any point. \(\frac{1}{V_s}\) is needed in order transform the canopy’s gravitational force into a force per unit stem volume.
In order to stress the stem, a constant wind profile was applied to the canopy. The crown sail area was assumed to be the upper half of an ellipse attached to the stem on the surface \(\Gamma_1\) (defined by Equation \ref{eq:T1}) a surface subregion of total surface \(\Gamma\). The common drag model presented in Equation \ref{eq:Tw} has been used previously \citep{spatz_basic_2000, rudnicki_wind_2004, mayer_windthrow_1989} and is used to approximate the wind load. It should be noted that more complex models are available \citep{coutts_wind_1995}. The drag coefficient \(\varsigma\) in Equation \ref{eq:DC} was produced from data reported by \citet{Mayhead_1973}, for Scotts pine as no data was available for radiata. The use of the \citet{Mayhead_1973} Scotts pine data set has previously been suggested as a suitable substitute \citep{moore2001relative}.
\begin{equation}\label{eq:DC}
\varsigma = e^{-0.377\omega-0.306}
\end{equation}
\begin{equation}\label{eq:Ac}
A_c = \int_{z_0}^{z_1}2r_c\sqrt{1-\frac{z^2}{(h-S_c)^2}}dz
\end{equation}
\begin{equation}\label{eq:As}
A_s = \frac{\pi}{2}\Big( r_{\|_{z_0}}\sqrt{z_0^2 + r_{\|_{z_0}}^2} - r_{\|_{z_1}}\sqrt{z_1^2 + r_{\|_{z_1}}^2} \Big)
\end{equation}
\begin{equation}\label{eq:Tw}
T_\omega |_{\boldsymbol{\Gamma_1}} = \frac{1}{2}\rho_{air}\varsigma \omega^2A_c \frac{1}{A_s}
\end{equation}
\begin{equation}\label{eq:T1}
\boldsymbol{\Gamma_1} = \{ (x,y,z) \in \boldsymbol{\Gamma}: x > 0, z > S_c \}
\end{equation}
\(T_w\) is the force induced on the stem via the canopy for a given wind speed \(\omega\). Air density \(\rho_{air}\) air is constant at 1.226 \(kg/m^3\) \citep{Mayhead_1973}. The canopy area \(A_c\) is calculated as per Equation \ref{eq:Ac} and stem area \(A_s\) is calculated as per Equation \ref{eq:As}. \(\frac{1}{A_s}\) is needed in order to transform the wind induced force into a force per unit stem surface area. \(T_{\omega}\) is imposed on the stem as a boundary force in the sub-domain \(\Gamma_1\).
\citet{papesch_mechanical_1997} reported statistical regressions (reproduced in Equations \ref{eq:Mc} and \ref{eq:DMc}) in order to predict the maximum bending moment and the angle of deflection at the maximum applied bending moment. Assuming the deflection when a stem first reaches its proportionality limit stress coincides with the angle of deflection at the maximum bending moment, by calculating the expected deflection at the maximum bending moment and comparing with the results the model produces gives insight into the general accuracy of the model. The force imposed by the canopy can also be converted into a bending moment at the first wind speed which breaks proportionality and compared to Equation \ref{eq:Mc}’s prediction for further insight.
\begin{equation}\label{eq:Mc}
ln(M_c) = 2.5578ln(h) − 3.18
\end{equation}
\begin{equation}\label{eq:DMc}
D_{M_c} = −0.5416h + 21.099
\end{equation}
Where \(M_c\) is the maximum bending moment, \(h\) is the height and \(D\) is deflection.
Growth stresses are represented in a simplified fashion. \citet{gillis_elastic_1979} provide Equation \ref{eq:GS} to describe growth stress profiles. The growth stress profile presented in Equations \ref{eq:rcore} to \ref{eq:GS} are imposed on the stem at each growth step, and the new growth added to the pre-stressed stem. Growth stresses are imposed by adding the growth stress vector directly onto the stress vector during its calculation (in a similar way to how temperature dependent stresses are often represented) the growth stress vector, \(\boldsymbol{\sigma_{gs}}\) is presented in Equation \ref{eq:stress_gs}. It should be noted that the growth stresses are only added into the longitudinal direction, it is assumed that the resulting deformation through the interaction of the stiffness matrix will result in an appropriate stress profile in the other material directions. At the periphery of the stem in a 15 year old tree the tensile growth stresses calculated here ranges from \(0.8\) to \(2 MP_a\) depending on the outerwood used, this is similar to the values reported by \citet{timell_compression_1986-3} for \textit{Pinus taeda} of \(1.2 MP_a\) in tension at the periphery, however they also report similar values for compression at the centre of the stem, which is somewhat lower than the assumption used here.
\(r_{core}\) is the radius which the growth stress model assumes the stress strain relationship is flat at the limit of proportionality under compression. \(Y_l\) is the yield limit set to \(0.1\) as presented by \citet{gillis_elastic_1979} and is the radial proportion of the stem at the compressive proportionality limit. \(R\) is the maximum radius of the stem for the height of the point being evaluated. When the point being evaluated is the apex, \(R\) is set to \(0\), as the apical growth from the current time step is assumed not to have gone through the maturation process.
\begin{equation}\label{eq:rcore}
r_{core} = Y_lR
\end{equation}
\begin{equation}\label{eq:ITstress}
IT_{stress} = \frac{Y_{cs}}{\frac{R}{r_{core}}\big(1 - \frac{r_{core}}{R}\big)}
\end{equation}
\(IT_{stress}\) is the initial longitudinal tensile stress and \(Y_{cs}\) is the compressive stress at the limit of proportionality calculated as per Equation \ref{eq:constants}.
\begin{equation}\label{eq:GS}
G_s = \begin{cases}
0,& \text{if } r_{core} = 0\\
-IT_{stress}\frac{R}{r_{core}}(1-\frac{r_{core}}{R}),& \text{else if } r \leq r_{core} \\
-IT_{stress}\frac{R}{r_{core}}\frac{(\frac{r_{core}}{r}-2\frac{r_{core}}{R}+\frac{r_{core}^2}{R^2})}{1-\frac{r_{core}}{R}}, & \text{otherwise}
\end{cases}
\end{equation}
\begin{equation}\label{eq:stress_gs}
\boldsymbol{\sigma_{gs}} = \begin{bmatrix}
0\\0\\G_s\\0\\0\\0
\end{bmatrix}
\end{equation}
Hooks law can be used to characterise elastic relationships in mathematical terms. Generalised Hook’s law (Equation \ref{eq:hooks}) is used here to describe the mechanical elastic characteristics of wood. Because the characterisation of the proportionality limit is the point at which the stress strain curve is no longer linear, stress is only proportional to the strain and not the strain rate. Above this point non-linear elastic and plastic effects need to be considered, which was beyond the scope of this study.
Stress, \(\boldsymbol{\sigma}\) is calculated from strain \(\boldsymbol{\epsilon}\) and the stiffness matrix \(\boldsymbol{C}\). When growth stresses were not considered Equation \ref{eq:hooks} was used, when growth stresses were considered \(\boldsymbol{\sigma_{gs}}\) was defined by Equation \ref{eq:stress_gs} and \(\boldsymbol{\sigma}\) was calculated through Equation \ref{eq:stress}.
\begin{equation}\label{eq:E}
\boldsymbol{\epsilon} = \frac{1}{2}(\bigtriangledown \boldsymbol{u}+\bigtriangledown \boldsymbol{u}^T)
\end{equation}
\begin{equation}\label{eq:hooks}
\boldsymbol{\sigma = C\epsilon}
\end{equation}
\begin{equation}\label{eq:stress}
\boldsymbol{\sigma = C\epsilon+\sigma_{gs}}
\end{equation}
Strain energy density can then be calculated
\begin{equation}\label{eq:W}
W = \frac{1}{2} \boldsymbol{\sigma\epsilon}
\end{equation}
and the total potential energy found
\begin{equation}\label{eq:psi}
\boldsymbol{\prod} = \boldsymbol{\int_\Omega} W \mathit{d}\boldsymbol{\Omega} - \boldsymbol{\int_\Omega B_s.u} \mathit{d}\boldsymbol{\Omega} - \boldsymbol{\int_{\Omega_1} B_c.u }\mathit{d}\boldsymbol{\Omega_1} - \boldsymbol{\int_{\Gamma_1}T_s.u} \mathit{d}\boldsymbol{\Gamma_1}
\end{equation}
\begin{equation*}
\boldsymbol{\Omega_1} = \{(x,y,z) \in \boldsymbol{\Omega}:z>S_c\}
\end{equation*}
\begin{equation*}
\boldsymbol{\Gamma_1} = \{(x,y,z) \in \boldsymbol{\Gamma}:x>0,z>S_c\}
\end{equation*}
By taking the directional derivative of \(\prod\) with respect to the change in \(u\) and setting it to zero the displacement field \(u\) can be calculated at the minimum potential energy.
\begin{equation}\label{eq:F_new}
F = \boldsymbol{\bigtriangledown_u} \prod(\boldsymbol{u}) = 0
\end{equation}
Subject to the Dirichlet boundary conditions, equations \ref{eq:BC1} and \ref{eq:BC2}.
\begin{equation}\label{eq:BC1}
\boldsymbol{u} |_{\boldsymbol{\Omega_{bc1}}} = \begin{bmatrix}
u_x\\
u_y\\
0
\end{bmatrix} \hspace{0.2cm}
\end{equation}
\begin{equation*}
\boldsymbol{\Omega_{bc1}} = \{(x,y,z) \in \boldsymbol{\Omega}:z>tol_z\}
\end{equation*}
and
\begin{equation}\label{eq:BC2}
\boldsymbol{u} |_{\boldsymbol{\Omega_{bc2}}} = \begin{bmatrix}
0\\
0\\
0
\end{bmatrix} \hspace{0.2cm}
\end{equation}
\begin{equation*}
\boldsymbol{\Omega_{bc2}} = \{(x,y,z) \in \boldsymbol{\Omega_{bc1}}:(x,y)-tol\}
\end{equation*}
Where \(tol_z\) is the tolerance required to select all points at the boundary between the base of the stem and the ground, and where \(tol\) is the tolerance sufficient to select only the points at the boundary between the base of the stem and the ground for the initial seedling's radius.\\
By solving Equation \ref{eq:F_new} for displacement \(\boldsymbol{u}\) subject to the boundary conditions (Equations \ref{eq:BC1} and \ref{eq:BC2}) and internal and external loadings (Equations \ref{eq:Bs}, \ref{eq:Bc1}, \ref{eq:Tw} and \ref{eq:GS}) the displacement field was obtained. Hooks law (Equation \ref{eq:hooks}) was used to find the stress. Solving of the system is achieved through the FEM.
Failure surfaces were created using Tsai and Wu’s failure criterion \citep{tsai_general_1971} and calculated for every point in the stem for all wind speeds. Each point is evaluated for its safety factor, where a factor of one is on the failure surface, with a lower than one factor being before the limit of proportionality and higher than one factor after the limit. The values are the observed stress over the proportional limit stress given the other five stress states. Due to the dependence of each direction of failure on the other directions, all but the variable in question are held constant at their modelled values. Once passed the proportional limit the linear stress strain curves were still assumed, this is not physically accurate. The proportional limit stresses were calculated from the linear interpolation described by Equation \ref{eq:constants}
From the experimental work presented in \citet{Davies_2016} and the relationships presented in Equation \ref{eq:constants} proportional limit surfaces for each point were calculated using \citet{tsai_general_1971}’s criterion. Because the proportional limit surfaces are defined in the local \(rtl\) coordinate system and the stresses provided by the model are in the global \(xyz\) system Equation \ref{eq:stress_tans} was used to convert each stress vector into the local system at its given global location. The failure criterion was applied through Equation \ref{eq:tw_app} and the maximum and minimum stress values which could be obtained without failure calculated for each point assuming all other stresses stay fixed (Equations \ref{eq:tw_opt_max} and \ref{eq:tw_opt_min}). With the maximum tensile and compressive stresses a safety factor was calculated as per Equations \ref{eq:SFten} and \ref{eq:SFcomp}. The stress bounds can be used to investigate how much redundant strength is available at failure by both position in the stem and direction of stress.
\begin{equation}\label{eq:tw_app}
\boldsymbol{\sigma^T q + \sigma^T P \sigma} - 1 = 0
\end{equation}
\begin{equation}\label{eq:tw_opt_max}
\boldsymbol{\sigma_{max}} = max(\sigma_i) \hspace{0.2cm} \text{subject to} \hspace{0.2cm} \boldsymbol{\sigma^T q + \sigma^T P \sigma} - 1 = 0
\end{equation}
\begin{equation}\label{eq:tw_opt_min}
\boldsymbol{\sigma_{min}} = min(\sigma_i) \hspace{0.2cm} \text{subject to} \hspace{0.2cm} \boldsymbol{\sigma^T q + \sigma^T P \sigma} - 1 = 0
\end{equation}
\begin{equation}\label{eq:SFten}
SF_{ten} = \boldsymbol{|\frac{\sigma}{\sigma_{max}}|}
\end{equation}
\begin{equation}\label{eq:SFcomp}
SF_{comp} = \boldsymbol{|\frac{\sigma}{\sigma_{min}}|}
\end{equation}
Where \(\sigma_i\) is a single entry in the stress vector \(\boldsymbol{\sigma}\), \(SF_{ten}\) is the safety factor in tension and \(SF_{comp}\) is the safety factor in compression.
