deletions | additions       diff --git a/method.tex b/method.tex  index 8ad48af..ea0bed0 100644  --- a/method.tex  +++ b/method.tex 
...
\section{Method}  Stems were aproximated as truncated cone meshes, represented using quadratic tetrahedral elements. Each (fully grown) mesh was devided into elements covering one fifth the radious, one fifth the height and one sixteeth of the circunfrence. Data presented in Waghorn et al. (2007a) \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} 
...
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 acomidate these time dependent changes, the mesh needs to be 'grown', here time is treated in descrete steps. Wind effects are not acounted 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 longatudinly, for more details see ---refs---. \citep{archer_growth_1986}.  As the pre-existing strucute 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 utalises 1 year steps, a tree with wind imposed at 15 years has 3 year time steps. 
...
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 Pinus radiata usually follows a non-linear increase for density and a non-linear decrease for MFA (Burdon et al., 2004), \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 
...
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: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 Waghorn et al. (2007a). \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}.
...
\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 Beets and Whitehead (1996) \citet{Beets_1996}  at a stocking of 741 stems per hectare. The canopy force due to gravity is only applied to the subdomain \(\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 (Spatz and Bruechert, 2000; Rudnicki et al., 2004; Mayer et al., 1989) \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 (Coutts and Grace, 1995). \citep{coutts_wind_1995}.  The drag coefficient \(\varsigma\) in Equation \ref{eq:DC} was produced from data reported by Mayhead (1973), \citet{Mayhead_1973},  for Scotts pine as no data was available for radiata. The use of the Mayhead (1973) \citet{Mayhead_1973}  Scotts pine data set has previously been suggested as a suitable substitute (Moore and Gardiner, 2001). \citep{moore2001relative}.  \begin{equation}\label{eq:DC}  \varsigma = e^{-0.377\omega-0.306} 
...
\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 (Mayhead, 1973). \(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\). Papesch et al. (1997) \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  
...
Where \(M_c\) is the maximum bending moment, \(h\) is the height and \(D\) is deflection.  Growth stresses are represented in a simplified fashion. Gillis and Hsu (1979) \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 Timell (1986b) \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 Gillis and Hsu (1979) \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 
...
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 (Tsai and Wu, 1971) \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,as demonstrated in --ref--, for this calculation  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 --ref-- and the relationships presented in equation \ref{eq:constants} proportional limit surfaces for each point were calculated using Tsai and Wu’s \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